diff --git a/tnqvm/TNQVM.hpp b/tnqvm/TNQVM.hpp
index a6f41d8..18b962d 100644
--- a/tnqvm/TNQVM.hpp
+++ b/tnqvm/TNQVM.hpp
@@ -65,9 +65,16 @@ class TNQVM : public Accelerator {
       vqeMode = config.get<bool>("vqe-mode");
     }
 
-    if (config.stringExists("tnqvm-visitor")) {
-      const auto requestedBackend = config.getString("tnqvm-visitor");
-      const auto& allVisitorServices = xacc::getServices<TNQVMVisitor>();
+    if (config.stringExists("tnqvm-visitor") ||
+        config.stringExists("backend")) {
+      // Get the specific TNQVM visitor, either using the `tnqvm-visitor` key
+      // or the `backend` key.
+      // e.g., when users use the convention "tnqvm::exatn", "exatn" will be
+      // parsed and passed in the "backend" field.
+      const auto requestedBackend = config.stringExists("tnqvm-visitor")
+                                        ? config.getString("tnqvm-visitor")
+                                        : config.getString("backend");
+      const auto &allVisitorServices = xacc::getServices<TNQVMVisitor>();
       // We must have at least one TNQVM service registered.
       assert(!allVisitorServices.empty());
       bool foundRequestedBackend = false;
@@ -133,7 +140,13 @@ class TNQVM : public Accelerator {
   const std::string& getVisitorName() const { return backendName; }
   
   virtual ~TNQVM() {}
-
+  
+  virtual HeterogeneousMap getExecutionInfo() const override { 
+    auto result = visitor->getExecutionInfo();
+    result.insert("visitor", visitor->name());
+    return result; 
+  }
+  
   int verbose() const { return __verbose; }
   void verbose(int level) { __verbose = level; }
   void set_verbose(int level) { __verbose = level; }
diff --git a/tnqvm/visitors/CMakeLists.txt b/tnqvm/visitors/CMakeLists.txt
index 9d763f0..17f6c2f 100644
--- a/tnqvm/visitors/CMakeLists.txt
+++ b/tnqvm/visitors/CMakeLists.txt
@@ -32,4 +32,5 @@ include_directories(${CMAKE_CURRENT_SOURCE_DIR})
 add_subdirectory(itensor)
 add_subdirectory(exatn)
 add_subdirectory(exatn-mps)
-add_subdirectory(exatn-mpo)
\ No newline at end of file
+add_subdirectory(exatn-mpo)
+add_subdirectory(exatn-dm)
\ No newline at end of file
diff --git a/tnqvm/visitors/TNQVMVisitor.hpp b/tnqvm/visitors/TNQVMVisitor.hpp
index 9e40cca..1c50b6a 100644
--- a/tnqvm/visitors/TNQVMVisitor.hpp
+++ b/tnqvm/visitors/TNQVMVisitor.hpp
@@ -71,9 +71,14 @@ class TNQVMVisitor : public AllGateVisitor, public OptionsProvider,
   // Does this visitor implementation support VQE mode execution?
   // i.e. ability to cache the state vector after simulating the ansatz.
   virtual bool supportVqeMode() const { return false; }
+  // Execution information that visitor wants to persist.
+  HeterogeneousMap getExecutionInfo() const { return executionInfo; }
+
 protected:
   std::shared_ptr<AcceleratorBuffer> buffer;
   HeterogeneousMap options;
+  // Visitor impl to set if need be.
+  HeterogeneousMap executionInfo;
 };
 
 } // namespace tnqvm
diff --git a/tnqvm/visitors/exatn-dm/CMakeLists.txt b/tnqvm/visitors/exatn-dm/CMakeLists.txt
new file mode 100644
index 0000000..42be8f5
--- /dev/null
+++ b/tnqvm/visitors/exatn-dm/CMakeLists.txt
@@ -0,0 +1,79 @@
+find_package(ExaTN QUIET)
+
+if (ExaTN_FOUND)
+   message(STATUS "Found ExaTN at ${EXATN_ROOT}")
+   set (PACKAGE_NAME "TNQVM ExaTN Density Matrix Visitor")
+   set (PACKAGE_DESCIPTION "TNQVM ExaTN Density Matrix backend")
+   set (LIBRARY_NAME tnqvm-exatn-dm)
+
+  if (TNQVM_MPI_ENABLED) 
+   find_package(MPI REQUIRED)
+   message(STATUS "Found a suitable MPI compiler ${MPI_CXX_COMPILER}.")
+   message(STATUS "Compiler vendor is [${CMAKE_CXX_COMPILER_ID}]")
+   message(STATUS "Include path: ${MPI_CXX_INCLUDE_DIRS}")
+   message(STATUS "Compile flags: ${MPI_CXX_COMPILE_FLAGS}")
+   message(STATUS "Link flags: ${MPI_CXX_LINK_FLAGS}")
+   message(STATUS "Libraries: ${MPI_CXX_LIBRARIES}")
+   include_directories(${MPI_CXX_INCLUDE_DIRS})
+   link_libraries(${MPI_CXX_LIBRARIES})
+   add_definitions(-DTNQVM_MPI_ENABLED)
+  endif()  
+   
+   set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -DTNQVM_HAS_EXATN")
+   set(CMAKE_INSTALL_RPATH "${EXATN_ROOT}/lib")
+   set(EXATN_VISITOR_CPP_FILE ExaTnDmVisitor.cpp)
+
+   if (EXATN_BLAS_LIB MATCHES MKL)
+      # Fix for bug #30
+      message(STATUS "Exatn built with MKL, updating our build: ${EXATN_MKL_PATH}")
+      set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -DTNQVM_EXATN_USES_MKL_BLAS")
+      configure_file("${CMAKE_CURRENT_SOURCE_DIR}/ExaTnDmVisitor.cpp"
+               "${CMAKE_BINARY_DIR}/tnqvm/visitors/exatn-mpo/ExaTnDmVisitor.cpp" @ONLY)
+      set(EXATN_VISITOR_CPP_FILE ${CMAKE_BINARY_DIR}/tnqvm/visitors/exatn-mpo/ExaTnDmVisitor.cpp)
+   endif()
+
+   file (GLOB HEADERS *.hpp)
+   file (GLOB SRC ${EXATN_VISITOR_CPP_FILE})
+
+   usFunctionGetResourceSource(TARGET ${LIBRARY_NAME} OUT SRC)
+   usFunctionGenerateBundleInit(TARGET ${LIBRARY_NAME} OUT SRC)
+
+   add_library(${LIBRARY_NAME} SHARED ${SRC})
+
+   set(_bundle_name tnqvm_exatn_dm)
+   set_target_properties(${LIBRARY_NAME} PROPERTIES
+   # This is required for every bundle
+   COMPILE_DEFINITIONS US_BUNDLE_NAME=${_bundle_name}
+   # This is for convenience, used by other CMake functions
+   US_BUNDLE_NAME ${_bundle_name}
+   )
+
+   # Embed meta-data from a manifest.json file
+   usFunctionEmbedResources(TARGET ${LIBRARY_NAME}
+      WORKING_DIRECTORY
+      ${CMAKE_CURRENT_SOURCE_DIR}
+      FILES
+      manifest.json
+   )
+
+   target_include_directories(${LIBRARY_NAME} PUBLIC . .. ${XACC_DIR}/include/xacc/)
+
+   # Links to ExaTN using its linker config flags.
+   target_link_libraries(${LIBRARY_NAME} PUBLIC xacc::xacc exatn::exatn xacc::quantum_gate)
+
+   if(APPLE)
+      set_target_properties(${LIBRARY_NAME} PROPERTIES INSTALL_RPATH "@loader_path/../lib")
+      set_target_properties(${LIBRARY_NAME} PROPERTIES LINK_FLAGS "-undefined dynamic_lookup")
+   else()
+      set_target_properties(${LIBRARY_NAME} PROPERTIES INSTALL_RPATH "$ORIGIN/../lib;${EXATN_ROOT}/lib;${BLAS_PATH}")
+      set_target_properties(${LIBRARY_NAME} PROPERTIES LINK_FLAGS "-shared")
+   endif()
+
+   install(TARGETS ${LIBRARY_NAME} DESTINATION ${CMAKE_INSTALL_PREFIX}/plugins)
+
+   if(TNQVM_BUILD_TESTS)
+	   add_subdirectory(tests)
+   endif()
+else()
+   message(STATUS "ExaTN not found, skipping ExaTNVisitor build")
+endif()
\ No newline at end of file
diff --git a/tnqvm/visitors/exatn-dm/ExaTnDmVisitor.cpp b/tnqvm/visitors/exatn-dm/ExaTnDmVisitor.cpp
new file mode 100644
index 0000000..60ce905
--- /dev/null
+++ b/tnqvm/visitors/exatn-dm/ExaTnDmVisitor.cpp
@@ -0,0 +1,858 @@
+#include "ExaTnDmVisitor.hpp"
+#include "tensor_basic.hpp"
+#include "talshxx.hpp"
+#include "utils/GateMatrixAlgebra.hpp"
+#include "base/Gates.hpp"
+#include "NoiseModel.hpp"
+#include "xacc_service.hpp"
+#include "xacc_plugin.hpp"
+#ifdef TNQVM_EXATN_USES_MKL_BLAS
+#include <dlfcn.h>
+#endif
+#include <bits/stdc++.h> 
+
+#define QUBIT_DIM 2
+
+namespace {
+std::vector<std::complex<double>> Q_ZERO_TENSOR_BODY(size_t in_volume) {
+  std::vector<std::complex<double>> body(in_volume, {0.0, 0.0});
+  body[0] = std::complex<double>(1.0, 0.0);
+  return body;
+}
+
+const std::string ROOT_TENSOR_NAME = "Root";
+
+void printDensityMatrix(const exatn::TensorNetwork &in_tensorNet,
+                        size_t in_nbQubit, bool in_checkTrace = true) {
+  exatn::TensorNetwork tempNetwork(in_tensorNet);
+  tempNetwork.rename("__TEMP__" + in_tensorNet.getName());
+  const bool evaledOk = exatn::evaluateSync(tempNetwork);
+  assert(evaledOk);
+  auto talsh_tensor =
+      exatn::getLocalTensor(tempNetwork.getTensor(0)->getName());
+  const auto expectedDensityMatrixVolume =
+      (1ULL << in_nbQubit) * (1ULL << in_nbQubit);
+  const auto nbRows = 1ULL << in_nbQubit;
+  if (talsh_tensor) {
+    const std::complex<double> *body_ptr;
+    const bool access_granted = talsh_tensor->getDataAccessHostConst(&body_ptr);
+    if (!access_granted) {
+      std::cout << "Failed to retrieve tensor data!!!\n";
+    } else {
+      assert(expectedDensityMatrixVolume == talsh_tensor->getVolume());
+      std::complex<double> traceVal{0.0, 0.0};
+      int rowId = -1;
+      for (int i = 0; i < talsh_tensor->getVolume(); ++i) {
+        if (i % nbRows == 0) {
+          std::cout << "\n";
+          rowId++;
+        }
+        if (i == (rowId + rowId * nbRows)) {
+          traceVal += body_ptr[i];
+        }
+        const auto &elem = body_ptr[i];
+        std::cout << elem << " ";
+      }
+      std::cout << "\n";
+      if (in_checkTrace) {
+        // Verify trace(density matrix) == 1.0
+        assert(std::abs(traceVal.real() - 1.0) < 1e-12);
+        assert(std::abs(traceVal.imag()) < 1e-12);
+      }
+    }
+  } else {
+    std::cout << "Failed to retrieve tensor data!!!\n";
+  }
+}
+
+bool checkDmTensorNetwork(const exatn::TensorNetwork &in_tensorNet,
+                          size_t in_nbQubit) {
+  exatn::TensorNetwork tempNetwork(in_tensorNet);
+  tempNetwork.rename("__TEMP__" + in_tensorNet.getName());
+  const bool evaledOk = exatn::evaluateSync(tempNetwork);
+  assert(evaledOk);
+  auto talsh_tensor =
+      exatn::getLocalTensor(tempNetwork.getTensor(0)->getName());
+  const auto expectedDensityMatrixVolume =
+      (1ULL << in_nbQubit) * (1ULL << in_nbQubit);
+  const auto nbRows = 1ULL << in_nbQubit;
+  if (talsh_tensor) {
+    const std::complex<double> *body_ptr;
+    const bool access_granted = talsh_tensor->getDataAccessHostConst(&body_ptr);
+    if (!access_granted) {
+      std::cout << "Failed to retrieve tensor data!!!\n";
+    } else {
+      assert(expectedDensityMatrixVolume == talsh_tensor->getVolume());
+      std::complex<double> traceVal{0.0, 0.0};
+      int rowId = -1;
+      for (int i = 0; i < talsh_tensor->getVolume(); ++i) {
+        if (i % nbRows == 0) {
+          rowId++;
+        }
+        if (i == (rowId + rowId * nbRows)) {
+          traceVal += body_ptr[i];
+        }
+        const auto &elem = body_ptr[i];
+      }
+      return (std::abs(traceVal.real() - 1.0) < 1e-12) &&
+             (std::abs(traceVal.imag()) < 1e-12);
+    }
+  } else {
+    std::cout << "Failed to retrieve tensor data!!!\n";
+  }
+  return false;
+}
+
+std::vector<std::complex<double>> flattenGateMatrix(
+    const std::vector<std::vector<std::complex<double>>> &in_gateMatrix) {
+  std::vector<std::complex<double>> resultVector;
+  resultVector.reserve(in_gateMatrix.size() * in_gateMatrix.size());
+  for (const auto &row : in_gateMatrix) {
+    for (const auto &entry : row) {
+      resultVector.emplace_back(entry);
+    }
+  }
+
+  return resultVector;
+}
+
+std::vector<std::vector<std::complex<double>>> conjugateMatrix(
+    const std::vector<std::vector<std::complex<double>>> &in_mat) {
+  auto result = in_mat;
+  const auto dim = in_mat.size();
+  for (size_t rowId = 0; rowId < dim; ++rowId) {
+    for (size_t colId = 0; colId < dim; ++colId) {
+      result[colId][rowId] = std::conj(in_mat[colId][rowId]);
+    }
+  }
+  return result;
+}
+
+std::vector<std::complex<double>>
+getGateMatrix(const xacc::Instruction &in_gate, bool in_dagger = false) {
+  using namespace tnqvm;
+  const auto gateEnum = GetGateType(in_gate.name());
+  const auto getMatrix = [&]() {
+    switch (gateEnum) {
+    case CommonGates::Rx:
+      return GetGateMatrix<CommonGates::Rx>(
+          in_gate.getParameter(0).as<double>());
+    case CommonGates::Ry:
+      return GetGateMatrix<CommonGates::Ry>(
+          in_gate.getParameter(0).as<double>());
+    case CommonGates::Rz:
+      return GetGateMatrix<CommonGates::Rz>(
+          in_gate.getParameter(0).as<double>());
+    case CommonGates::I:
+      return GetGateMatrix<CommonGates::I>();
+    case CommonGates::H:
+      return GetGateMatrix<CommonGates::H>();
+    case CommonGates::X:
+      return GetGateMatrix<CommonGates::X>();
+    case CommonGates::Y:
+      return GetGateMatrix<CommonGates::Y>();
+    case CommonGates::Z:
+      return GetGateMatrix<CommonGates::Z>();
+    case CommonGates::T:
+      return GetGateMatrix<CommonGates::T>();
+    case CommonGates::Tdg:
+      return GetGateMatrix<CommonGates::Tdg>();
+    case CommonGates::CNOT:
+      return GetGateMatrix<CommonGates::CNOT>();
+    case CommonGates::Swap:
+      return GetGateMatrix<CommonGates::Swap>();
+    case CommonGates::iSwap:
+      return GetGateMatrix<CommonGates::iSwap>();
+    case CommonGates::fSim:
+      return GetGateMatrix<CommonGates::fSim>(
+          in_gate.getParameter(0).as<double>(),
+          in_gate.getParameter(1).as<double>());
+    default:
+      return GetGateMatrix<CommonGates::I>();
+    }
+  };  
+
+  const auto gateMatrix = getMatrix();
+  return in_dagger ? flattenGateMatrix(conjugateMatrix(gateMatrix))
+                   : flattenGateMatrix(gateMatrix);
+}
+
+void recursiveFindAllCombinations(std::vector<std::vector<unsigned int>>& io_result,
+const std::vector<unsigned int> &arr,
+                std::vector<unsigned int> &data, int start, int end, int index,
+                int r) {
+  // Current combination is ready
+  if (index == r) {
+    std::vector<unsigned int> combination;
+    for (int j = 0; j < r; j++) {
+      combination.emplace_back(data[j]);
+    }
+    io_result.emplace_back(combination);
+    return;
+  }
+  for (int i = start; i <= end && end - i + 1 >= r - index; i++) {
+    data[index] = arr[i];
+    recursiveFindAllCombinations(io_result, arr, data, i + 1, end, index + 1, r);
+  }
+}
+
+std::vector<std::vector<unsigned int>>
+findAllPermutations(const std::vector<unsigned int> &in_elements) {
+  auto a = in_elements;
+  std::sort(a.begin(), a.end());
+  std::vector<std::vector<unsigned int>> result;
+  do {
+    result.emplace_back(a);
+
+  } while (std::next_permutation(a.begin(), a.end()));
+  return result;
+}
+
+std::vector<unsigned int>
+filterList(const std::vector<unsigned int> &in_totalList,
+           const std::vector<unsigned int> &in_filter) {
+  std::vector<unsigned int> result;
+  for (const auto &el : in_totalList) {
+    if (!xacc::container::contains(in_filter, el)) {
+      result.emplace_back(el);
+    }
+  }
+  return result;
+}
+
+// Helper to find all input->output leg connection permutation.
+// This is *only* used for development purposes to figure out the
+// correct connection configs to preseve the order of the root density matrix
+// tensor.
+std::vector<std::vector<std::pair<unsigned int, unsigned int>>>
+findAllGateTensorPairings(int in_tensorRank) {
+  assert(in_tensorRank % 2 == 0);
+  std::vector<unsigned int> legIds;
+  std::vector<unsigned int> tempVec(in_tensorRank / 2);
+  for (unsigned int i = 0; i < in_tensorRank; ++i) {
+    legIds.emplace_back(i);
+  }
+  // Find all possible input leg combination:
+  std::vector<std::vector<unsigned int>> inputCombinations;
+  recursiveFindAllCombinations(inputCombinations, legIds, tempVec, 0,
+                               in_tensorRank - 1, 0, in_tensorRank / 2);
+
+  std::vector<std::vector<std::pair<unsigned int, unsigned int>>> result;
+  for (const auto &inputLegs : inputCombinations) {
+    const auto inputLegPerm = findAllPermutations(inputLegs);
+    const auto outputLegPerm =
+        findAllPermutations(filterList(legIds, inputLegs));
+    for (const auto &inputPerm : inputLegPerm) {
+      for (const auto &outputPerm : outputLegPerm) {
+        assert(inputPerm.size() == outputPerm.size());
+        std::vector<std::pair<unsigned int, unsigned int>> combo;
+        for (size_t i = 0; i < inputPerm.size(); ++i) {
+          combo.emplace_back(std::make_pair(inputPerm[i], outputPerm[i]));
+        }
+        result.emplace_back(combo);
+      }
+    }
+  }
+
+  return result;
+}
+
+void checkKrausTensorConfig(const exatn::TensorNetwork &in_tensorNet,
+                            unsigned int in_q0, unsigned int in_q1,
+                            unsigned int in_numQubits,
+                            const std::string &in_krausTensorName) {
+  auto noiseTensor = exatn::getTensor(in_krausTensorName);
+  std::cout << "Checking Kraus tensor connections:\n";
+  auto test = findAllGateTensorPairings(noiseTensor->getDimExtents().size());
+  for (const auto &config : test) {
+    const std::vector<
+        std::pair<unsigned int, std::pair<unsigned int, unsigned int>>>
+        noisePairing{{in_q0, config[0]},
+                     {in_q0 + in_numQubits, config[1]},
+                     {in_q1, config[2]},
+                     {in_q1 + in_numQubits, config[3]}};
+
+    // Append the tensor for this gate to the network
+    auto temp = in_tensorNet;
+    const auto tensorIdCounter = temp.getMaxTensorId() + 1;
+    const bool appended = temp.appendTensorGateGeneral(
+        tensorIdCounter, noiseTensor, noisePairing);
+    assert(appended);
+    if (checkDmTensorNetwork(temp, in_numQubits)) {
+      std::cout << "============================== \n Valid: ";
+      for (const auto &[iLeg, oLeg] : config) {
+        std::cout << "(" << iLeg << "-->" << oLeg << ")";
+      }
+      printDensityMatrix(temp, in_numQubits);
+      std::cout << "==============================\n";
+    }
+  }
+}
+} // namespace
+
+namespace tnqvm {
+ExaTnDmVisitor::ExaTnDmVisitor() {
+  // TODO
+}
+
+void ExaTnDmVisitor::initialize(std::shared_ptr<AcceleratorBuffer> buffer,
+                                int nbShots) {
+  // Initialize ExaTN (if not already initialized)
+  if (!exatn::isInitialized()) {
+#ifdef TNQVM_EXATN_USES_MKL_BLAS
+    // Fix for TNQVM bug #30
+    void *core_handle =
+        dlopen("@EXATN_MKL_PATH@/libmkl_core@CMAKE_SHARED_LIBRARY_SUFFIX@",
+               RTLD_LAZY | RTLD_GLOBAL);
+    if (core_handle == nullptr) {
+      std::string err = std::string(dlerror());
+      xacc::error("Could not load mkl_core - " + err);
+    }
+
+    void *thread_handle = dlopen(
+        "@EXATN_MKL_PATH@/libmkl_gnu_thread@CMAKE_SHARED_LIBRARY_SUFFIX@",
+        RTLD_LAZY | RTLD_GLOBAL);
+    if (thread_handle == nullptr) {
+      std::string err = std::string(dlerror());
+      xacc::error("Could not load mkl_gnu_thread - " + err);
+    }
+#endif
+    // If exaTN has not been initialized, do it now.
+    exatn::initialize();
+    // ExaTN and XACC logging levels are always in-synced.
+    exatn::resetRuntimeLoggingLevel(xacc::verbose ? xacc::getLoggingLevel()
+                                                  : 0);
+    xacc::subscribeLoggingLevel([](int level) {
+      exatn::resetRuntimeLoggingLevel(xacc::verbose ? level : 0);
+    });
+  }
+  m_buffer = buffer;
+  m_tensorNetwork = buildInitialNetwork(buffer->size());
+  m_tensorIdCounter = m_tensorNetwork.getMaxTensorId();
+  if (options.pointerLikeExists<xacc::NoiseModel>("noise-model")) {
+    m_noiseConfig = xacc::as_shared_ptr(
+        options.getPointerLike<xacc::NoiseModel>("noise-model"));
+  }
+}
+
+exatn::TensorNetwork
+ExaTnDmVisitor::buildInitialNetwork(size_t in_nbQubits) const {
+  for (int i = 0; i < in_nbQubits; ++i) {
+    const std::string tensorName = "Q" + std::to_string(i);
+    auto tensor = std::make_shared<exatn::Tensor>(
+            tensorName, exatn::TensorShape{QUBIT_DIM, QUBIT_DIM});
+    const bool created =
+        exatn::createTensorSync(tensor, exatn::TensorElementType::COMPLEX64);
+    assert(created);
+    const bool initialized = exatn::initTensorDataSync(
+        tensorName, Q_ZERO_TENSOR_BODY(tensor->getVolume()));
+    assert(initialized);
+  }
+  exatn::TensorNetwork qubitTensorNet("Tensor_Network");
+  // Append the qubit tensors to the tensor network
+  size_t tensorIdCounter = 0;
+  for (int i = 0; i < in_nbQubits; ++i) {
+    tensorIdCounter++;
+    const std::string tensorName = "Q" + std::to_string(i);
+    qubitTensorNet.appendTensor(
+        tensorIdCounter, exatn::getTensor(tensorName),
+        std::vector<std::pair<unsigned int, unsigned int>>{});
+  }
+  // std::cout << "Qubit Tensor network:\n";
+  // qubitTensorNet.printIt();
+  return qubitTensorNet;
+}
+
+void ExaTnDmVisitor::finalize() {
+  executionInfo.clear();
+  // Max number of qubits that we allow for a full density matrix retrieval.
+  // For more qubits, only expectation contraction is supported.
+  constexpr size_t MAX_SIZE_TO_COLLAPSE_DM = 10; 
+  
+  if (m_buffer->size() <= MAX_SIZE_TO_COLLAPSE_DM) {
+    xacc::ExecutionInfo::DensityMatrixType densityMatrix;
+    const auto nbRows = 1ULL << m_buffer->size();
+    densityMatrix.reserve(1 << m_buffer->size());
+    exatn::TensorNetwork tempNetwork(m_tensorNetwork);
+    tempNetwork.rename("__TEMP__" + m_tensorNetwork.getName());
+    const bool evaledOk = exatn::evaluateSync(tempNetwork);
+    assert(evaledOk);
+    auto talsh_tensor =
+        exatn::getLocalTensor(tempNetwork.getTensor(0)->getName());
+    const auto expectedDensityMatrixVolume = nbRows * nbRows;
+    if (talsh_tensor) {
+      const std::complex<double> *body_ptr;
+      const bool access_granted =
+          talsh_tensor->getDataAccessHostConst(&body_ptr);
+      if (!access_granted) {
+        std::cout << "Failed to retrieve tensor data!!!\n";
+      } else {
+        assert(expectedDensityMatrixVolume == talsh_tensor->getVolume());
+        xacc::ExecutionInfo::WaveFuncType dmRow;
+        dmRow.reserve(nbRows);
+        for (int i = 0; i < talsh_tensor->getVolume(); ++i) {
+          if ((i != 0) && (i % nbRows == 0)) {
+            assert(dmRow.size() == nbRows);
+            densityMatrix.emplace_back(std::move(dmRow));
+            dmRow.clear();
+            dmRow.reserve(nbRows);
+          }
+          const auto &elem = body_ptr[i];
+          dmRow.emplace_back(elem);
+        }
+        assert(dmRow.size() == nbRows);
+        densityMatrix.emplace_back(std::move(dmRow));
+        assert(densityMatrix.size() == nbRows);
+      }
+    }
+    
+    executionInfo.insert(ExecutionInfo::DmKey, std::make_shared<ExecutionInfo::DensityMatrixType>(std::move(densityMatrix)));
+  }
+
+  // Expectation value calculation:
+  // Exp = Trace(rho * Op)
+  if (!m_measuredBits.empty()) {
+    auto tensorIdCounter = m_tensorIdCounter;
+    auto expValTensorNet = m_tensorNetwork;
+    const std::string measZTensorName = "MEAS_Z";
+    {
+      xacc::quantum::Z zGate(0);
+      const auto gateMatrix = getGateMatrix(zGate);
+      assert(gateMatrix.size() == 4);
+      // Create the tensor
+      const bool created = exatn::createTensorSync(
+          measZTensorName, exatn::TensorElementType::COMPLEX64,
+          exatn::TensorShape{2, 2});
+      assert(created);
+      // Init tensor body data
+      const bool initialized =
+          exatn::initTensorDataSync(measZTensorName, gateMatrix);
+      assert(initialized);
+    }
+    // Add Z tensors for measurement
+    for (const auto &measBit : m_measuredBits) {
+      tensorIdCounter++;
+      const std::vector<unsigned int> gatePairing{
+          static_cast<unsigned int>(measBit)};
+      // Append the tensor for this gate to the network
+      const bool appended = expValTensorNet.appendTensorGate(
+          tensorIdCounter,
+          // Get the gate tensor data which must have been initialized.
+          exatn::getTensor(measZTensorName),
+          // which qubits that the gate is acting on
+          gatePairing);
+      assert(appended);
+    }
+    // DEBUG:
+    // std::cout << "TENSOR NETWORK TO COMPUTE THE TRACE:\n";
+    // printDensityMatrix(expValTensorNet, m_buffer->size(), false);
+    // Compute the trace, closing the tensor network:
+    const std::string idTensor = "ID_TRACE";
+    {
+      xacc::quantum::Identity idGate(0);
+      const auto idGateMatrix = getGateMatrix(idGate);
+      // Create the tensor
+      const bool created =
+          exatn::createTensorSync(idTensor, exatn::TensorElementType::COMPLEX64,
+                                  exatn::TensorShape{2, 2});
+      assert(created);
+      // Init tensor body data
+      const bool initialized =
+          exatn::initTensorDataSync(idTensor, idGateMatrix);
+      assert(initialized);
+    }
+
+    for (size_t qId = 0; qId < m_buffer->size(); ++qId) {
+      tensorIdCounter++;
+      const std::vector<std::pair<unsigned int, unsigned int>> tracePairing{
+          {static_cast<unsigned int>(0), 0},
+          {static_cast<unsigned int>(m_buffer->size() - qId), 1}};
+      const bool appended = expValTensorNet.appendTensor(
+          tensorIdCounter, exatn::getTensor(idTensor), tracePairing);
+      assert(appended);
+    }
+
+    // Evaluate the trace by contraction:
+    const bool evaledOk = exatn::evaluateSync(expValTensorNet);
+    assert(evaledOk);
+    auto talsh_tensor =
+        exatn::getLocalTensor(expValTensorNet.getTensor(0)->getName());
+    if (talsh_tensor) {
+      const std::complex<double> *body_ptr;
+      const bool access_granted =
+          talsh_tensor->getDataAccessHostConst(&body_ptr);
+      if (!access_granted) {
+        xacc::error("Failed to retrieve tensor data!");
+      } else {
+        // Should only have 1 element:
+        assert(talsh_tensor->getVolume() == 1);
+        const std::complex<double> traceVal = *body_ptr;
+        assert(traceVal.imag() < 1e-3);
+        const double expValZ = traceVal.real();
+        // std::cout << "Exp-val Z = " << expValZ << "\n";
+        m_buffer->addExtraInfo("exp-val-z", expValZ);
+      }
+    }
+    bool destroyed = exatn::destroyTensor(measZTensorName);
+    assert(destroyed);
+    destroyed = exatn::destroyTensor(idTensor);
+    assert(destroyed);
+  }
+
+  std::unordered_set<std::string> tensorList;
+  for (auto iter = m_tensorNetwork.cbegin(); iter != m_tensorNetwork.cend();
+       ++iter) {
+    const auto &tensorName = iter->second.getTensor()->getName();
+    // Not a root tensor
+    if (!tensorName.empty() && tensorName[0] != '_') {
+      tensorList.emplace(iter->second.getTensor()->getName());
+    }
+  }
+
+  for (const auto &tensorName : tensorList) {
+    const bool destroyed = exatn::destroyTensor(tensorName);
+    assert(destroyed);
+  }
+
+  m_buffer.reset();
+  m_noiseConfig.reset();
+}
+
+void ExaTnDmVisitor::applySingleQubitGate(
+    xacc::quantum::Gate &in_gateInstruction) {
+  {
+    m_tensorIdCounter++;
+    assert(in_gateInstruction.bits().size() == 1);
+    const auto gateMatrix = getGateMatrix(in_gateInstruction);
+    // std::cout << "Gate mattrix:\n";
+    // for (const auto &el : gateMatrix) {
+    //   std::cout << el << "\n";
+    // }
+    assert(gateMatrix.size() == 4);
+    const std::string uniqueGateName =
+        in_gateInstruction.name() + std::to_string(m_tensorIdCounter);
+    // Create the tensor
+    const bool created = exatn::createTensorSync(
+        uniqueGateName, exatn::TensorElementType::COMPLEX64,
+        exatn::TensorShape{2, 2});
+    assert(created);
+    // Init tensor body data
+    const bool initialized =
+        exatn::initTensorDataSync(uniqueGateName, gateMatrix);
+    assert(initialized);
+
+    const std::vector<unsigned int> gatePairing{
+        static_cast<unsigned int>(in_gateInstruction.bits()[0])};
+    // Append the tensor for this gate to the network
+    const bool appended = m_tensorNetwork.appendTensorGate(
+        m_tensorIdCounter,
+        // Get the gate tensor data which must have been initialized.
+        exatn::getTensor(uniqueGateName),
+        // which qubits that the gate is acting on
+        gatePairing);
+    assert(appended);
+  }
+
+  {
+    // Append the dagger gate
+    m_tensorIdCounter++;
+    const auto gateMatrix = getGateMatrix(in_gateInstruction, true);
+    // std::cout << "Gate mattrix dagger:\n";
+    // for (const auto &el : gateMatrix) {
+    //   std::cout << el << "\n";
+    // }
+    assert(gateMatrix.size() == 4);
+    const std::string uniqueGateName =
+        in_gateInstruction.name() + "_CONJ_" + std::to_string(m_tensorIdCounter);
+    // Create the tensor
+    const bool created = exatn::createTensorSync(
+        uniqueGateName, exatn::TensorElementType::COMPLEX64,
+        exatn::TensorShape{2, 2});
+    assert(created);
+    // Init tensor body data
+    const bool initialized =
+        exatn::initTensorDataSync(uniqueGateName, gateMatrix);
+    assert(initialized);
+    const std::vector<unsigned int> gatePairingConj{static_cast<unsigned int>(
+        m_buffer->size() + in_gateInstruction.bits()[0])};
+    const bool conjAppended = m_tensorNetwork.appendTensorGate(
+        m_tensorIdCounter,
+        // Get the gate tensor data which must have been initialized.
+        exatn::getTensor(uniqueGateName),
+        // which qubits that the gate is acting on
+        gatePairingConj);
+    assert(conjAppended);
+  }
+
+  // DEBUG
+  // std::cout << "Before noise:\n";
+  // m_tensorNetwork.printIt();
+
+  // Adding noise tensors.
+  applyNoise(in_gateInstruction);
+}
+
+void ExaTnDmVisitor::applyTwoQubitGate(
+    xacc::quantum::Gate &in_gateInstruction) {
+  {
+    m_tensorIdCounter++;
+    assert(in_gateInstruction.bits().size() == 2);
+    const auto gateMatrix = getGateMatrix(in_gateInstruction);
+    assert(gateMatrix.size() == 16);
+    const std::string uniqueGateName =
+        in_gateInstruction.name() + std::to_string(m_tensorIdCounter);
+    // Create the tensor
+    const bool created = exatn::createTensorSync(
+        uniqueGateName, exatn::TensorElementType::COMPLEX64,
+        exatn::TensorShape{2, 2, 2, 2});
+    assert(created);
+    // Init tensor body data
+    const bool initialized =
+        exatn::initTensorDataSync(uniqueGateName, gateMatrix);
+    assert(initialized);
+
+    const std::vector<unsigned int> gatePairing{
+        static_cast<unsigned int>(in_gateInstruction.bits()[1]),
+        static_cast<unsigned int>(in_gateInstruction.bits()[0])};
+    // Append the tensor for this gate to the network
+    const bool appended = m_tensorNetwork.appendTensorGate(
+        m_tensorIdCounter,
+        // Get the gate tensor data which must have been initialized.
+        exatn::getTensor(uniqueGateName),
+        // which qubits that the gate is acting on
+        gatePairing);
+    assert(appended);
+  }
+  {
+    m_tensorIdCounter++;
+    const auto gateMatrix = getGateMatrix(in_gateInstruction, true);
+    assert(gateMatrix.size() == 16);
+    const std::string uniqueGateName =
+        in_gateInstruction.name() + "_CONJ_" + std::to_string(m_tensorIdCounter);
+    // Create the tensor
+    const bool created = exatn::createTensorSync(
+        uniqueGateName, exatn::TensorElementType::COMPLEX64,
+        exatn::TensorShape{2, 2, 2, 2});
+    assert(created);
+    // Init tensor body data
+    const bool initialized =
+        exatn::initTensorDataSync(uniqueGateName, gateMatrix);
+    assert(initialized);
+    const std::vector<unsigned int> gatePairingConj{
+        static_cast<unsigned int>(m_buffer->size() +
+                                  in_gateInstruction.bits()[1]),
+        static_cast<unsigned int>(m_buffer->size() +
+                                  in_gateInstruction.bits()[0])};
+    const bool conjAppended = m_tensorNetwork.appendTensorGate(
+        m_tensorIdCounter,
+        // Get the gate tensor data which must have been initialized.
+        exatn::getTensor(uniqueGateName),
+        // which qubits that the gate is acting on
+        gatePairingConj);
+    assert(conjAppended);
+  }
+
+  // Adding noise tensors.
+  applyNoise(in_gateInstruction);
+}
+
+void ExaTnDmVisitor::applyNoise(xacc::quantum::Gate &in_gateInstruction) {
+  if (!m_noiseConfig) {
+    return;
+  }
+  // std::cout << "Before noise:\n";
+  // printDensityMatrix(m_tensorNetwork, m_buffer->size());
+
+  const auto noiseChannels =
+      m_noiseConfig->getNoiseChannels(in_gateInstruction);
+  auto noiseUtils = xacc::getService<NoiseModelUtils>("default");
+
+  for (auto &channel : noiseChannels) {
+    auto noiseMat = noiseUtils->krausToChoi(channel.mats);
+
+    // for (const auto &row : noiseMat) {
+    //   for (const auto &el : row) {
+    //     std::cout << el << " ";
+    //   }
+    //   std::cout << "\n";
+    // }
+
+    m_tensorIdCounter++;
+    const std::string noiseTensorName = in_gateInstruction.name() + "_Noise_" +
+                                        std::to_string(m_tensorIdCounter);
+    if (channel.noise_qubits.size() == 1) {
+      // Create the tensor
+      const bool created = exatn::createTensorSync(
+          noiseTensorName, exatn::TensorElementType::COMPLEX64,
+          exatn::TensorShape{2, 2, 2, 2});
+      assert(created);
+      // Init tensor body data
+      const bool initialized = exatn::initTensorDataSync(
+          noiseTensorName, flattenGateMatrix(noiseMat));
+      assert(initialized);
+
+      const std::vector<
+          std::pair<unsigned int, std::pair<unsigned int, unsigned int>>>
+          noisePairing{
+              {static_cast<unsigned int>(channel.noise_qubits[0]), {1, 0}},
+              {static_cast<unsigned int>(channel.noise_qubits[0] +
+                                         m_buffer->size()),
+               {3, 2}}};
+      // Append the tensor for this gate to the network
+      const bool appended = m_tensorNetwork.appendTensorGateGeneral(
+          m_tensorIdCounter, exatn::getTensor(noiseTensorName),
+          // which qubits that the gate is acting on
+          noisePairing);
+      assert(appended);
+    } else if (channel.noise_qubits.size() == 2) {
+      // Create the tensor
+      const bool created = exatn::createTensorSync(
+          noiseTensorName, exatn::TensorElementType::COMPLEX64,
+          exatn::TensorShape{2, 2, 2, 2, 2, 2, 2, 2});
+      assert(created);
+      // Init tensor body data
+      const bool initialized = exatn::initTensorDataSync(
+          noiseTensorName, flattenGateMatrix(noiseMat));
+      assert(initialized);
+
+      // DEBUG: Using this to check all possible permutations.
+      // checkKrausTensorConfig(m_tensorNetwork, channel.noise_qubits[0],
+      //                        channel.noise_qubits[1], m_buffer->size(),
+      //                        noiseTensorName);
+
+      const std::vector<
+          std::pair<unsigned int, std::pair<unsigned int, unsigned int>>>
+          noisePairingMsb{
+              {static_cast<unsigned int>(channel.noise_qubits[0]), {5, 6}},
+              {static_cast<unsigned int>(channel.noise_qubits[0] +
+                                         m_buffer->size()),
+               {1, 2}},
+              {static_cast<unsigned int>(channel.noise_qubits[1]), {4, 7}},
+              {static_cast<unsigned int>(channel.noise_qubits[1] +
+                                         m_buffer->size()),
+               {0, 3}}};
+
+      const std::vector<
+          std::pair<unsigned int, std::pair<unsigned int, unsigned int>>>
+          noisePairingLsb{
+              {static_cast<unsigned int>(channel.noise_qubits[1]), {5, 6}},
+              {static_cast<unsigned int>(channel.noise_qubits[1] +
+                                         m_buffer->size()),
+               {1, 2}},
+              {static_cast<unsigned int>(channel.noise_qubits[0]), {4, 7}},
+              {static_cast<unsigned int>(channel.noise_qubits[0] +
+                                         m_buffer->size()),
+               {0, 3}}};
+      const auto noisePairing = (channel.bit_order == KrausMatBitOrder::MSB)
+                                    ? noisePairingMsb
+                                    : noisePairingLsb;
+      // Append the tensor for this gate to the network
+      const bool appended = m_tensorNetwork.appendTensorGateGeneral(
+          m_tensorIdCounter, exatn::getTensor(noiseTensorName),
+          // which qubits that the gate is acting on
+          noisePairing);
+      assert(appended);
+    } else {
+      xacc::error("Unsupported noise data.");
+    }
+  }
+
+  // DEBUG:
+  // std::cout << "Apply Noise: " << in_gateInstruction.toString() << "\n";
+  // m_tensorNetwork.printIt();
+  // printDensityMatrix(m_tensorNetwork, m_buffer->size());
+}
+
+void ExaTnDmVisitor::visit(Identity &in_IdentityGate) {
+  applySingleQubitGate(in_IdentityGate);
+}
+
+void ExaTnDmVisitor::visit(Hadamard &in_HadamardGate) {
+  applySingleQubitGate(in_HadamardGate);
+  // DEBUG:
+  // std::cout << "Apply: " << in_HadamardGate.toString() << "\n";
+  // printDensityMatrix(m_tensorNetwork, m_buffer->size());
+}
+
+void ExaTnDmVisitor::visit(X &in_XGate) {
+  applySingleQubitGate(in_XGate);
+  // DEBUG:
+  // std::cout << "Apply: " << in_XGate.toString() << "\n";
+  // m_tensorNetwork.printIt();
+  // printDensityMatrix(m_tensorNetwork, m_buffer->size());
+}
+
+void ExaTnDmVisitor::visit(Y &in_YGate) { applySingleQubitGate(in_YGate); }
+
+void ExaTnDmVisitor::visit(Z &in_ZGate) { applySingleQubitGate(in_ZGate); }
+
+void ExaTnDmVisitor::visit(Rx &in_RxGate) { 
+  applySingleQubitGate(in_RxGate); 
+  // std::cout << "Apply: " << in_RxGate.toString() << "\n";
+  // m_tensorNetwork.printIt();
+  // printDensityMatrix(m_tensorNetwork, m_buffer->size());
+}
+
+void ExaTnDmVisitor::visit(Ry &in_RyGate) { 
+  applySingleQubitGate(in_RyGate); 
+  // std::cout << "Apply: " << in_RyGate.toString() << "\n";
+  // m_tensorNetwork.printIt();
+  // printDensityMatrix(m_tensorNetwork, m_buffer->size());
+}
+
+void ExaTnDmVisitor::visit(Rz &in_RzGate) { applySingleQubitGate(in_RzGate); }
+
+void ExaTnDmVisitor::visit(T &in_TGate) { 
+  applySingleQubitGate(in_TGate); 
+  // DEBUG:
+  // std::cout << "Apply: " << in_TGate.toString() << "\n";
+  // m_tensorNetwork.printIt();
+  // printDensityMatrix(m_tensorNetwork, m_buffer->size());
+}
+
+void ExaTnDmVisitor::visit(Tdg &in_TdgGate) {
+  applySingleQubitGate(in_TdgGate);
+}
+
+// others
+void ExaTnDmVisitor::visit(Measure &in_MeasureGate) {
+  m_measuredBits.emplace_back(in_MeasureGate.bits()[0]);
+}
+
+void ExaTnDmVisitor::visit(U &in_UGate) { applySingleQubitGate(in_UGate); }
+
+// two-qubit gates:
+void ExaTnDmVisitor::visit(CNOT &in_CNOTGate) {
+  applyTwoQubitGate(in_CNOTGate);
+  // DEBUG:
+  // std::cout << "Apply: " << in_CNOTGate.toString() << "\n";
+  // m_tensorNetwork.printIt();
+  // printDensityMatrix(m_tensorNetwork, m_buffer->size());
+}
+
+void ExaTnDmVisitor::visit(Swap &in_SwapGate) {
+  applyTwoQubitGate(in_SwapGate);
+}
+
+void ExaTnDmVisitor::visit(CZ &in_CZGate) { applyTwoQubitGate(in_CZGate); }
+
+void ExaTnDmVisitor::visit(CPhase &in_CPhaseGate) {
+  applyTwoQubitGate(in_CPhaseGate);
+}
+
+void ExaTnDmVisitor::visit(iSwap &in_iSwapGate) {
+  applyTwoQubitGate(in_iSwapGate);
+}
+
+void ExaTnDmVisitor::visit(fSim &in_fsimGate) {
+  applyTwoQubitGate(in_fsimGate);
+}
+
+const double ExaTnDmVisitor::getExpectationValueZ(
+    std::shared_ptr<CompositeInstruction> in_function) {
+  return 0.0;
+}
+} // namespace tnqvm
+
+// Register with CppMicroservices
+REGISTER_PLUGIN(tnqvm::ExaTnDmVisitor, tnqvm::TNQVMVisitor)
diff --git a/tnqvm/visitors/exatn-dm/ExaTnDmVisitor.hpp b/tnqvm/visitors/exatn-dm/ExaTnDmVisitor.hpp
new file mode 100644
index 0000000..e8f1cda
--- /dev/null
+++ b/tnqvm/visitors/exatn-dm/ExaTnDmVisitor.hpp
@@ -0,0 +1,99 @@
+#pragma once
+
+#include "TNQVMVisitor.hpp"
+#include "tensor_network.hpp"
+#include "exatn.hpp"
+
+// Full density matrix noisy visitor:
+// Name: "exatn-dm"
+// Supported initialization keys:
+// +-----------------------------+------------------------------------------------------------------------+-------------+--------------------------+
+// |  Initialization Parameter   |                  Parameter Description                                 |    type     |         default          |
+// +=============================+========================================================================+=============+==========================+
+// | backend-json                | Backend configuration JSON to estimate the noise model from.           |    string   | None                     |
+// +-----------------------------+------------------------------------------------------------------------+-------------+--------------------------+
+// | backend                     | Name of the IBMQ backend to query the backend configuration.           |    string   | None                     |
+// +-----------------------------+------------------------------------------------------------------------+-------------+--------------------------+
+// If either `backend-json` or `backend` is provided, the `exatn-dm` simulator will simulate the backend noise associated with each quantum gate.
+
+namespace xacc {
+// Forward declaration
+class NoiseModel;
+struct NoiseChannelKraus;
+} // namespace xacc
+namespace tnqvm {
+// Single-site Kraus Op (this Accelerator can only model this)
+struct KrausOp {
+  size_t qubit;
+  // Choi matrix
+  std::vector<std::vector<std::complex<double>>> mats;
+};
+
+class ExaTnDmVisitor : public TNQVMVisitor
+{
+public:
+    // Constructor
+    ExaTnDmVisitor();
+
+    // Virtual function impls:        
+    virtual void initialize(std::shared_ptr<AcceleratorBuffer> buffer, int nbShots) override;
+    virtual void finalize() override;
+
+    // Service name as defined in manifest.json
+    virtual const std::string name() const override { return "exatn-dm"; }
+    virtual const std::string description() const override { return "ExaTN Density Matrix Visitor"; }
+    virtual std::shared_ptr<TNQVMVisitor> clone() override { return std::make_shared<ExaTnDmVisitor>(); }
+
+    // one-qubit gates
+    virtual void visit(Identity& in_IdentityGate) override;
+    virtual void visit(Hadamard& in_HadamardGate) override;
+    virtual void visit(X& in_XGate) override;
+    virtual void visit(Y& in_YGate) override;
+    virtual void visit(Z& in_ZGate) override;
+    virtual void visit(Rx& in_RxGate) override;
+    virtual void visit(Ry& in_RyGate) override;
+    virtual void visit(Rz& in_RzGate) override;
+    virtual void visit(T& in_TGate) override;
+    virtual void visit(Tdg& in_TdgGate) override;
+    virtual void visit(CPhase& in_CPhaseGate) override;
+    virtual void visit(U& in_UGate) override;
+    // two-qubit gates
+    virtual void visit(CNOT& in_CNOTGate) override;
+    virtual void visit(Swap& in_SwapGate) override;
+    virtual void visit(CZ& in_CZGate) override;
+    virtual void visit(iSwap& in_iSwapGate) override;
+    virtual void visit(fSim& in_fsimGate) override;
+    // others
+    virtual void visit(Measure& in_MeasureGate) override;
+
+    virtual const double getExpectationValueZ(std::shared_ptr<CompositeInstruction> in_function) override;
+    class ExaTnTensorFunctor : public talsh::TensorFunctor<exatn::Identifiable>
+    {
+    public:
+        ExaTnTensorFunctor(const std::function<int(talsh::Tensor& in_tensor)>& in_func) :
+            m_func(in_func)
+        {}
+
+        const std::string name() const override { return "default"; }
+        const std::string description() const override { return ""; }
+        virtual void pack(BytePacket& packet) override {};
+        virtual void unpack(BytePacket& packet) override {};
+        virtual int apply(talsh::Tensor& local_tensor) override { return m_func(local_tensor); }
+    private:
+        std::function<int(talsh::Tensor& in_tensor)> m_func;
+    }; 
+
+    [[nodiscard]] exatn::TensorNetwork buildInitialNetwork(size_t in_nbQubits) const;
+    void applySingleQubitGate(xacc::quantum::Gate& in_gateInstruction);
+    void applyTwoQubitGate(xacc::quantum::Gate& in_gateInstruction);
+    void applyNoise(xacc::quantum::Gate &in_gateInstruction);
+    
+  private:
+    exatn::TensorNetwork m_tensorNetwork;
+    std::shared_ptr<AcceleratorBuffer> m_buffer;
+    std::vector<size_t> m_measuredBits;
+    int m_nbShots;
+    int m_tensorIdCounter;
+    std::shared_ptr<xacc::NoiseModel> m_noiseConfig;
+};
+} // namespace tnqvm
\ No newline at end of file
diff --git a/tnqvm/visitors/exatn-dm/manifest.json b/tnqvm/visitors/exatn-dm/manifest.json
new file mode 100644
index 0000000..4379f7b
--- /dev/null
+++ b/tnqvm/visitors/exatn-dm/manifest.json
@@ -0,0 +1,6 @@
+{
+  "bundle.symbolic_name" : "tnqvm_exatn_dm",
+  "bundle.activator" : true,
+  "bundle.name" : "XACC ExaTN Density Matrix Contract backend",
+  "bundle.description" : "This bundle provides a full density matrix tensor network Accelerator for Gate Model QC."
+}
\ No newline at end of file
diff --git a/tnqvm/visitors/exatn-dm/tests/CMakeLists.txt b/tnqvm/visitors/exatn-dm/tests/CMakeLists.txt
new file mode 100644
index 0000000..c2a154e
--- /dev/null
+++ b/tnqvm/visitors/exatn-dm/tests/CMakeLists.txt
@@ -0,0 +1,6 @@
+# Forcing use of XACC gtest install... :/ 
+include_directories(${XACC_ROOT}/include/gtest)
+
+add_executable(DensityMatrixSimTester JsonNoiseModelTester.cpp)
+target_link_libraries(DensityMatrixSimTester PRIVATE ${XACC_ROOT}/lib/libgtest.so ${XACC_ROOT}/lib/libgtest_main.so tnqvm-exatn)
+add_test(NAME DensityMatrixSimTester COMMAND DensityMatrixSimTester)
diff --git a/tnqvm/visitors/exatn-dm/tests/JsonNoiseModelTester.cpp b/tnqvm/visitors/exatn-dm/tests/JsonNoiseModelTester.cpp
new file mode 100644
index 0000000..812e7e5
--- /dev/null
+++ b/tnqvm/visitors/exatn-dm/tests/JsonNoiseModelTester.cpp
@@ -0,0 +1,535 @@
+#include <gtest/gtest.h>
+#include "xacc.hpp"
+#include "xacc_service.hpp"
+#include "NoiseModel.hpp"
+#include "Optimizer.hpp"
+#include "xacc_observable.hpp"
+#include "Algorithm.hpp"
+
+namespace {
+// A sample Json for testing
+// Single-qubit depolarizing:
+const std::string depol_json =
+    R"({"gate_noise": [{"gate_name": "X", "register_location": ["0"], "noise_channels": [{"matrix": [[[[0.99498743710662, 0.0], [0.0, 0.0]], [[0.0, 0.0], [0.99498743710662, 0.0]]], [[[0.0, 0.0], [0.05773502691896258, 0.0]], [[0.05773502691896258, 0.0], [0.0, 0.0]]], [[[0.0, 0.0], [0.0, -0.05773502691896258]], [[0.0, 0.05773502691896258], [0.0, 0.0]]], [[[0.05773502691896258, 0.0], [0.0, 0.0]], [[0.0, 0.0], [-0.05773502691896258, 0.0]]]]}]}], "bit_order": "MSB"})";
+// Single-qubit amplitude damping (25% rate):
+const std::string ad_json =
+    R"({"gate_noise": [{"gate_name": "X", "register_location": ["0"], "noise_channels": [{"matrix": [[[[1.0, 0.0], [0.0, 0.0]], [[0.0, 0.0], [0.8660254037844386, 0.0]]], [[[0.0, 0.0], [0.5, 0.0]], [[0.0, 0.0], [0.0, 0.0]]]]}]}], "bit_order": "MSB"})";
+// Two-qubit noise channel (on a CNOT gate) in MSB and LSB order convention
+// (matrix representation)
+const std::string msb_noise_model =
+    R"({"gate_noise": [{"gate_name": "CNOT", "register_location": ["0", "1"], "noise_channels": [{"matrix": [[[[0.99498743710662, 0.0], [0.0, 0.0], [0.0, 0.0], [0.0, 0.0]], [[0.0, 0.0], [0.99498743710662, 0.0], [0.0, 0.0], [0.0, 0.0]], [[0.0, 0.0], [0.0, 0.0], [0.99498743710662, 0.0], [0.0, 0.0]], [[0.0, 0.0], [0.0, 0.0], [0.0, 0.0], [0.99498743710662, 0.0]]], [[[0.0, 0.0], [0.05773502691896258, 0.0], [0.0, 0.0], [0.0, 0.0]], [[0.05773502691896258, 0.0], [0.0, 0.0], [0.0, 0.0], [0.0, 0.0]], [[0.0, 0.0], [0.0, 0.0], [0.0, 0.0], [0.05773502691896258, 0.0]], [[0.0, 0.0], [0.0, 0.0], [0.05773502691896258, 0.0], [0.0, 0.0]]], [[[0.0, 0.0], [0.0, -0.05773502691896258], [0.0, 0.0], [0.0, 0.0]], [[0.0, 0.05773502691896258], [0.0, 0.0], [0.0, 0.0], [0.0, 0.0]], [[0.0, 0.0], [0.0, 0.0], [0.0, 0.0], [0.0, -0.05773502691896258]], [[0.0, 0.0], [0.0, 0.0], [0.0, 0.05773502691896258], [0.0, 0.0]]], [[[0.05773502691896258, 0.0], [0.0, 0.0], [0.0, 0.0], [0.0, 0.0]], [[0.0, 0.0], [-0.05773502691896258, 0.0], [0.0, 0.0], [-0.0, 0.0]], [[0.0, 0.0], [0.0, 0.0], [0.05773502691896258, 0.0], [0.0, 0.0]], [[0.0, 0.0], [-0.0, 0.0], [0.0, 0.0], [-0.05773502691896258, 0.0]]]]}]}], "bit_order": "MSB"})";
+const std::string lsb_noise_model =
+    R"({"gate_noise": [{"gate_name": "CNOT", "register_location": ["0", "1"], "noise_channels": [{"matrix": [[[[0.99498743710662, 0.0], [0.0, 0.0], [0.0, 0.0], [0.0, 0.0]], [[0.0, 0.0], [0.99498743710662, 0.0], [0.0, 0.0], [0.0, 0.0]], [[0.0, 0.0], [0.0, 0.0], [0.99498743710662, 0.0], [0.0, 0.0]], [[0.0, 0.0], [0.0, 0.0], [0.0, 0.0], [0.99498743710662, 0.0]]], [[[0.0, 0.0], [0.0, 0.0], [0.05773502691896258, 0.0], [0.0, 0.0]], [[0.0, 0.0], [0.0, 0.0], [0.0, 0.0], [0.05773502691896258, 0.0]], [[0.05773502691896258, 0.0], [0.0, 0.0], [0.0, 0.0], [0.0, 0.0]], [[0.0, 0.0], [0.05773502691896258, 0.0], [0.0, 0.0], [0.0, 0.0]]], [[[0.0, 0.0], [0.0, 0.0], [0.0, -0.05773502691896258], [0.0, 0.0]], [[0.0, 0.0], [0.0, 0.0], [0.0, 0.0], [0.0, -0.05773502691896258]], [[0.0, 0.05773502691896258], [0.0, 0.0], [0.0, 0.0], [0.0, 0.0]], [[0.0, 0.0], [0.0, 0.05773502691896258], [0.0, 0.0], [0.0, 0.0]]], [[[0.05773502691896258, 0.0], [0.0, 0.0], [0.0, 0.0], [0.0, 0.0]], [[0.0, 0.0], [0.05773502691896258, 0.0], [0.0, 0.0], [0.0, 0.0]], [[0.0, 0.0], [0.0, 0.0], [-0.05773502691896258, 0.0], [-0.0, 0.0]], [[0.0, 0.0], [0.0, 0.0], [-0.0, 0.0], [-0.05773502691896258, 0.0]]]]}]}], "bit_order": "LSB"})";
+// Noise model that only has readout errors for validation:
+// P(1|0) = 0.1; P(0|1) = 0.2
+const std::string ro_error_noise_model =
+    R"({"gate_noise": [], "bit_order": "MSB", "readout_errors": [{"register_location": "0", "prob_meas0_prep1": 0.2, "prob_meas1_prep0": 0.1}]})";
+
+// Multiple qubit depolarizing:
+const std::string noise_model_multi_depol =
+    R"({"gate_noise": [{"gate_name": "CNOT", "register_location": ["0", "1"], "noise_channels": [{"matrix": [[[[0.99498743710662, 0.0], [0.0, 0.0], [0.0, 0.0], [0.0, 0.0]], [[0.0, 0.0], [0.99498743710662, 0.0], [0.0, 0.0], [0.0, 0.0]], [[0.0, 0.0], [0.0, 0.0], [0.99498743710662, 0.0], [0.0, 0.0]], [[0.0, 0.0], [0.0, 0.0], [0.0, 0.0], [0.99498743710662, 0.0]]], [[[0.0, 0.0], [0.0, 0.0], [0.0, 0.0], [0.05773502691896258, 0.0]], [[0.0, 0.0], [0.0, 0.0], [0.05773502691896258, 0.0], [0.0, 0.0]], [[0.0, 0.0], [0.05773502691896258, 0.0], [0.0, 0.0], [0.0, 0.0]], [[0.05773502691896258, 0.0], [0.0, 0.0], [0.0, 0.0], [0.0, 0.0]]], [[[0.0, 0.0], [0.0, 0.0], [0.0, 0.0], [-0.05773502691896258, 0.0]], [[0.0, 0.0], [0.0, 0.0], [0.05773502691896258, 0.0], [0.0, 0.0]], [[0.0, 0.0], [0.05773502691896258, 0.0], [0.0, 0.0], [0.0, 0.0]], [[-0.05773502691896258, 0.0], [0.0, 0.0], [0.0, 0.0], [0.0, 0.0]]], [[[0.05773502691896258, 0.0], [0.0, 0.0], [0.0, 0.0], [0.0, 0.0]], [[0.0, 0.0], [-0.05773502691896258, 0.0], [0.0, 0.0], [0.0, 0.0]], [[0.0, 0.0], [0.0, 0.0], [-0.05773502691896258, 0.0], [0.0, 0.0]], [[0.0, 0.0], [0.0, 0.0], [0.0, 0.0], [0.05773502691896258, 0.0]]]]}]}], "bit_order": "MSB"})";
+
+// Extremely-weak noise, just to check DM tensor network contraction.
+const std::string noise_model_weak_noise =
+    R"({"gate_noise": [{"gate_name": "H", "register_location": ["0"], "noise_channels": [{"matrix": [[[[0.9999999999995, 0.0], [0.0, 0.0]], [[0.0, 0.0], [0.9999999999995, 0.0]]], [[[0.0, 0.0], [5.773502691896258e-07, 0.0]], [[5.773502691896258e-07, 0.0], [0.0, 0.0]]], [[[0.0, 0.0], [0.0, -5.773502691896258e-07]], [[0.0, 5.773502691896258e-07], [0.0, 0.0]]], [[[5.773502691896258e-07, 0.0], [0.0, 0.0]], [[0.0, 0.0], [-5.773502691896258e-07, 0.0]]]]}]}, {"gate_name": "X", "register_location": ["0"], "noise_channels": [{"matrix": [[[[0.9999999999995, 0.0], [0.0, 0.0]], [[0.0, 0.0], [0.9999999999995, 0.0]]], [[[0.0, 0.0], [5.773502691896258e-07, 0.0]], [[5.773502691896258e-07, 0.0], [0.0, 0.0]]], [[[0.0, 0.0], [0.0, -5.773502691896258e-07]], [[0.0, 5.773502691896258e-07], [0.0, 0.0]]], [[[5.773502691896258e-07, 0.0], [0.0, 0.0]], [[0.0, 0.0], [-5.773502691896258e-07, 0.0]]]]}]}, {"gate_name": "H", "register_location": ["1"], "noise_channels": [{"matrix": [[[[0.9999999999995, 0.0], [0.0, 0.0]], [[0.0, 0.0], [0.9999999999995, 0.0]]], [[[0.0, 0.0], [5.773502691896258e-07, 0.0]], [[5.773502691896258e-07, 0.0], [0.0, 0.0]]], [[[0.0, 0.0], [0.0, -5.773502691896258e-07]], [[0.0, 5.773502691896258e-07], [0.0, 0.0]]], [[[5.773502691896258e-07, 0.0], [0.0, 0.0]], [[0.0, 0.0], [-5.773502691896258e-07, 0.0]]]]}]}, {"gate_name": "X", "register_location": ["1"], "noise_channels": [{"matrix": [[[[0.9999999999995, 0.0], [0.0, 0.0]], [[0.0, 0.0], [0.9999999999995, 0.0]]], [[[0.0, 0.0], [5.773502691896258e-07, 0.0]], [[5.773502691896258e-07, 0.0], [0.0, 0.0]]], [[[0.0, 0.0], [0.0, -5.773502691896258e-07]], [[0.0, 5.773502691896258e-07], [0.0, 0.0]]], [[[5.773502691896258e-07, 0.0], [0.0, 0.0]], [[0.0, 0.0], [-5.773502691896258e-07, 0.0]]]]}]}, {"gate_name": "H", "register_location": ["2"], "noise_channels": [{"matrix": [[[[0.9999999999995, 0.0], [0.0, 0.0]], [[0.0, 0.0], [0.9999999999995, 0.0]]], [[[0.0, 0.0], [5.773502691896258e-07, 0.0]], [[5.773502691896258e-07, 0.0], [0.0, 0.0]]], [[[0.0, 0.0], [0.0, -5.773502691896258e-07]], [[0.0, 5.773502691896258e-07], [0.0, 0.0]]], [[[5.773502691896258e-07, 0.0], [0.0, 0.0]], [[0.0, 0.0], [-5.773502691896258e-07, 0.0]]]]}]}, {"gate_name": "X", "register_location": ["2"], "noise_channels": [{"matrix": [[[[0.9999999999995, 0.0], [0.0, 0.0]], [[0.0, 0.0], [0.9999999999995, 0.0]]], [[[0.0, 0.0], [5.773502691896258e-07, 0.0]], [[5.773502691896258e-07, 0.0], [0.0, 0.0]]], [[[0.0, 0.0], [0.0, -5.773502691896258e-07]], [[0.0, 5.773502691896258e-07], [0.0, 0.0]]], [[[5.773502691896258e-07, 0.0], [0.0, 0.0]], [[0.0, 0.0], [-5.773502691896258e-07, 0.0]]]]}]}, {"gate_name": "H", "register_location": ["3"], "noise_channels": [{"matrix": [[[[0.9999999999995, 0.0], [0.0, 0.0]], [[0.0, 0.0], [0.9999999999995, 0.0]]], [[[0.0, 0.0], [5.773502691896258e-07, 0.0]], [[5.773502691896258e-07, 0.0], [0.0, 0.0]]], [[[0.0, 0.0], [0.0, -5.773502691896258e-07]], [[0.0, 5.773502691896258e-07], [0.0, 0.0]]], [[[5.773502691896258e-07, 0.0], [0.0, 0.0]], [[0.0, 0.0], [-5.773502691896258e-07, 0.0]]]]}]}, {"gate_name": "X", "register_location": ["3"], "noise_channels": [{"matrix": [[[[0.9999999999995, 0.0], [0.0, 0.0]], [[0.0, 0.0], [0.9999999999995, 0.0]]], [[[0.0, 0.0], [5.773502691896258e-07, 0.0]], [[5.773502691896258e-07, 0.0], [0.0, 0.0]]], [[[0.0, 0.0], [0.0, -5.773502691896258e-07]], [[0.0, 5.773502691896258e-07], [0.0, 0.0]]], [[[5.773502691896258e-07, 0.0], [0.0, 0.0]], [[0.0, 0.0], [-5.773502691896258e-07, 0.0]]]]}]}], "bit_order": "MSB"})";
+
+// Complex noise model: 4 qubits, X and H and CNOT
+const std::string noise_model_4q =
+    R"({"gate_noise": [{"gate_name": "H", "register_location": ["0"], "noise_channels": [{"matrix": [[[[0.99498743710662, 0.0], [0.0, 0.0]], [[0.0, 0.0], [0.99498743710662, 0.0]]], [[[0.0, 0.0], [0.05773502691896258, 0.0]], [[0.05773502691896258, 0.0], [0.0, 0.0]]], [[[0.0, 0.0], [0.0, -0.05773502691896258]], [[0.0, 0.05773502691896258], [0.0, 0.0]]], [[[0.05773502691896258, 0.0], [0.0, 0.0]], [[0.0, 0.0], [-0.05773502691896258, 0.0]]]]}]}, {"gate_name": "X", "register_location": ["0"], "noise_channels": [{"matrix": [[[[0.99498743710662, 0.0], [0.0, 0.0]], [[0.0, 0.0], [0.99498743710662, 0.0]]], [[[0.0, 0.0], [0.05773502691896258, 0.0]], [[0.05773502691896258, 0.0], [0.0, 0.0]]], [[[0.0, 0.0], [0.0, -0.05773502691896258]], [[0.0, 0.05773502691896258], [0.0, 0.0]]], [[[0.05773502691896258, 0.0], [0.0, 0.0]], [[0.0, 0.0], [-0.05773502691896258, 0.0]]]]}]}, {"gate_name": "H", "register_location": ["1"], "noise_channels": [{"matrix": [[[[0.99498743710662, 0.0], [0.0, 0.0]], [[0.0, 0.0], [0.99498743710662, 0.0]]], [[[0.0, 0.0], [0.05773502691896258, 0.0]], [[0.05773502691896258, 0.0], [0.0, 0.0]]], [[[0.0, 0.0], [0.0, -0.05773502691896258]], [[0.0, 0.05773502691896258], [0.0, 0.0]]], [[[0.05773502691896258, 0.0], [0.0, 0.0]], [[0.0, 0.0], [-0.05773502691896258, 0.0]]]]}]}, {"gate_name": "X", "register_location": ["1"], "noise_channels": [{"matrix": [[[[0.99498743710662, 0.0], [0.0, 0.0]], [[0.0, 0.0], [0.99498743710662, 0.0]]], [[[0.0, 0.0], [0.05773502691896258, 0.0]], [[0.05773502691896258, 0.0], [0.0, 0.0]]], [[[0.0, 0.0], [0.0, -0.05773502691896258]], [[0.0, 0.05773502691896258], [0.0, 0.0]]], [[[0.05773502691896258, 0.0], [0.0, 0.0]], [[0.0, 0.0], [-0.05773502691896258, 0.0]]]]}]}, {"gate_name": "H", "register_location": ["2"], "noise_channels": [{"matrix": [[[[0.99498743710662, 0.0], [0.0, 0.0]], [[0.0, 0.0], [0.99498743710662, 0.0]]], [[[0.0, 0.0], [0.05773502691896258, 0.0]], [[0.05773502691896258, 0.0], [0.0, 0.0]]], [[[0.0, 0.0], [0.0, -0.05773502691896258]], [[0.0, 0.05773502691896258], [0.0, 0.0]]], [[[0.05773502691896258, 0.0], [0.0, 0.0]], [[0.0, 0.0], [-0.05773502691896258, 0.0]]]]}]}, {"gate_name": "X", "register_location": ["2"], "noise_channels": [{"matrix": [[[[0.99498743710662, 0.0], [0.0, 0.0]], [[0.0, 0.0], [0.99498743710662, 0.0]]], [[[0.0, 0.0], [0.05773502691896258, 0.0]], [[0.05773502691896258, 0.0], [0.0, 0.0]]], [[[0.0, 0.0], [0.0, -0.05773502691896258]], [[0.0, 0.05773502691896258], [0.0, 0.0]]], [[[0.05773502691896258, 0.0], [0.0, 0.0]], [[0.0, 0.0], [-0.05773502691896258, 0.0]]]]}]}, {"gate_name": "H", "register_location": ["3"], "noise_channels": [{"matrix": [[[[0.99498743710662, 0.0], [0.0, 0.0]], [[0.0, 0.0], [0.99498743710662, 0.0]]], [[[0.0, 0.0], [0.05773502691896258, 0.0]], [[0.05773502691896258, 0.0], [0.0, 0.0]]], [[[0.0, 0.0], [0.0, -0.05773502691896258]], [[0.0, 0.05773502691896258], [0.0, 0.0]]], [[[0.05773502691896258, 0.0], [0.0, 0.0]], [[0.0, 0.0], [-0.05773502691896258, 0.0]]]]}]}, {"gate_name": "X", "register_location": ["3"], "noise_channels": [{"matrix": [[[[0.99498743710662, 0.0], [0.0, 0.0]], [[0.0, 0.0], [0.99498743710662, 0.0]]], [[[0.0, 0.0], [0.05773502691896258, 0.0]], [[0.05773502691896258, 0.0], [0.0, 0.0]]], [[[0.0, 0.0], [0.0, -0.05773502691896258]], [[0.0, 0.05773502691896258], [0.0, 0.0]]], [[[0.05773502691896258, 0.0], [0.0, 0.0]], [[0.0, 0.0], [-0.05773502691896258, 0.0]]]]}]}, {"gate_name": "CNOT", "register_location": ["0", "1"], "noise_channels": [{"matrix": [[[[0.99498743710662, 0.0], [0.0, 0.0], [0.0, 0.0], [0.0, 0.0]], [[0.0, 0.0], [0.99498743710662, 0.0], [0.0, 0.0], [0.0, 0.0]], [[0.0, 0.0], [0.0, 0.0], [0.99498743710662, 0.0], [0.0, 0.0]], [[0.0, 0.0], [0.0, 0.0], [0.0, 0.0], [0.99498743710662, 0.0]]], [[[0.0, 0.0], [0.0, 0.0], [0.0, 0.0], [0.05773502691896258, 0.0]], [[0.0, 0.0], [0.0, 0.0], [0.05773502691896258, 0.0], [0.0, 0.0]], [[0.0, 0.0], [0.05773502691896258, 0.0], [0.0, 0.0], [0.0, 0.0]], [[0.05773502691896258, 0.0], [0.0, 0.0], [0.0, 0.0], [0.0, 0.0]]], [[[0.0, 0.0], [0.0, 0.0], [0.0, 0.0], [-0.05773502691896258, 0.0]], [[0.0, 0.0], [0.0, 0.0], [0.05773502691896258, 0.0], [0.0, 0.0]], [[0.0, 0.0], [0.05773502691896258, 0.0], [0.0, 0.0], [0.0, 0.0]], [[-0.05773502691896258, 0.0], [0.0, 0.0], [0.0, 0.0], [0.0, 0.0]]], [[[0.05773502691896258, 0.0], [0.0, 0.0], [0.0, 0.0], [0.0, 0.0]], [[0.0, 0.0], [-0.05773502691896258, 0.0], [0.0, 0.0], [0.0, 0.0]], [[0.0, 0.0], [0.0, 0.0], [-0.05773502691896258, 0.0], [0.0, 0.0]], [[0.0, 0.0], [0.0, 0.0], [0.0, 0.0], [0.05773502691896258, 0.0]]]]}]}, {"gate_name": "CNOT", "register_location": ["0", "2"], "noise_channels": [{"matrix": [[[[0.99498743710662, 0.0], [0.0, 0.0], [0.0, 0.0], [0.0, 0.0]], [[0.0, 0.0], [0.99498743710662, 0.0], [0.0, 0.0], [0.0, 0.0]], [[0.0, 0.0], [0.0, 0.0], [0.99498743710662, 0.0], [0.0, 0.0]], [[0.0, 0.0], [0.0, 0.0], [0.0, 0.0], [0.99498743710662, 0.0]]], [[[0.0, 0.0], [0.0, 0.0], [0.0, 0.0], [0.05773502691896258, 0.0]], [[0.0, 0.0], [0.0, 0.0], [0.05773502691896258, 0.0], [0.0, 0.0]], [[0.0, 0.0], [0.05773502691896258, 0.0], [0.0, 0.0], [0.0, 0.0]], [[0.05773502691896258, 0.0], [0.0, 0.0], [0.0, 0.0], [0.0, 0.0]]], [[[0.0, 0.0], [0.0, 0.0], [0.0, 0.0], [-0.05773502691896258, 0.0]], [[0.0, 0.0], [0.0, 0.0], [0.05773502691896258, 0.0], [0.0, 0.0]], [[0.0, 0.0], [0.05773502691896258, 0.0], [0.0, 0.0], [0.0, 0.0]], [[-0.05773502691896258, 0.0], [0.0, 0.0], [0.0, 0.0], [0.0, 0.0]]], [[[0.05773502691896258, 0.0], [0.0, 0.0], [0.0, 0.0], [0.0, 0.0]], [[0.0, 0.0], [-0.05773502691896258, 0.0], [0.0, 0.0], [0.0, 0.0]], [[0.0, 0.0], [0.0, 0.0], [-0.05773502691896258, 0.0], [0.0, 0.0]], [[0.0, 0.0], [0.0, 0.0], [0.0, 0.0], [0.05773502691896258, 0.0]]]]}]}, {"gate_name": "CNOT", "register_location": ["0", "3"], "noise_channels": [{"matrix": [[[[0.99498743710662, 0.0], [0.0, 0.0], [0.0, 0.0], [0.0, 0.0]], [[0.0, 0.0], [0.99498743710662, 0.0], [0.0, 0.0], [0.0, 0.0]], [[0.0, 0.0], [0.0, 0.0], [0.99498743710662, 0.0], [0.0, 0.0]], [[0.0, 0.0], [0.0, 0.0], [0.0, 0.0], [0.99498743710662, 0.0]]], [[[0.0, 0.0], [0.0, 0.0], [0.0, 0.0], [0.05773502691896258, 0.0]], [[0.0, 0.0], [0.0, 0.0], [0.05773502691896258, 0.0], [0.0, 0.0]], [[0.0, 0.0], [0.05773502691896258, 0.0], [0.0, 0.0], [0.0, 0.0]], [[0.05773502691896258, 0.0], [0.0, 0.0], [0.0, 0.0], [0.0, 0.0]]], [[[0.0, 0.0], [0.0, 0.0], [0.0, 0.0], [-0.05773502691896258, 0.0]], [[0.0, 0.0], [0.0, 0.0], [0.05773502691896258, 0.0], [0.0, 0.0]], [[0.0, 0.0], [0.05773502691896258, 0.0], [0.0, 0.0], [0.0, 0.0]], [[-0.05773502691896258, 0.0], [0.0, 0.0], [0.0, 0.0], [0.0, 0.0]]], [[[0.05773502691896258, 0.0], [0.0, 0.0], [0.0, 0.0], [0.0, 0.0]], [[0.0, 0.0], [-0.05773502691896258, 0.0], [0.0, 0.0], [0.0, 0.0]], [[0.0, 0.0], [0.0, 0.0], [-0.05773502691896258, 0.0], [0.0, 0.0]], [[0.0, 0.0], [0.0, 0.0], [0.0, 0.0], [0.05773502691896258, 0.0]]]]}]}, {"gate_name": "CNOT", "register_location": ["1", "0"], "noise_channels": [{"matrix": [[[[0.99498743710662, 0.0], [0.0, 0.0], [0.0, 0.0], [0.0, 0.0]], [[0.0, 0.0], [0.99498743710662, 0.0], [0.0, 0.0], [0.0, 0.0]], [[0.0, 0.0], [0.0, 0.0], [0.99498743710662, 0.0], [0.0, 0.0]], [[0.0, 0.0], [0.0, 0.0], [0.0, 0.0], [0.99498743710662, 0.0]]], [[[0.0, 0.0], [0.0, 0.0], [0.0, 0.0], [0.05773502691896258, 0.0]], [[0.0, 0.0], [0.0, 0.0], [0.05773502691896258, 0.0], [0.0, 0.0]], [[0.0, 0.0], [0.05773502691896258, 0.0], [0.0, 0.0], [0.0, 0.0]], [[0.05773502691896258, 0.0], [0.0, 0.0], [0.0, 0.0], [0.0, 0.0]]], [[[0.0, 0.0], [0.0, 0.0], [0.0, 0.0], [-0.05773502691896258, 0.0]], [[0.0, 0.0], [0.0, 0.0], [0.05773502691896258, 0.0], [0.0, 0.0]], [[0.0, 0.0], [0.05773502691896258, 0.0], [0.0, 0.0], [0.0, 0.0]], [[-0.05773502691896258, 0.0], [0.0, 0.0], [0.0, 0.0], [0.0, 0.0]]], [[[0.05773502691896258, 0.0], [0.0, 0.0], [0.0, 0.0], [0.0, 0.0]], [[0.0, 0.0], [-0.05773502691896258, 0.0], [0.0, 0.0], [0.0, 0.0]], [[0.0, 0.0], [0.0, 0.0], [-0.05773502691896258, 0.0], [0.0, 0.0]], [[0.0, 0.0], [0.0, 0.0], [0.0, 0.0], [0.05773502691896258, 0.0]]]]}]}, {"gate_name": "CNOT", "register_location": ["1", "2"], "noise_channels": [{"matrix": [[[[0.99498743710662, 0.0], [0.0, 0.0], [0.0, 0.0], [0.0, 0.0]], [[0.0, 0.0], [0.99498743710662, 0.0], [0.0, 0.0], [0.0, 0.0]], [[0.0, 0.0], [0.0, 0.0], [0.99498743710662, 0.0], [0.0, 0.0]], [[0.0, 0.0], [0.0, 0.0], [0.0, 0.0], [0.99498743710662, 0.0]]], [[[0.0, 0.0], [0.0, 0.0], [0.0, 0.0], [0.05773502691896258, 0.0]], [[0.0, 0.0], [0.0, 0.0], [0.05773502691896258, 0.0], [0.0, 0.0]], [[0.0, 0.0], [0.05773502691896258, 0.0], [0.0, 0.0], [0.0, 0.0]], [[0.05773502691896258, 0.0], [0.0, 0.0], [0.0, 0.0], [0.0, 0.0]]], [[[0.0, 0.0], [0.0, 0.0], [0.0, 0.0], [-0.05773502691896258, 0.0]], [[0.0, 0.0], [0.0, 0.0], [0.05773502691896258, 0.0], [0.0, 0.0]], [[0.0, 0.0], [0.05773502691896258, 0.0], [0.0, 0.0], [0.0, 0.0]], [[-0.05773502691896258, 0.0], [0.0, 0.0], [0.0, 0.0], [0.0, 0.0]]], [[[0.05773502691896258, 0.0], [0.0, 0.0], [0.0, 0.0], [0.0, 0.0]], [[0.0, 0.0], [-0.05773502691896258, 0.0], [0.0, 0.0], [0.0, 0.0]], [[0.0, 0.0], [0.0, 0.0], [-0.05773502691896258, 0.0], [0.0, 0.0]], [[0.0, 0.0], [0.0, 0.0], [0.0, 0.0], [0.05773502691896258, 0.0]]]]}]}, {"gate_name": "CNOT", "register_location": ["1", "3"], "noise_channels": [{"matrix": [[[[0.99498743710662, 0.0], [0.0, 0.0], [0.0, 0.0], [0.0, 0.0]], [[0.0, 0.0], [0.99498743710662, 0.0], [0.0, 0.0], [0.0, 0.0]], [[0.0, 0.0], [0.0, 0.0], [0.99498743710662, 0.0], [0.0, 0.0]], [[0.0, 0.0], [0.0, 0.0], [0.0, 0.0], [0.99498743710662, 0.0]]], [[[0.0, 0.0], [0.0, 0.0], [0.0, 0.0], [0.05773502691896258, 0.0]], [[0.0, 0.0], [0.0, 0.0], [0.05773502691896258, 0.0], [0.0, 0.0]], [[0.0, 0.0], [0.05773502691896258, 0.0], [0.0, 0.0], [0.0, 0.0]], [[0.05773502691896258, 0.0], [0.0, 0.0], [0.0, 0.0], [0.0, 0.0]]], [[[0.0, 0.0], [0.0, 0.0], [0.0, 0.0], [-0.05773502691896258, 0.0]], [[0.0, 0.0], [0.0, 0.0], [0.05773502691896258, 0.0], [0.0, 0.0]], [[0.0, 0.0], [0.05773502691896258, 0.0], [0.0, 0.0], [0.0, 0.0]], [[-0.05773502691896258, 0.0], [0.0, 0.0], [0.0, 0.0], [0.0, 0.0]]], [[[0.05773502691896258, 0.0], [0.0, 0.0], [0.0, 0.0], [0.0, 0.0]], [[0.0, 0.0], [-0.05773502691896258, 0.0], [0.0, 0.0], [0.0, 0.0]], [[0.0, 0.0], [0.0, 0.0], [-0.05773502691896258, 0.0], [0.0, 0.0]], [[0.0, 0.0], [0.0, 0.0], [0.0, 0.0], [0.05773502691896258, 0.0]]]]}]}, {"gate_name": "CNOT", "register_location": ["2", "0"], "noise_channels": [{"matrix": [[[[0.99498743710662, 0.0], [0.0, 0.0], [0.0, 0.0], [0.0, 0.0]], [[0.0, 0.0], [0.99498743710662, 0.0], [0.0, 0.0], [0.0, 0.0]], [[0.0, 0.0], [0.0, 0.0], [0.99498743710662, 0.0], [0.0, 0.0]], [[0.0, 0.0], [0.0, 0.0], [0.0, 0.0], [0.99498743710662, 0.0]]], [[[0.0, 0.0], [0.0, 0.0], [0.0, 0.0], [0.05773502691896258, 0.0]], [[0.0, 0.0], [0.0, 0.0], [0.05773502691896258, 0.0], [0.0, 0.0]], [[0.0, 0.0], [0.05773502691896258, 0.0], [0.0, 0.0], [0.0, 0.0]], [[0.05773502691896258, 0.0], [0.0, 0.0], [0.0, 0.0], [0.0, 0.0]]], [[[0.0, 0.0], [0.0, 0.0], [0.0, 0.0], [-0.05773502691896258, 0.0]], [[0.0, 0.0], [0.0, 0.0], [0.05773502691896258, 0.0], [0.0, 0.0]], [[0.0, 0.0], [0.05773502691896258, 0.0], [0.0, 0.0], [0.0, 0.0]], [[-0.05773502691896258, 0.0], [0.0, 0.0], [0.0, 0.0], [0.0, 0.0]]], [[[0.05773502691896258, 0.0], [0.0, 0.0], [0.0, 0.0], [0.0, 0.0]], [[0.0, 0.0], [-0.05773502691896258, 0.0], [0.0, 0.0], [0.0, 0.0]], [[0.0, 0.0], [0.0, 0.0], [-0.05773502691896258, 0.0], [0.0, 0.0]], [[0.0, 0.0], [0.0, 0.0], [0.0, 0.0], [0.05773502691896258, 0.0]]]]}]}, {"gate_name": "CNOT", "register_location": ["2", "1"], "noise_channels": [{"matrix": [[[[0.99498743710662, 0.0], [0.0, 0.0], [0.0, 0.0], [0.0, 0.0]], [[0.0, 0.0], [0.99498743710662, 0.0], [0.0, 0.0], [0.0, 0.0]], [[0.0, 0.0], [0.0, 0.0], [0.99498743710662, 0.0], [0.0, 0.0]], [[0.0, 0.0], [0.0, 0.0], [0.0, 0.0], [0.99498743710662, 0.0]]], [[[0.0, 0.0], [0.0, 0.0], [0.0, 0.0], [0.05773502691896258, 0.0]], [[0.0, 0.0], [0.0, 0.0], [0.05773502691896258, 0.0], [0.0, 0.0]], [[0.0, 0.0], [0.05773502691896258, 0.0], [0.0, 0.0], [0.0, 0.0]], [[0.05773502691896258, 0.0], [0.0, 0.0], [0.0, 0.0], [0.0, 0.0]]], [[[0.0, 0.0], [0.0, 0.0], [0.0, 0.0], [-0.05773502691896258, 0.0]], [[0.0, 0.0], [0.0, 0.0], [0.05773502691896258, 0.0], [0.0, 0.0]], [[0.0, 0.0], [0.05773502691896258, 0.0], [0.0, 0.0], [0.0, 0.0]], [[-0.05773502691896258, 0.0], [0.0, 0.0], [0.0, 0.0], [0.0, 0.0]]], [[[0.05773502691896258, 0.0], [0.0, 0.0], [0.0, 0.0], [0.0, 0.0]], [[0.0, 0.0], [-0.05773502691896258, 0.0], [0.0, 0.0], [0.0, 0.0]], [[0.0, 0.0], [0.0, 0.0], [-0.05773502691896258, 0.0], [0.0, 0.0]], [[0.0, 0.0], [0.0, 0.0], [0.0, 0.0], [0.05773502691896258, 0.0]]]]}]}, {"gate_name": "CNOT", "register_location": ["2", "3"], "noise_channels": [{"matrix": [[[[0.99498743710662, 0.0], [0.0, 0.0], [0.0, 0.0], [0.0, 0.0]], [[0.0, 0.0], [0.99498743710662, 0.0], [0.0, 0.0], [0.0, 0.0]], [[0.0, 0.0], [0.0, 0.0], [0.99498743710662, 0.0], [0.0, 0.0]], [[0.0, 0.0], [0.0, 0.0], [0.0, 0.0], [0.99498743710662, 0.0]]], [[[0.0, 0.0], [0.0, 0.0], [0.0, 0.0], [0.05773502691896258, 0.0]], [[0.0, 0.0], [0.0, 0.0], [0.05773502691896258, 0.0], [0.0, 0.0]], [[0.0, 0.0], [0.05773502691896258, 0.0], [0.0, 0.0], [0.0, 0.0]], [[0.05773502691896258, 0.0], [0.0, 0.0], [0.0, 0.0], [0.0, 0.0]]], [[[0.0, 0.0], [0.0, 0.0], [0.0, 0.0], [-0.05773502691896258, 0.0]], [[0.0, 0.0], [0.0, 0.0], [0.05773502691896258, 0.0], [0.0, 0.0]], [[0.0, 0.0], [0.05773502691896258, 0.0], [0.0, 0.0], [0.0, 0.0]], [[-0.05773502691896258, 0.0], [0.0, 0.0], [0.0, 0.0], [0.0, 0.0]]], [[[0.05773502691896258, 0.0], [0.0, 0.0], [0.0, 0.0], [0.0, 0.0]], [[0.0, 0.0], [-0.05773502691896258, 0.0], [0.0, 0.0], [0.0, 0.0]], [[0.0, 0.0], [0.0, 0.0], [-0.05773502691896258, 0.0], [0.0, 0.0]], [[0.0, 0.0], [0.0, 0.0], [0.0, 0.0], [0.05773502691896258, 0.0]]]]}]}, {"gate_name": "CNOT", "register_location": ["3", "0"], "noise_channels": [{"matrix": [[[[0.99498743710662, 0.0], [0.0, 0.0], [0.0, 0.0], [0.0, 0.0]], [[0.0, 0.0], [0.99498743710662, 0.0], [0.0, 0.0], [0.0, 0.0]], [[0.0, 0.0], [0.0, 0.0], [0.99498743710662, 0.0], [0.0, 0.0]], [[0.0, 0.0], [0.0, 0.0], [0.0, 0.0], [0.99498743710662, 0.0]]], [[[0.0, 0.0], [0.0, 0.0], [0.0, 0.0], [0.05773502691896258, 0.0]], [[0.0, 0.0], [0.0, 0.0], [0.05773502691896258, 0.0], [0.0, 0.0]], [[0.0, 0.0], [0.05773502691896258, 0.0], [0.0, 0.0], [0.0, 0.0]], [[0.05773502691896258, 0.0], [0.0, 0.0], [0.0, 0.0], [0.0, 0.0]]], [[[0.0, 0.0], [0.0, 0.0], [0.0, 0.0], [-0.05773502691896258, 0.0]], [[0.0, 0.0], [0.0, 0.0], [0.05773502691896258, 0.0], [0.0, 0.0]], [[0.0, 0.0], [0.05773502691896258, 0.0], [0.0, 0.0], [0.0, 0.0]], [[-0.05773502691896258, 0.0], [0.0, 0.0], [0.0, 0.0], [0.0, 0.0]]], [[[0.05773502691896258, 0.0], [0.0, 0.0], [0.0, 0.0], [0.0, 0.0]], [[0.0, 0.0], [-0.05773502691896258, 0.0], [0.0, 0.0], [0.0, 0.0]], [[0.0, 0.0], [0.0, 0.0], [-0.05773502691896258, 0.0], [0.0, 0.0]], [[0.0, 0.0], [0.0, 0.0], [0.0, 0.0], [0.05773502691896258, 0.0]]]]}]}, {"gate_name": "CNOT", "register_location": ["3", "1"], "noise_channels": [{"matrix": [[[[0.99498743710662, 0.0], [0.0, 0.0], [0.0, 0.0], [0.0, 0.0]], [[0.0, 0.0], [0.99498743710662, 0.0], [0.0, 0.0], [0.0, 0.0]], [[0.0, 0.0], [0.0, 0.0], [0.99498743710662, 0.0], [0.0, 0.0]], [[0.0, 0.0], [0.0, 0.0], [0.0, 0.0], [0.99498743710662, 0.0]]], [[[0.0, 0.0], [0.0, 0.0], [0.0, 0.0], [0.05773502691896258, 0.0]], [[0.0, 0.0], [0.0, 0.0], [0.05773502691896258, 0.0], [0.0, 0.0]], [[0.0, 0.0], [0.05773502691896258, 0.0], [0.0, 0.0], [0.0, 0.0]], [[0.05773502691896258, 0.0], [0.0, 0.0], [0.0, 0.0], [0.0, 0.0]]], [[[0.0, 0.0], [0.0, 0.0], [0.0, 0.0], [-0.05773502691896258, 0.0]], [[0.0, 0.0], [0.0, 0.0], [0.05773502691896258, 0.0], [0.0, 0.0]], [[0.0, 0.0], [0.05773502691896258, 0.0], [0.0, 0.0], [0.0, 0.0]], [[-0.05773502691896258, 0.0], [0.0, 0.0], [0.0, 0.0], [0.0, 0.0]]], [[[0.05773502691896258, 0.0], [0.0, 0.0], [0.0, 0.0], [0.0, 0.0]], [[0.0, 0.0], [-0.05773502691896258, 0.0], [0.0, 0.0], [0.0, 0.0]], [[0.0, 0.0], [0.0, 0.0], [-0.05773502691896258, 0.0], [0.0, 0.0]], [[0.0, 0.0], [0.0, 0.0], [0.0, 0.0], [0.05773502691896258, 0.0]]]]}]}, {"gate_name": "CNOT", "register_location": ["3", "2"], "noise_channels": [{"matrix": [[[[0.99498743710662, 0.0], [0.0, 0.0], [0.0, 0.0], [0.0, 0.0]], [[0.0, 0.0], [0.99498743710662, 0.0], [0.0, 0.0], [0.0, 0.0]], [[0.0, 0.0], [0.0, 0.0], [0.99498743710662, 0.0], [0.0, 0.0]], [[0.0, 0.0], [0.0, 0.0], [0.0, 0.0], [0.99498743710662, 0.0]]], [[[0.0, 0.0], [0.0, 0.0], [0.0, 0.0], [0.05773502691896258, 0.0]], [[0.0, 0.0], [0.0, 0.0], [0.05773502691896258, 0.0], [0.0, 0.0]], [[0.0, 0.0], [0.05773502691896258, 0.0], [0.0, 0.0], [0.0, 0.0]], [[0.05773502691896258, 0.0], [0.0, 0.0], [0.0, 0.0], [0.0, 0.0]]], [[[0.0, 0.0], [0.0, 0.0], [0.0, 0.0], [-0.05773502691896258, 0.0]], [[0.0, 0.0], [0.0, 0.0], [0.05773502691896258, 0.0], [0.0, 0.0]], [[0.0, 0.0], [0.05773502691896258, 0.0], [0.0, 0.0], [0.0, 0.0]], [[-0.05773502691896258, 0.0], [0.0, 0.0], [0.0, 0.0], [0.0, 0.0]]], [[[0.05773502691896258, 0.0], [0.0, 0.0], [0.0, 0.0], [0.0, 0.0]], [[0.0, 0.0], [-0.05773502691896258, 0.0], [0.0, 0.0], [0.0, 0.0]], [[0.0, 0.0], [0.0, 0.0], [-0.05773502691896258, 0.0], [0.0, 0.0]], [[0.0, 0.0], [0.0, 0.0], [0.0, 0.0], [0.05773502691896258, 0.0]]]]}]}], "bit_order": "MSB"})";
+
+bool validateDensityMatrix(
+    const std::vector<std::vector<std::complex<double>>> &in_dm) {
+  double sum_diag = 0.0;
+  const auto nbRows = in_dm.size();
+  for (size_t i = 0; i < nbRows; ++i) {
+    if (in_dm[i].size() != nbRows) {
+      return false;
+    }
+    const auto diag_elem = in_dm[i][i];
+    if (std::abs(diag_elem.imag()) > 1e-3) {
+      return false;
+    }
+    sum_diag += diag_elem.real();
+  }
+  return std::abs(1.0 - sum_diag) < 1e-3;
+}
+
+} // namespace
+
+TEST(JsonNoiseModelTester, checkNoNoise) {
+  {
+    auto accelerator =
+        xacc::getAccelerator("tnqvm", {{"tnqvm-visitor", "exatn-dm"}});
+    auto xasmCompiler = xacc::getCompiler("xasm");
+    auto program = xasmCompiler
+                       ->compile(R"(__qpu__ void testX(qbit q) {
+        X(q[0]);
+        Measure(q[0]);
+      })",
+                                 accelerator)
+                       ->getComposite("testX");
+    auto buffer = xacc::qalloc(1);
+    accelerator->execute(buffer, program);
+    auto exeInfo = accelerator->getExecutionInfo();
+    auto dm = accelerator
+                  ->getExecutionInfo<xacc::ExecutionInfo::DensityMatrixPtrType>(
+                      xacc::ExecutionInfo::DmKey);
+    std::cout << "Density matrix\n";
+    for (const auto &row : *dm) {
+      for (const auto &elem : row) {
+        std::cout << elem << " ";
+      }
+      std::cout << "\n";
+    }
+    EXPECT_EQ(dm->size(), 2);
+    EXPECT_TRUE(validateDensityMatrix(*dm));
+  }
+
+  {
+    auto accelerator =
+        xacc::getAccelerator("tnqvm", {{"tnqvm-visitor", "exatn-dm"}});
+    auto xasmCompiler = xacc::getCompiler("xasm");
+    auto program = xasmCompiler
+                       ->compile(R"(__qpu__ void testBell(qbit q) {
+        H(q[0]);
+        CX(q[0], q[1]);
+        Measure(q[0]);
+      })",
+                                 accelerator)
+                       ->getComposite("testBell");
+    auto buffer = xacc::qalloc(2);
+    accelerator->execute(buffer, program);
+    buffer->print();
+    auto exeInfo = accelerator->getExecutionInfo();
+    auto dm = accelerator
+                  ->getExecutionInfo<xacc::ExecutionInfo::DensityMatrixPtrType>(
+                      xacc::ExecutionInfo::DmKey);
+    // std::cout << "Density matrix\n";
+    // for (const auto &row : *dm) {
+    //   for (const auto &elem : row) {
+    //     std::cout << elem << " ";
+    //   }
+    //   std::cout << "\n";
+    // }
+    EXPECT_EQ(dm->size(), 4);
+    EXPECT_TRUE(validateDensityMatrix(*dm));
+  }
+
+  {
+    // Check exp-val by trace (direct tensor network contraction)
+    auto accelerator =
+        xacc::getAccelerator("tnqvm", {{"tnqvm-visitor", "exatn-dm"}});
+    auto xasmCompiler = xacc::getCompiler("xasm");
+    auto program = xasmCompiler
+                       ->compile(R"(__qpu__ void testRX(qbit q, double t) {
+        Rx(q[0], t);
+        Measure(q[0]);
+      })",
+                                 accelerator)
+                       ->getComposite("testRX");
+    const auto angles = xacc::linspace(-M_PI, M_PI, 20);
+    for (const auto &angle : angles) {
+      auto evaled = program->operator()({angle});
+      const double expResult =
+          -1.0 + 2.0 * std::cos(angle / 2.0) * std::cos(angle / 2.0);
+      std::cout << "Angle = " << angle << ": " << expResult << "\n";
+      auto buffer = xacc::qalloc(1);
+      accelerator->execute(buffer, evaled);
+      // buffer->print();
+      EXPECT_NEAR(buffer->getExpectationValueZ(), expResult, 1e-3);
+    }
+  }
+
+  {
+    // Check exp-val by trace (direct tensor network contraction)
+    auto accelerator =
+        xacc::getAccelerator("tnqvm", {{"tnqvm-visitor", "exatn-dm"}});
+    auto xasmCompiler = xacc::getCompiler("xasm");
+    auto program = xasmCompiler
+                       ->compile(R"(__qpu__ void testRY(qbit q, double t) {
+        Ry(q[0], t);
+        Measure(q[0]);
+      })",
+                                 accelerator)
+                       ->getComposite("testRY");
+    const auto angles = xacc::linspace(-M_PI, M_PI, 20);
+    for (const auto &angle : angles) {
+      auto evaled = program->operator()({angle});
+      const double expResult =
+          -1.0 + 2.0 * std::cos(angle / 2.0) * std::cos(angle / 2.0);
+      std::cout << "Angle = " << angle << ": " << expResult << "\n";
+      auto buffer = xacc::qalloc(1);
+      accelerator->execute(buffer, evaled);
+      // buffer->print();
+      EXPECT_NEAR(buffer->getExpectationValueZ(), expResult, 1e-3);
+    }
+  }
+
+  // Two-qubit expectation calculation by tracing DM tensor network
+  {
+    auto accelerator =
+        xacc::getAccelerator("tnqvm", {{"tnqvm-visitor", "exatn-dm"}});
+    auto xasmCompiler = xacc::getCompiler("xasm");
+    auto ir = xasmCompiler->compile(R"(__qpu__ void ansatz(qbit q, double t) {
+      X(q[0]);
+      Ry(q[1], t);
+      CX(q[1], q[0]);
+      H(q[0]);
+      H(q[1]);
+      Measure(q[0]);
+      Measure(q[1]);
+    })",
+                                    accelerator);
+
+    auto program = ir->getComposite("ansatz");
+    // Expected results from deuteron_2qbit_xasm_X0X1
+    const std::vector<double> expectedResults{
+        0.0,       -0.324699, -0.614213, -0.837166, -0.9694,
+        -0.996584, -0.915773, -0.735724, -0.475947, -0.164595,
+        0.164595,  0.475947,  0.735724,  0.915773,  0.996584,
+        0.9694,    0.837166,  0.614213,  0.324699,  0.0};
+
+    const auto angles =
+        xacc::linspace(-xacc::constants::pi, xacc::constants::pi, 20);
+    for (size_t i = 0; i < angles.size(); ++i) {
+      auto buffer = xacc::qalloc(2);
+      auto evaled = program->operator()({angles[i]});
+      accelerator->execute(buffer, evaled);
+      std::cout << "Angle = " << angles[i]
+                << ": Expected = " << expectedResults[i]
+                << "; ExaTN = " << buffer->getExpectationValueZ() << "\n";
+      EXPECT_NEAR(buffer->getExpectationValueZ(), expectedResults[i], 1e-6);
+    }
+  }
+}
+
+TEST(JsonNoiseModelTester, checkSimple) {
+  // Check depolarizing channels
+  {
+    auto noiseModel = xacc::getService<xacc::NoiseModel>("json");
+    noiseModel->initialize({{"noise-model", depol_json}});
+    auto accelerator = xacc::getAccelerator(
+        "tnqvm", {{"tnqvm-visitor", "exatn-dm"}, {"noise-model", noiseModel}});
+    auto xasmCompiler = xacc::getCompiler("xasm");
+    auto program = xasmCompiler
+                       ->compile(R"(__qpu__ void testX(qbit q) {
+        X(q[0]);
+        Measure(q[0]);
+      })",
+                                 accelerator)
+                       ->getComposite("testX");
+    auto buffer = xacc::qalloc(1);
+    accelerator->execute(buffer, program);
+    buffer->print();
+    auto dmPtr =
+        accelerator
+            ->getExecutionInfo<xacc::ExecutionInfo::DensityMatrixPtrType>(
+                xacc::ExecutionInfo::DmKey);
+    auto densityMatrix = *dmPtr;
+    EXPECT_TRUE(validateDensityMatrix(densityMatrix));
+    EXPECT_EQ(densityMatrix.size(), 2);
+    // Check trace
+    EXPECT_NEAR(densityMatrix[0][0].real() + densityMatrix[1][1].real(), 1.0,
+                1e-6);
+    // Expected result:
+    // 0.00666667+0.j 0.        +0.j
+    // 0.        +0.j 0.99333333+0.j
+    // Check real part
+    EXPECT_NEAR(densityMatrix[0][0].real(), 0.00666667, 1e-6);
+    EXPECT_NEAR(densityMatrix[1][0].real(), 0.0, 1e-6);
+    EXPECT_NEAR(densityMatrix[0][1].real(), 0.0, 1e-6);
+    EXPECT_NEAR(densityMatrix[1][1].real(), 0.99333333, 1e-6);
+  }
+
+  // Check amplitude damping channels
+  {
+    auto noiseModel = xacc::getService<xacc::NoiseModel>("json");
+    noiseModel->initialize({{"noise-model", ad_json}});
+    auto accelerator = xacc::getAccelerator(
+        "tnqvm", {{"tnqvm-visitor", "exatn-dm"}, {"noise-model", noiseModel}});
+    auto xasmCompiler = xacc::getCompiler("xasm");
+    auto program = xasmCompiler
+                       ->compile(R"(__qpu__ void testX_ad(qbit q) {
+        X(q[0]);
+        Measure(q[0]);
+      })",
+                                 accelerator)
+                       ->getComposites()[0];
+    auto buffer = xacc::qalloc(1);
+    accelerator->execute(buffer, program);
+    auto dmPtr =
+        accelerator
+            ->getExecutionInfo<xacc::ExecutionInfo::DensityMatrixPtrType>(
+                xacc::ExecutionInfo::DmKey);
+    auto densityMatrix = *dmPtr;
+    EXPECT_TRUE(validateDensityMatrix(densityMatrix));
+    // Verify the distribution (25% amplitude damping)
+    EXPECT_NEAR(densityMatrix[0][0].real(), 0.25, 1e-6);
+    EXPECT_NEAR(densityMatrix[1][1].real(), 0.75, 1e-6);
+  }
+}
+
+TEST(JsonNoiseModelTester, checkBitOrdering) {
+  auto xasmCompiler = xacc::getCompiler("xasm");
+  auto program = xasmCompiler
+                     ->compile(R"(__qpu__ void testCX(qbit q) {
+        CX(q[0], q[1]);
+        Measure(q[0]);
+        Measure(q[1]);
+      })",
+                               nullptr)
+                     ->getComposites()[0];
+
+  std::vector<std::vector<std::complex<double>>> densityMatrix_msb,
+      densityMatrix_lsb;
+  // Check MSB:
+  {
+    auto noiseModel = xacc::getService<xacc::NoiseModel>("json");
+    noiseModel->initialize({{"noise-model", msb_noise_model}});
+    auto accelerator = xacc::getAccelerator(
+        "tnqvm", {{"tnqvm-visitor", "exatn-dm"}, {"noise-model", noiseModel}});
+
+    auto buffer = xacc::qalloc(2);
+    accelerator->execute(buffer, program);
+    auto dmPtr =
+        accelerator
+            ->getExecutionInfo<xacc::ExecutionInfo::DensityMatrixPtrType>(
+                xacc::ExecutionInfo::DmKey);
+    densityMatrix_msb = *dmPtr;
+    EXPECT_TRUE(validateDensityMatrix(densityMatrix_msb));
+    std::cout << "Density matrix MSB\n";
+    for (const auto &row : densityMatrix_msb) {
+      for (const auto &elem : row) {
+        std::cout << elem << " ";
+      }
+      std::cout << "\n";
+    }
+  }
+
+  // Check LSB:
+  {
+    auto noiseModel = xacc::getService<xacc::NoiseModel>("json");
+    noiseModel->initialize({{"noise-model", lsb_noise_model}});
+    auto accelerator = xacc::getAccelerator(
+        "tnqvm", {{"tnqvm-visitor", "exatn-dm"}, {"noise-model", noiseModel}});
+    auto buffer = xacc::qalloc(2);
+    accelerator->execute(buffer, program);
+    auto dmPtr =
+        accelerator
+            ->getExecutionInfo<xacc::ExecutionInfo::DensityMatrixPtrType>(
+                xacc::ExecutionInfo::DmKey);
+    densityMatrix_lsb = *dmPtr;
+    EXPECT_TRUE(validateDensityMatrix(densityMatrix_lsb));
+    std::cout << "Density matrix LSB\n";
+    for (const auto &row : densityMatrix_lsb) {
+      for (const auto &elem : row) {
+        std::cout << elem << " ";
+      }
+      std::cout << "\n";
+    }
+  }
+
+  for (int row = 0; row < 4; ++row) {
+    for (int col = 0; col < 4; ++col) {
+      EXPECT_NEAR(densityMatrix_msb[row][col].real(),
+                  densityMatrix_lsb[row][col].real(), 1e-9);
+      EXPECT_NEAR(densityMatrix_msb[row][col].imag(),
+                  densityMatrix_lsb[row][col].imag(), 1e-9);
+    }
+  }
+}
+
+TEST(JsonNoiseModelTester, checkMultiDepol) {
+  auto xasmCompiler = xacc::getCompiler("xasm");
+  auto program = xasmCompiler
+                     ->compile(R"(__qpu__ void testCX2(qbit q) {
+        CX(q[0], q[1]);
+        Measure(q[0]);
+        Measure(q[1]);
+      })",
+                               nullptr)
+                     ->getComposites()[0];
+
+  auto noiseModel = xacc::getService<xacc::NoiseModel>("json");
+  noiseModel->initialize({{"noise-model", noise_model_multi_depol}});
+  auto accelerator = xacc::getAccelerator(
+      "tnqvm", {{"tnqvm-visitor", "exatn-dm"}, {"noise-model", noiseModel}});
+
+  auto buffer = xacc::qalloc(2);
+  accelerator->execute(buffer, program);
+  auto dmPtr =
+      accelerator->getExecutionInfo<xacc::ExecutionInfo::DensityMatrixPtrType>(
+          xacc::ExecutionInfo::DmKey);
+  auto densityMatrix = *dmPtr;
+  EXPECT_TRUE(validateDensityMatrix(densityMatrix));
+  std::cout << "Density matrix:\n";
+  // for (const auto &row : densityMatrix) {
+  //   for (const auto &elem : row) {
+  //     std::cout << elem << " ";
+  //   }
+  //   std::cout << "\n";
+  // }
+}
+
+// Test the ordering of root tensor legs after contraction.
+TEST(JsonNoiseModelTester, checkOutputLegOrder) {
+  std::vector<std::vector<std::complex<double>>> densityMatrix_no_noise,
+      densityMatrix_weak_noise;
+  auto xasmCompiler = xacc::getCompiler("xasm");
+  auto program = xasmCompiler
+                     ->compile(R"(__qpu__ void testLegOrdering(qbit q) {
+        H(q[0]);
+      })",
+                               nullptr)
+                     ->getComposites()[0];
+
+  {
+    // No noise
+    auto accelerator =
+        xacc::getAccelerator("tnqvm", {{"tnqvm-visitor", "exatn-dm"}});
+
+    auto buffer = xacc::qalloc(4);
+    accelerator->execute(buffer, program);
+    auto dmPtr =
+        accelerator
+            ->getExecutionInfo<xacc::ExecutionInfo::DensityMatrixPtrType>(
+                xacc::ExecutionInfo::DmKey);
+    densityMatrix_no_noise = *dmPtr;
+    EXPECT_TRUE(validateDensityMatrix(densityMatrix_no_noise));
+    std::cout << "Density matrix:\n";
+    for (const auto &row : densityMatrix_no_noise) {
+      for (const auto &elem : row) {
+        std::cout << elem << " ";
+      }
+      std::cout << "\n";
+    }
+  }
+  {
+    // Extremely weak noise:
+    auto noiseModel = xacc::getService<xacc::NoiseModel>("json");
+    noiseModel->initialize({{"noise-model", noise_model_weak_noise}});
+    auto accelerator = xacc::getAccelerator(
+        "tnqvm", {{"tnqvm-visitor", "exatn-dm"}, {"noise-model", noiseModel}});
+
+    auto buffer = xacc::qalloc(4);
+    accelerator->execute(buffer, program);
+    auto dmPtr =
+        accelerator
+            ->getExecutionInfo<xacc::ExecutionInfo::DensityMatrixPtrType>(
+                xacc::ExecutionInfo::DmKey);
+    densityMatrix_weak_noise = *dmPtr;
+    EXPECT_TRUE(validateDensityMatrix(densityMatrix_weak_noise));
+    std::cout << "Density matrix:\n";
+    for (const auto &row : densityMatrix_weak_noise) {
+      for (const auto &elem : row) {
+        std::cout << elem << " ";
+      }
+      std::cout << "\n";
+    }
+  }
+
+  // 4 qubits
+  const auto nbRows = 1ULL << 4;
+  for (int row = 0; row < nbRows; ++row) {
+    for (int col = 0; col < nbRows; ++col) {
+      EXPECT_NEAR(densityMatrix_no_noise[row][col].real(),
+                  densityMatrix_weak_noise[row][col].real(), 1e-3);
+      EXPECT_NEAR(densityMatrix_no_noise[row][col].imag(),
+                  densityMatrix_weak_noise[row][col].imag(), 1e-3);
+    }
+  }
+}
+
+TEST(JsonNoiseModelTester, checkGHZ) {
+  auto xasmCompiler = xacc::getCompiler("xasm");
+  auto program = xasmCompiler
+                     ->compile(R"(__qpu__ void testGHZ(qbit q) {
+        H(q[0]);
+        CX(q[0], q[1]);
+        CX(q[1], q[2]);
+        CX(q[2], q[3]);
+        Measure(q[0]);
+        Measure(q[1]);
+        Measure(q[2]);
+        Measure(q[3]);
+      })",
+                               nullptr)
+                     ->getComposites()[0];
+
+  auto noiseModel = xacc::getService<xacc::NoiseModel>("json");
+  noiseModel->initialize({{"noise-model", noise_model_4q}});
+  auto accelerator = xacc::getAccelerator(
+      "tnqvm", {{"tnqvm-visitor", "exatn-dm"}, {"noise-model", noiseModel}});
+
+  auto buffer = xacc::qalloc(4);
+  accelerator->execute(buffer, program);
+  auto dmPtr =
+      accelerator->getExecutionInfo<xacc::ExecutionInfo::DensityMatrixPtrType>(
+          xacc::ExecutionInfo::DmKey);
+  auto densityMatrix = *dmPtr;
+  EXPECT_TRUE(validateDensityMatrix(densityMatrix));
+  std::cout << "Density matrix:\n";
+  for (const auto &row : densityMatrix) {
+    for (const auto &elem : row) {
+      std::cout << elem << " ";
+    }
+    std::cout << "\n";
+  }
+
+  // Expect the first and last element of the matrix are near 0.5 (GHZ state)
+  const auto dmDimension = 1ULL << buffer->size();
+  EXPECT_NEAR(densityMatrix[0][0].real(), 0.5, 0.1);
+  EXPECT_NEAR(densityMatrix[dmDimension - 1][dmDimension - 1].real(), 0.5, 0.1);
+  EXPECT_NEAR(densityMatrix[0][0].imag(), 0.0, 1e-6);
+  EXPECT_NEAR(densityMatrix[dmDimension - 1][dmDimension - 1].imag(), 0.0,
+              1e-6);
+}
+
+// Test VQE with noise
+TEST(JsonNoiseModelTester, testDeuteronVqeH2) {
+  auto noiseModel = xacc::getService<xacc::NoiseModel>("json");
+  noiseModel->initialize({{"noise-model", noise_model_4q}});
+  auto accelerator = xacc::getAccelerator(
+      "tnqvm", {{"tnqvm-visitor", "exatn-dm"}, {"noise-model", noiseModel}});
+
+  // Create the N=2 deuteron Hamiltonian
+  auto H_N_2 = xacc::quantum::getObservable(
+      "pauli", std::string("5.907 - 2.1433 X0X1 "
+                           "- 2.1433 Y0Y1"
+                           "+ .21829 Z0 - 6.125 Z1"));
+
+  auto optimizer = xacc::getOptimizer("nlopt");
+  xacc::qasm(R"(
+        .compiler xasm
+        .circuit deuteron_ansatz1
+        .parameters theta
+        .qbit q
+        X(q[0]);
+        Ry(q[1], theta);
+        CNOT(q[1],q[0]);
+    )");
+  auto ansatz = xacc::getCompiled("deuteron_ansatz1");
+
+  // Get the VQE Algorithm and initialize it
+  auto vqe = xacc::getAlgorithm("vqe");
+  vqe->initialize({std::make_pair("ansatz", ansatz),
+                   std::make_pair("observable", H_N_2),
+                   std::make_pair("accelerator", accelerator),
+                   std::make_pair("optimizer", optimizer)});
+
+  // Allocate some qubits and execute
+  auto buffer = xacc::qalloc(2);
+  vqe->execute(buffer);
+  std::cout << "Energy = " << (*buffer)["opt-val"].as<double>() << "\n";
+  // Expected result: -1.74886
+  // Relax the limit due to noise.
+  EXPECT_NEAR((*buffer)["opt-val"].as<double>(), -1.74886, 0.5);
+}
+
+int main(int argc, char **argv) {
+  xacc::Initialize(argc, argv);
+  ::testing::InitGoogleTest(&argc, argv);
+  auto ret = RUN_ALL_TESTS();
+  xacc::Finalize();
+  return ret;
+}
\ No newline at end of file