Barthel research group,  Duke Physics condensed matter theory and numerics

• Dynamic response of strongly-correlated quantum systems
• Tensor network state methods for strongly-correlated systems in D≥2 dimensions
• Investigating quantum matter with entanglement renormalization on quantum computers
• Scaling and evolution of entanglement in many-body systems
• Nonequilibrium dynamics and driven-dissipative systems
• Stochastic dynamics in network systems and epidemic outbreaks

Dynamic response. - The response of a quantum system to a perturbation allows us to probe its low-energy physics, in particular, the nature and spectrum of excitations. Typical examples of response functions are spectral functions in electronic systems and dynamic structure factors in spin systems. They are the natural observables in scattering experiments. We have worked on tensor network state techniques for the computation of precise response functions and apply them to investigate quantum magnets and ultracold atom systems.

Method development. - In [1], we showed how to compute dynamic structure factors of strongly correlated systems at finite temperatures using matrix-product-state (MPS) purifications and linear prediction. Due to the increase of entanglement with time, the precise simulations are usually limited to finite time windows. Linear prediction helps to extrapolate to longer times and substantially improve spectral resolution. Furthermore, the maximum reachable times can be increased by at least a factor of two by rearranging time evolution steps to reduce entanglement growth in the purifications [2,3]. We also showed how infinite boundary conditions, i.e., a simulation in the thermodynamic limit using infinite MPS, can be used to reduce the number of time evolution runs to just two [4,5]. Lastly, as an alternative to purifications, one can employ a stochastic approach called METTS [White 2009], which has advantages for 2D systems [6,7]. Generally, the big advantage of these techniques is that they work directly in the real-time/real-frequency domain and hence do not require the (usually uncontrolled) analytical continuation employed in corresponding quantum Monte Carlo simulations [1].

Applications. - We have applied these techniques to several systems. For example, we worked with experimentalists of the Helmholtz-Zentrum Berlin and Oak Ridge National Lab to study thermal spin structure factors of quasi-1D spin-1/2 magnets as obtained by neutron scattering [8]. We have also determined the universal scaling function for bosonic spectral functions in the z=2 quantum critical regime [2]. In [5], we investigated bilinear-biquadratic spin-1 chains, which describe quantum magnets like CsNiCl₃, NENP, or LiVGe₂O₆. The Haldane phase, which features symmetry-protected topological order, is usually described by the O(3) nonlinear sigma model. However, we found qualitative discrepancies of this approximate field theoretical description. The evolution of multisoliton continua in the critical phase of the spin-1 chains can be understood qualitatively through the renormalization group analysis of the corresponding level-one SU(3) Wess-Zumino-Witten (WZW) model. Numerically, we found previously unnoticed stable excitations above the contracting multisoliton continua.

		[1] "Spectral functions in one-dimensional quantum systems at finite temperature using the density matrix renormalization group" T. Barthel, U. Schollwöck, and S. R. White arXiv:0901.2342, pdf, Phys. Rev. B 79, 245101 (2009) We present for the first time time-dependent density-matrix renormalization-group simulations (t-DMRG) at finite temperatures. It is demonstrated how a combination of finite-temperature t-DMRG and time-series prediction allows for an easy and very accurate calculation of spectral functions in one-dimensional quantum systems, irrespective of their statistics, for arbitrary temperatures. This is illustrated with spin structure factors of XX and XXX spin-1/2 chains. For the XX model we can compare against an exact solution and for the XXX model (Heisenberg antiferromagnet) against a Bethe Ansatz solution and quantum Monte Carlo data.
		[2] "Scaling of the thermal spectral function for quantum critical bosons in one dimension" T. Barthel, U. Schollwöck, and S. Sachdev arXiv:1212.3570, pdf We present an improved scheme for the precise evaluation of finite-temperature response functions of strongly correlated systems in the framework of the time-dependent density matrix renormalization group. The maximum times that we can reach at finite temperatures T are typically increased by a factor of two, when compared against the earlier approaches. This novel scheme, complemented with linear prediction, allows us now to evaluate dynamic correlators for interacting bosons in one dimension. We demonstrate that the considered spectral function in the quantum critical regime with dynamic critical exponent z=2 is captured by the universal scaling form S(k,omega)=(1/T)*Phi(k/sqrt(T),omega/T) and calculate the scaling function precisely.
		[3] "Precise evaluation of thermal response functions by optimized density matrix renormalization group schemes" T. Barthel arXiv:1301.2246, pdf, New J. Phys. 15, 073010 (2013) This paper provides a study and discussion of earlier as well as novel more efficient schemes for the precise evaluation of finite-temperature response functions of strongly correlated quantum systems in the framework of the time-dependent density matrix renormalization group (tDMRG). The computational costs and bond dimensions as functions of time and temperature are examined at the example of the spin-1/2 XXZ Heisenberg chain in the critical XY phase and the gapped Néel phase. The matrix product state purifications occurring in the algorithms are in one-to-one relation with corresponding matrix product operators. This notational simplification elucidates implications of quasi-locality on the computational costs. Based on the observation that there is considerable freedom in designing efficient tDMRG schemes for the calculation of dynamical correlators at finite temperatures, a new class of optimizable schemes, as recently suggested in arXiv:1212.3570, is explained and analyzed numerically. A specific novel near-optimal scheme that requires no additional optimization reaches maximum times that are typically increased by a factor of two, when compared against earlier approaches. These increased reachable times make many more physical applications accessible. For each of the described tDMRG schemes, one can devise a corresponding transfer matrix renormalization group (TMRG) variant.
		[4] "Infinite boundary conditions for response functions and limit cycles in iDMRG, demonstrated for bilinear-biquadratic spin-1 chains" M. Binder and T. Barthel arXiv:1804.09163, pdf, Phys. Rev. B 98, 235114 (2018) Response functions <A(x,t)B(y,0)> for one-dimensional strongly correlated quantum many-body systems can be computed with matrix product state (MPS) techniques. Especially, when one is interested in spectral functions or dynamic structure factors of translation-invariant systems, the response for some range |x-y|<L is needed. We demonstrate how the number of required time-evolution runs can be reduced substantially: (a) If finite-system simulations are employed, the number of time-evolution runs can be reduced from L to 2sqrt(L). (b) To go beyond, one can employ infinite MPS (iMPS) such that two evolution runs suffice. To this purpose, iMPS that are heterogeneous only around the causal cone of the perturbation are evolved in time, i.e., the simulation is done with so-called infinite boundary conditions. Computing overlaps of these states, spatially shifted relative to each other, yields the response functions for all distances |x-y|. As a specific application, we compute the dynamic structure factor for ground states of bilinear-biquadratic spin-1 chains with very high resolution and explain the underlying low-energy physics. To determine the initial uniform iMPS for such simulations, infinite-system density-matrix renormalization group (iDMRG) can be employed. We discuss that, depending on the system and chosen bond dimension, iDMRG with a cell size n may converge to a non-trivial limit cycle of length m. This then corresponds to an iMPS with an enlarged unit cell of size mn.
		[5] "Low-energy physics of isotropic spin-1 chains in the critical and Haldane phases" M. Binder and T. Barthel arXiv:2005.03643, pdf, Phys. Rev. B 102, 014447 (2020) Using a matrix product state algorithm with infinite boundary conditions, we compute high-resolution dynamic spin and quadrupolar structure factors in the thermodynamic limit to explore the low-energy excitations of isotropic bilinear-biquadratic spin-1 chains. Haldane mapped the spin-1 Heisenberg antiferromagnet to a continuum field theory, the nonlinear sigma model (NLSM). We find that the NLSM fails to capture the influence of the biquadratic term and provides only an unsatisfactory description of the Haldane phase physics. But several features in the Haldane phase can be explained by noninteracting multimagnon states. The physics at the Uimin-Lai-Sutherland point is characterized by multisoliton continua. Moving into the extended critical phase, we find that these excitation continua contract, which we explain using a field-theoretic description. New excitations emerge at higher energies and, in the vicinity of the purely biquadratic point, they show simple cosine dispersions. Using block fidelities, we identify them as elementary one-particle excitations and relate them to the integrable Temperley-Lieb chain.
		[6] "Minimally entangled typical thermal states versus matrix product purifications for the simulation of equilibrium states and time evolution" M. Binder and T. Barthel arXiv:1411.3033, pdf, Phys. Rev. B 92, 125119 (2015) For the simulation of equilibrium states and finite-temperature response functions of strongly-correlated quantum many-body systems, we compare the efficiencies of two different approaches in the framework of the density matrix renormalization group (DMRG). The first is based on matrix product purifications. The second, more recent one, is based on so-called minimally entangled typical thermal states (METTS). For the latter, we highlight the interplay of statistical and DMRG truncation errors, discuss the use of self-averaging effects, and describe schemes for the computation of response functions. For critical as well as gapped phases of the spin-1/2 XXZ chain and the one-dimensional Bose-Hubbard model, we assess computation costs and accuracies of the two methods at different temperatures. For almost all considered cases, we find that, for the same computation cost, purifications yield more accurate results than METTS -- often by orders of magnitude. The METTS algorithm becomes more efficient only for temperatures well below the system's energy gap. The exponential growth of the computation cost in the evaluation of response functions limits the attainable timescales in both methods and we find that in this regard, METTS do not outperform purifications.
		[7] "Symmetric minimally entangled typical thermal states for canonical and grand-canonical ensembles" M. Binder and T. Barthel arXiv:1701.03872, pdf, Phys. Rev. B 95, 195148 (2017) Based on the density matrix renormalization group (DMRG), strongly correlated quantum many-body systems at finite temperatures can be simulated by sampling over a certain class of pure matrix product states (MPS) called minimally entangled typical thermal states (METTS). When a system features symmetries, these can be utilized to substantially reduce MPS computation costs. It is conceptually straightforward to simulate canonical ensembles using symmetric METTS. In practice, it is important to alternate between different symmetric collapse bases to decrease autocorrelations in the Markov chain of METTS. To this purpose, we introduce symmetric Fourier and Haar-random block bases that are efficiently mixing. We also show how grand-canonical ensembles can be simulated efficiently with symmetric METTS. We demonstrate these approaches for spin-1/2 XXZ chains and discuss how the choice of the collapse bases influences autocorrelations as well as the distribution of measurement values and, hence, convergence speeds.
		[8] "Multispinon continua at zero and finite temperature in a near-ideal Heisenberg chain" B. Lake, D. A. Tennant, J.-S. Caux, T. Barthel, U. Schollwöck, S. E. Nagler, and C. D. Frost arXiv:1307.4071, pdf, Phys. Rev. Lett. 111, 137205 (2013) The space- and time-dependent response of many-body quantum systems is the most informative aspect of their emergent behaviour. The dynamical structure factor, experimentally measurable using neutron scattering, can map this response in wavevector and energy with great detail, allowing theories to be quantitatively tested to high accuracy. Here, we present a comparison between neutron scattering measurements on the one-dimensional spin-1/2 Heisenberg antiferromagnet KCuF3, and recent state-of-the-art theoretical methods based on integrability and density matrix renormalization group simulations. The unprecedented quantitative agreement shows that precise descriptions of strongly correlated states at all distance, time and temperature scales are now possible, and highlights the need to apply these novel techniques to other problems in low-dimensional magnetism.

Strong correlations. - Due to their complexity, several systems of high technological relevance like high-temperature superconductors, frustrated magnets, and systems with topological order are still unsatisfactorily understood. This is due to the fact that the Hilbert space dimension grows exponentially in the system size and that perturbative, mean-field, as well as effective single-particle methods (like DFT) fail for strongly correlated systems. Hope comes from the analysis of entanglement properties, which shows that nature typically only explores a small corner of the exponentially big Hilbert space, and we may hope to capture it with an appropriately chosen reduced set of effective degrees of freedom.

Tensor network states. - We are working on tensor network state (TNS) methods, which have become a dominant numerical tool for quantum many-body physics. TNS are designed such that the number of effective degrees of freedom can be tuned (and hence accuracy and computation costs) but does not need to grow exponentially with the system size. In particular, the many-body quantum state is approximated by a network of partially contracted tensors, with bond dimensions of these tensors controlling the accuracy.

Development of TNS techniques. - Some of our contributions concern the efficient TNS-based evaluation of dynamic response functions. We developed a generic and efficient algorithm that makes it possible to apply TNS to fermionic systems in D≥2 dimensions with only marginal computational overhead by avoiding a global Jordan-Wigner transformation; see [1,2] and references therein. This gives TNS a major advantage for the investigation of fermionic systems in D≥2 dimensions and frustrated magnets because it avoids the negative-sign problem that hampers quantum Monte Carlo for such systems. TNS have been applied in countless studies and produced important new insights. However, due to the scaling of computation costs in the bond dimension m (e.g., O(m¹²) for 2D PEPS), the overwhelming majority of studies address one-dimensional (1D) and quasi-1D systems, and practicable m are rather small for systems in D≥2 dimensions. We are pursuing different routes to substantially reduce TNS computation costs while maintaining expressiveness: (a) In [3], we found that this can be achieved by imposing constraints on the canonical polyadic (CP) rank of the tensors. We have also adapted the approach for machine learning tasks, optimizing tree tensor network classifiers with rank-constrained tensors for image classification [4]. The CP rank constraints make it possible to work with networks of high vertex degree. (b) In [5], we investigated how simulations based on the multi-scale entanglement renormalization ansatz (MERA) [Vidal 2007] can be made more efficient through tensor Trotterization (choosing the tensors as circuits of two-qubit gates) and/or a stochastic evaluation of energy gradients (variational Monte Carlo). (c) With the goal of realistically capturing microscopic material details, we are working on a new TNS impurity solver for dynamical mean-field theory (DMFT). DMFT can be interpreted as an extension of mean-field theory, where the mean field is now a dynamical single-particle Green's function. A translation-invariant system is split into a small impurity (a cluster), and the rest of the lattice is interpreted as the impurity's environment. The full-lattice problem is reduced to a self-consistent impurity problem, where interactions are restricted to the impurity, and the environment is modeled by a non-interacting bath.

		[1] "Contraction of fermionic operator circuits and the simulation of strongly correlated fermions" T. Barthel, C. Pineda, and J. Eisert arXiv:0907.3689, pdf, Phys. Rev. A 80, 042333 (2009), also in Virtual Journal of Nanoscale Science and Technology 20, Issue 20 (2009), and in Virtual Journal of Quantum Information 9, Issue 11 (2009) A fermionic operator circuit is a product of fermionic operators of usually different and partially overlapping support. Further elements of fermionic operator circuits (FOCs) are partial traces and partial projections. The presented framework allows for the introduction of fermionic versions of known qudit operator circuits (QUOC), important for the simulation of strongly correlated d-dimensional systems: the multiscale entanglement renormalization ansaetze (MERA), tree tensor networks (TTN), projected entangled pair states (PEPS), or their infinite-size versions (iPEPS etc.). After the definition of a FOC, we present a method to contract it with the same computation and memory requirements as a corresponding QUOC, for which all fermionic operators are replaced by qudit operators of identical dimension. A given scheme for contracting the QUOC relates to an analogous scheme for the corresponding fermionic circuit, where additional marginal computational costs arise only from reordering of modes for operators occurring in intermediate stages of the contraction. Our result hence generalizes efficient schemes for the simulation of d-dimensional spin systems, as MERA, TTN, or PEPS to the fermionic case.
		[2] "Unitary circuits for strongly correlated fermions" C. Pineda, T. Barthel, and J. Eisert arXiv:0905.0669, pdf, Phys. Rev. A 81, 050303(R) (2010) We introduce a scheme for efficiently describing pure states of strongly correlated fermions in higher dimensions using unitary circuits. A local way of computing local expectation values is presented. We formulate a dynamical reordering scheme, corresponding to time-adaptive Jordan-Wigner transformation that avoids non-local string operators and only keeps suitably ordered the causal cone. Primitives of such a reordering scheme are highlighted. Fermionic unitary circuits can be contracted with the same complexity as in the spin case. The scheme gives rise to a variational description of fermionic models that does not suffer from a sign problem. We present a numerical example on a 9x9 fermionic lattice model to show the functioning of the approach.
		[3] "Tensor network states with low-rank tensors" H. Chen and T. Barthel arXiv:2205.15296, pdf Tensor networks are used to efficiently approximate states of strongly-correlated quantum many-body systems. More generally, tensor network approximations may allow to reduce the costs for operating on an order-N tensor from exponential to polynomial in N, and this has become a popular approach for machine learning. We introduce the idea of imposing constraints on the canonical polyadic rank of the tensors that compose the tensor network. With this modification, the time and space complexities for the network optimization can be substantially reduced while maintaining high accuracy. We detail this idea for tree tensor network states (TTNS) and projected entangled-pair states. Simulations of spin models on Cayley trees with low-rank TTNS exemplify the effect of rank constraints on the expressive power. We find that choosing the tensor rank r to be on the order of the bond dimension m, is sufficient to obtain high-accuracy groundstate approximations and to substantially outperform standard TTNS computations. Thus low-rank tensor networks are a promising route for the simulation of quantum matter and machine learning on large data sets.
		[4] "Machine learning with tree tensor networks, CP rank constraints, and tensor dropout" H. Chen and T. Barthel arXiv:2305.19440, pdf, IEEE Trans. Pattern Anal. Mach. Intell. 1 (2024) Tensor networks approximate order-N tensors with a reduced number of degrees of freedom that is only polynomial in N and arranged as a network of partially contracted smaller tensors. As suggested in [arXiv:2205.15296] in the context of quantum many-body physics, computation costs can be further substantially reduced by imposing constraints on the canonical polyadic (CP) rank of the tensors in such networks. Here we demonstrate how tree tensor networks (TTN) with CP rank constraints and tensor dropout can be used in machine learning. The approach is found to outperform other tensor-network based methods in Fashion-MNIST image classification. A low-rank TTN classifier with branching ratio 4 reaches test set accuracy of 90.3 percent with low computation costs. Consisting of mostly linear elements, tensor network classifiers avoid the vanishing gradient problem of deep neural networks. The CP rank constraints have additional advantages: The number of parameters can be decreased and tuned more freely to control overfitting, improve generalization properties, and reduce computation costs. They allow us to employ trees with large branching ratios which substantially improves the representation power.
		[5] "Scaling of contraction costs for entanglement renormalization algorithms including tensor Trotterization and variational Monte Carlo" T. Barthel and Q. Miao arXiv:2407.21006, pdf The multi-scale entanglement renormalization ansatz (MERA) is a hierarchical class of tensor network states motivated by the real-space renormalization group. It is used to simulate strongly correlated quantum many-body systems. For prominent MERA structures in one and two spatial dimensions, we determine the optimal scaling of contraction costs as well as corresponding contraction sequences and algorithmic phase diagrams. This is motivated by recent efforts to employ MERA in hybrid quantum-classical algorithms, where the MERA tensors are Trotterized, i.e., chosen as circuits of quantum gates, and observables as well as energy gradients are evaluated by sampling causal-cone states. We investigate whether tensor Trotterization and/or variational Monte Carlo (VMC) sampling can lead to quantum-inspired classical MERA algorithms that perform better than the traditional optimization of full MERA based on the exact evaluation of energy gradients. Algorithmic phase diagrams indicate the best MERA method depending on the scaling of the energy accuracy and the number of Trotter steps with the bond dimension. The results suggest substantial gains due to VMC for two-dimensional systems.

VQA with Trotterized MERA. - So far, we have strong evidence that equilibrium properties of most strongly correlated systems of practical interest can be simulated efficiently using tensor network states (TNS). However, TNS simulations in D≥2 dimensions are quite expensive. On classical computers, the computational costly steps are tensor contractions which scale with a high power in the TNS bond dimension. Recently, we have developed a resource-efficient and noise-resilient variational quantum-classical algorithm (VQA) for the investigation of strongly-correlated quantum matter on quantum computers [1]. It is based on multi-scale entanglement renormalization tensor networks (MERA) and gradient-based optimization. Due to the narrow causal cones of MERA, the algorithm can be implemented on noisy intermediate-scale (NISQ) devices with a relatively small number of qubits and still describe very large systems. The number of required qubits is system-size independent. Translation invariance can be used to make computation costs square-logarithmic in the system size and describe the thermodynamic limit. The basic idea is to leverage the power of quantum computers to implement the MERA tensor contractions. To this purpose all tensors are Trotterized, i.e., realized as brickwall circuits of two-qubit gates [2]. We have established a polynomial quantum advantage for different critical spin models [3] and are working on various improvements and experimental demonstrations.

Absence of barren plateaus for isometric TNS. - A common obstacle for VQA are so-called barren plateaus, where the average energy-gradient amplitude decreases exponentially in the system size. Barren plateaus are prevalent in various VQA like quantum neural networks. VQA with barren plateaus are not trainable as the inability to precisely estimate exponentially small gradients will result in random walks on a basically flat energy landscape. In the recent contributions [4,5], we proved that the variational optimizations of isometric TNS (MPS, TTNS, and MERA) for extensive Hamiltonians with finite-range interactions are free of barren plateaus. The variance of the energy gradient, evaluated by taking the Haar average over the TNS tensors, has a leading system-size independent term and decreases according to a power law in the bond dimension. For a hierarchical TNS (TTNS or MERA), the variance of the gradient decays exponentially in the layer index of the considered tensor. These results also motivate certain improved initialization schemes. Lastly, we found a direct connection between gradient vanishing and cost-function concentration for quantum circuits [6].

		[1] "Quantum-classical eigensolver using multiscale entanglement renormalization" Q. Miao and T. Barthel arXiv:2108.13401, pdf, Phys. Rev. Research 5, 033141 (2023) We propose a variational quantum eigensolver (VQE) for the simulation of strongly correlated quantum matter based on a multiscale entanglement renormalization ansatz (MERA) and gradient-based optimization. This MERA quantum eigensolver can have substantially lower computation costs than corresponding classical algorithms. Due to its narrow causal cone, the algorithm can be implemented on noisy intermediate-scale quantum (NISQ) devices and still describe large systems. It is particularly attractive for ion-trap devices with ion-shuttling capabilities. The number of required qubits is system-size independent and increases only to a logarithmic scaling when using quantum amplitude estimation to speed up gradient evaluations. Translation invariance can be used to make computation costs square-logarithmic in the system size and describe the thermodynamic limit. We demonstrate the approach numerically for a MERA with Trotterized disentanglers and isometries. With a few Trotter steps, one recovers the accuracy of the full MERA.
		[2] "Scaling of contraction costs for entanglement renormalization algorithms including tensor Trotterization and variational Monte Carlo" T. Barthel and Q. Miao arXiv:2407.21006, pdf The multi-scale entanglement renormalization ansatz (MERA) is a hierarchical class of tensor network states motivated by the real-space renormalization group. It is used to simulate strongly correlated quantum many-body systems. For prominent MERA structures in one and two spatial dimensions, we determine the optimal scaling of contraction costs as well as corresponding contraction sequences and algorithmic phase diagrams. This is motivated by recent efforts to employ MERA in hybrid quantum-classical algorithms, where the MERA tensors are Trotterized, i.e., chosen as circuits of quantum gates, and observables as well as energy gradients are evaluated by sampling causal-cone states. We investigate whether tensor Trotterization and/or variational Monte Carlo (VMC) sampling can lead to quantum-inspired classical MERA algorithms that perform better than the traditional optimization of full MERA based on the exact evaluation of energy gradients. Algorithmic phase diagrams indicate the best MERA method depending on the scaling of the energy accuracy and the number of Trotter steps with the bond dimension. The results suggest substantial gains due to VMC for two-dimensional systems.
		[3] "Convergence and quantum advantage of Trotterized MERA for strongly-correlated systems" Q. Miao and T. Barthel arXiv:2303.08910, pdf Strongly-correlated quantum many-body systems are difficult to study and simulate classically. We recently proposed a variational quantum eigensolver (VQE) based on the multiscale entanglement renormalization ansatz (MERA) with tensors constrained to certain Trotter circuits. Here, we extend the theoretical analysis, testing different initialization and convergence schemes, determining the scaling of computation costs for various critical spin models, and establishing a quantum advantage. For the Trotter circuits being composed of single-qubit and two-qubit rotations, it is experimentally advantageous to have small rotation angles. We find that the average angle amplitude can be reduced substantially with negligible effect on the energy accuracy. Benchmark simulations show that choosing TMERA tensors as brick-wall circuits or parallel random-pair circuits yields very similar energy accuracies.
		[4] "Absence of barren plateaus and scaling of gradients in the energy optimization of isometric tensor network states" T. Barthel and Q. Miao arXiv:2304.00161, pdf Vanishing gradients can pose substantial obstacles for high-dimensional optimization problems. Here we consider energy minimization problems for quantum many-body systems with extensive Hamiltonians and finite-range interactions, which can be studied on classical computers or in the form of variational quantum eigensolvers on quantum computers. Barren plateaus correspond to scenarios where the average amplitude of the energy gradient decreases exponentially with increasing system size. This occurs, for example, for quantum neural networks and for brickwall quantum circuits when the depth increases polynomially in the system size. Here we prove that the variational optimization problems for matrix product states, tree tensor networks, and the multiscale entanglement renormalization ansatz are free of barren plateaus. The derived scaling properties for the gradient variance provide an analytical guarantee for the trainability of randomly initialized tensor network states (TNS) and motivate certain initialization schemes. In a suitable representation, unitary tensors that parametrize the TNS are sampled according to the uniform Haar measure. We employ a Riemannian formulation of the gradient based optimizations which simplifies the analytical evaluation.
		[5] "Isometric tensor network optimization for extensive Hamiltonians is free of barren plateaus" Q. Miao and T. Barthel arXiv:2304.14320, pdf, Phys. Rev. A 109, L050402 (2024) We explain why and numerically confirm that there are no barren plateaus in the energy optimization of isometric tensor network states (TNS) for extensive Hamiltonians with finite-range interactions. Specifically, we consider matrix product states, tree tensor network states, and the multiscale entanglement renormalization ansatz. The variance of the energy gradient, evaluated by taking the Haar average over the TNS tensors, has a leading system-size independent term and decreases according to a power law in the bond dimension. For a hierarchical TNS with branching ratio b, the variance of the gradient with respect to a tensor in layer t scales as (br)^t, where r is the second largest eigenvalue of the Haar-average doubled layer-transition channel and decreases algebraically with increasing bond dimension. The observed scaling properties of the gradient variance bear implications for efficient initialization procedures.
		[6] "Equivalence of cost concentration and gradient vanishing for quantum circuits: An elementary proof in the Riemannian formulation" Q. Miao and T. Barthel arXiv:2402.07883, pdf, Quantum Sci. Technol. 9, 045039 (2024) The optimization of quantum circuits can be hampered by a decay of average gradient amplitudes with increasing system size. When the decay is exponential, this is called the barren plateau problem. Considering explicit circuit parametrizations (in terms of rotation angles), it has been shown in Arrasmith et al (2022 Quantum Sci. Technol. 7 045015) that barren plateaus are equivalent to an exponential decay of the variance of cost-function differences. We show that the issue is particularly simple in the (parametrization-free) Riemannian formulation of such optimization problems and obtain a tighter bound for the cost-function variance. An elementary derivation shows that the single-gate variance of the cost function is strictly equal to half the variance of the Riemannian single-gate gradient, where we sample variable gates according to the uniform Haar measure. The total variances of the cost function and its gradient are then both bounded from above by the sum of single-gate variances and, conversely, bound single-gate variances from above. So, decays of gradients and cost-function variations go hand in hand, and barren plateau problems cannot be resolved by avoiding gradient-based in favor of gradient-free optimization methods.

Entanglement entropies quantify the complexity of quantum matter and its utility for quantum information processing. The computation costs in tensor network simulations are intimately related to entanglement properties of the studied systems. Also, one can use the entanglement structure to optimize the network geometry. In addition to such practical issues, the information-theoretic lens on quantum many-body systems leads to a deeper understanding of the physics. We are trying to gain a comprehensive understanding of the entanglement scaling in many-body systems:

• In [1], we derived a conformal field theory (CFT) formula for subsystem Rényi entropies and the boundary entropy in ground states of critical 1D systems.
• In [2], we established the log-area law for critical fermionic systems with a finite Fermi surface, where the proportionality constant is predicted by the Widom conjecture. We showed that there is now logarithmic correction for 2D critical fermions with a zero-dimensional Fermi surface, and confirmed the area law for critical 2D bosons.
• For collective spin models in the thermodynamic limit, we derived the logarithmic scaling of the groundstate entanglement entropy at criticality and developed a perturbative method for 1/N expansions of entanglement entropies [3,4].
• In [5], we elucidated entanglement properties of excited states in Bethe-ansatz integrable systems.
• Employing the eigenstate thermalization hypothesis (ETH), we found that the entanglement entropy of (almost) all energy eigenstates is described by a single crossover function [6,7,8]. These functions capture the full crossover from the groundstate entanglement regime at low energies and small subsystem size (area or log-area law) to the extensive volume-law regime at high energies or large subsystem size. For quantum-critical regimes, the crossover functions have universal scaling properties.
• In the context of tensor network states (TNS), we proved, for example, that MERA states in D≥22 dimensions obey the entanglement area law and are, actually a subclass of PEPS [9]. In [10], I determined the scaling of entanglement properties in thermofield double states for typical 1D systems at finite temperatures. Thermofield double states are employed in TNS simulations for at nonzero temperatures. The results of field-theoretical arguments and numerical examples, show that the cost for the classical simulation at finite temperatures grows only polynomially with the inverse temperature and is system-size independent -- even for quantum critical systems.

		[1] "Entanglement and boundary critical phenomena" H.-Q. Zhou, T. Barthel, J. O. Fjærestad, and U. Schollwöck arXiv:cond-mat/0511732, pdf, Phys. Rev. A 74, 050305(R) (2006), also in Virtual Journal of Nanoscale Science and Technology 14, Issue 23 (2006) We investigate boundary critical phenomena from a quantum information perspective. Bipartite entanglement in the ground state of one-dimensional quantum systems is quantified using the Renyi entropy S_alpha, which includes the von Neumann entropy (alpha=1) and the single-copy entanglement (alpha=infinity) as special cases. We identify the contribution from the boundary entropy to the Renyi entropy, and show that there is an entanglement loss along boundary renormalization group (RG) flows. This property, which is intimately related to the Affleck-Ludwig g-theorem, can be regarded as a consequence of majorization relations between the spectra of the reduced density matrix along the boundary RG flows. We also point out that the bulk contribution to the single-copy entanglement is half of that to the von Neumann entropy, whereas the boundary contribution is the same.
		[2] "Entanglement scaling in critical two-dimensional fermionic and bosonic systems" T. Barthel, M.-C. Chung, and U. Schollwöck arXiv:cond-mat/0602077, pdf, Phys. Rev. A 74, 022329 (2006) We relate the reduced density matrices of quadratic fermionic and bosonic models to their Green's function matrices in a unified way and calculate the scaling of the entanglement entropy of finite systems in an infinite universe exactly. For critical fermionic two-dimensional (2D) systems at $T=0$, two regimes of scaling are identified: generically, we find a logarithmic correction to the area law with a prefactor dependence on the chemical potential that confirms earlier predictions based on the Widom conjecture. If, however, the Fermi surface of the critical system is zero-dimensional, then we find an area law with a sublogarithmic correction. For a critical bosonic 2D array of coupled oscillators at $T=0$, our results show that the entanglement entropy follows the area law without corrections.
		[3] "Entanglement entropy beyond the free case" T. Barthel, S. Dusuel, and J. Vidal arXiv:cond-mat/0606436, pdf, Phys. Rev. Lett. 97, 220402 (2006), also in Virtual Journal of Nanoscale Science and Technology 14, Issue 24 (2006) We present a perturbative method to compute the ground state entanglement entropy for interacting systems. We apply it to a collective model of mutually interacting spins in a magnetic field. At the quantum critical point, the entanglement entropy scales logarithmically with the subsystem size, the system size, and the anisotropy parameter. We determine the corresponding scaling prefactors and evaluate the leading finite-size correction to the entropy. Our analytical predictions are in perfect agreement with numerical results.
		[4] "Entanglement entropy in collective models" J. Vidal, S. Dusuel, and T. Barthel arXiv:cond-mat/0610833, pdf, J. Stat. Mech. P01015 (2006) We discuss the behaviour of the entanglement entropy of the ground state in various collective systems. Results for general quadratic two-mode boson models are given, yielding the relation between quantum phase transitions of the system (signalled by a divergence of the entanglement entropy) and the excitation energies. Such systems naturally arise when expanding collective spin Hamiltonians at leading order via the Holstein-Primakoff mapping. In a second step, we analyse several such models (the Dicke model, the two-level Bardeen-Cooper-Schrieffer model, the Lieb-Mattis model and the Lipkin-Meshkov-Glick model) and investigate the properties of the entanglement entropy over the whole parameter range. We show that when the system contains gapless excitations the entanglement entropy of the ground state diverges with increasing system size. We derive and classify the scaling behaviours that can be met.
		[5] "Bound states and entanglement in the excited states of quantum spin chains" J. Mölter, T. Barthel, U. Schollwöck, and V. Alba arXiv:1407.0066, pdf, J. Stat. Mech. P10029 (2014), Special Issue "Quantum Entanglement in Condensed Matter Physics" J. Stat. Mech. 2014 We investigate the entanglement properties of the excited states of the spin-1/2 Heisenberg (XXX) chain with isotropic antiferromagnetic interactions, by exploiting the Bethe ansatz solution of the model. We consider eigenstates obtained from both real and complex solutions ('strings') of the Bethe equations. Physically, the former are states of interacting magnons, whereas the latter contain bound states of groups of particles. We first focus on the situation with few particles in the chain. Using exact results and semiclassical arguments, we derive an upper bound Smax for the entanglement entropy. This exhibits an intermediate behaviour between logarithmic and extensive, and it is saturated for highly-entangled states. As a function of the eigenstate energy, the entanglement entropy is organized in bands. Their number depends on the number of blocks of contiguous Bethe-Takahashi quantum numbers. In the presence of bound states a significant reduction in the entanglement entropy occurs, reflecting that a group of bound particles behaves effectively as a single particle. Interestingly, the associated entanglement spectrum shows edge-related levels. At a finite particle density, the semiclassical bound Smax becomes inaccurate. For highly-entangled states the subsystem entropy is proprtional to the chord length, signalling the crossover to extensive entanglement. Finally, we consider eigenstates containing a single pair of bound particles. No significant entanglement reduction occurs, in contrast with the few-particle case.
		[6] "Eigenstate entanglement: Crossover from the ground state to volume laws" Q. Miao and T. Barthel arXiv:1905.07760, pdf, Phys. Rev. Lett. 127, 040603 (2021) For the typical quantum many-body systems that obey the eigenstate thermalization hypothesis (ETH) we argue that the entanglement entropy of (almost) all energy eigenstates is described by a single crossover function. The ETH implies that the crossover functions can be deduced from subsystem entropies of thermal ensembles and have universal properties. These functions capture the full crossover from the ground-state entanglement regime at low energies and small subsystem size (area or log-area law) to the extensive volume-law regime at high energies or large subsystem size. For critical one-dimensional systems, a universal scaling function follows from conformal field theory and can be adapted for nonlinear dispersions. We use it to also deduce the crossover scaling function for Fermi liquids in d>1 dimensions. The analytical results are complemented by numerics for large noninteracting systems of fermions in1
		[7] "Scaling functions for eigenstate entanglement crossovers in harmonic lattices" T. Barthel and Q. Miao arXiv:1912.10045, pdf, Phys. Rev. A 104, 022414 (2021) For quantum matter, eigenstate entanglement entropies obey an area law or log-area law at low energies and small subsystem sizes and cross over to volume laws for high energies and large subsystems. The crossover can be described by scaling functions. We demonstrate this for the harmonic lattice model which describes quantized lattice vibrations and is a regularization for free scalar field theories, describing, e.g., spin-0 bosonic particles. In one dimension, the groundstate entanglement obeys a log-area law. For dimensions d>1, it displays area laws, even at criticality, because excitation energies vanish only at a single point in momentum space. We sample excited states. The distribution of their entanglement entropies are sharply peaked around subsystem entropies of corresponding thermodynamic ensembles in accordance with the eigenstate thermalization hypothesis. In this way, we determine entanglement scaling functions which correspond to the quantum critical regime of the model. We show how infrared singularities of the system can be handled and how to access the thermodynamic limit using a perturbative trick for the covariance matrix. Eigenstates for quasi-free bosonic systems are not Gaussian. We resolve this problem by considering appropriate squeezed states instead. For these, entanglement entropies can be evaluated efficiently.
		[8] "Eigenstate entanglement scaling for critical interacting spin chains" Q. Miao and T. Barthel arXiv:2010.07265, pdf, Quantum 6, 642 (2022) With increasing subsystem size and energy, bipartite entanglement entropies of energy eigenstates cross over from the groundstate scaling to a volume law. In previous work, we pointed out that, when strong or weak eigenstate thermalization (ETH) applies, the entanglement of all or, respectively, almost all eigenstates follow universal scaling functions which are determined by the subsystem entropy of thermal states. This was demonstrated by field-theoretical arguments and by analysis of large systems of non-interacting fermions and bosons. Here, we further substantiate such scaling properties for integrable and non-integrable interacting spin-1/2 chains at criticality using exact diagonalization. In particular, we analyze XXZ and transverse-field Ising models with and without next-nearest-neighbor interactions. We first confirm that the crossover for subsystem entropies in thermal ensembles can be described by a universal scaling function following from conformal field theory. Then, we analyze the validity of ETH for entanglement in these models. Even for the relatively small system sizes that can be simulated, the distributions of eigenstate entanglement entropies are sharply peaked around the subsystem entropies of the corresponding thermal ensembles.
		[9] "Real-space renormalization yields finitely correlated states" T. Barthel, M. Kliesch, and J. Eisert arXiv:1003.2319, pdf, Phys. Rev. Lett. 105, 010502 (2010) Real-space renormalization approaches for quantum lattice systems generate certain hierarchical classes of states that are subsumed by the multi-scale entanglement renormalization ansatz (MERA). It is shown that, with the exception of one dimension, MERA states can be efficiently mapped to finitely correlated states, also known as projected entangled pair states (PEPS), with a bond dimension independent of the system size. Hence, MERA states form an efficiently contractible class of PEPS and obey an area law for the entanglement entropy. It is shown further that there exist other efficiently contractible schemes violating the area law.
		[10] "Typical one-dimensional quantum systems at finite temperatures can be simulated efficiently on classical computers" T. Barthel arXiv:1708.09349, pdf It is by now well-known that ground states of gapped one-dimensional (1d) quantum-many body systems with short-range interactions can be studied efficiently using classical computers and matrix product state techniques. A corresponding result for finite temperatures was missing. For 1d systems that can be described by a local 1+1d field theory, it is shown here that the cost for the classical simulation at finite temperatures grows in fact only polynomially with the inverse temperature and is system-size independent -- even for quantum critical systems. In particular, we show that the thermofield double state (TDS), a purification of the equilibrium density operator, can be obtained efficiently in matrix-product form. The argument is based on the scaling behavior of Rényi entanglement entropies in the TDS. At finite temperatures, they obey the area law. For quantum critical systems, the entanglement is found to grow only logarithmically with inverse temperature, S~log(beta). The field-theoretical results are confirmed by quasi-exact numerical simulations of quantum magnets and interacting bosons.
		[11] "Fundamental limitations for measurements in quantum many-body systems" T. Barthel and J. Lu arXiv:1802.04378, pdf, Phys. Rev. Lett. 121, 080406 (2018) Dynamical measurement schemes are an important tool for the investigation of quantum many-body systems, especially in the age of quantum simulation. Here, we address the question whether generic measurements can be implemented efficiently if we have access to a certain set of experimentally realizable measurements and can extend it through time evolution. For the latter, two scenarios are considered (a) evolution according to unitary circuits and (b) evolution due to Hamiltonians that we can control in a time-dependent fashion. We find that the time needed to realize a certain measurement to a predefined accuracy scales in general exponentially with the system size -- posing a fundamental limitation. The argument is based, on the construction of epsilon-packings for manifolds of observables with identical spectra and a comparison of their cardinalities to those of epsilon-coverings for quantum circuits and unitary time-evolution operators. The former is related to the study of Grassmann manifolds.

Nonequilibrium dynamics in closed systems. - Fundamental questions concerning the quantum dynamics of many-body systems are about the equilibration and thermalization after the system is driven out of equilibrium. In an early contribution, I derived general preconditions under which subsystems of large quasi-free systems converge to steady states [1]. These steady states are derived from generalized Gibbs ensembles and obey entanglement area laws. The criteria also lead to simple counter-examples, where systems fail to equilibrate. We have also studied out-of-equilibrium dynamics for non-integrable systems of ultracold atoms in optical lattices, e.g., concerning expansion dynamics and localization [2], quench dynamics and equilibration in systems with multiple atom species [3], as well as domain-wall melting in the presence of experimentally relevant hole and spin-flip defects [4]. We are currently working on new tensor network state methods to extend reachable times in the study of non-equilibrium dynamics, starting from low-entangled states.

Decoherence, criticality, and phase transitions in driven-dissipative systems. - In several experimental frameworks, a high level of control on complex quantum systems has been accomplished. Prominent examples are circuits of superconducting qubits, coupled quantum dots, ultracold atoms in optical lattices or tweezers, which can also be excited into Rydberg states to achieve strong interactions, ions in electromagnetic traps, and polaritons in circuit-QED or semiconductor-microcavity systems. Due to practical constraints and our aim of manipulating these systems efficiently, they are inevitably open in the sense that they are coupled to the environment. Such environment couplings naturally lead to dissipation and decoherence, which pose challenges for modern quantum technology. On the other hand, driving and dissipation in open systems could be designed to stabilize useful states or to drive the system into (novel) phases of matter that may not be accessible in equilibrium. In [5], we analyzed for example the stabilization of Bose-Einstein condensates through driving and dissipation. So far, our work is focused on Markovian quantum systems. Their dynamics is governed by a Lindblad master equation which, in addition to the unitary Hamiltonian part, comprises Lindblad operators that capture the environment couplings:

• We have used Trotter decompositions and Lieb-Robinson bound techniques to prove that Markovian dynamics can be simulated efficiently on unitary quantum computers and with system-size independent costs on classical computers [6,7]. We employ tensor network methods that rest on these findings to simulate large open quantum systems.
• We showed, for example, that the internal interactions of a dissipative system can hinder decoherence. And, under certain circumstances, the interplay of the internal interactions and dissipation processes can lead to divergent decoherence times, i.e., a critical system with an algebraic coherence decay [8].
• In [9], we discussed how block-triangular structures in the Liouvillian super-operator, corresponding to (possibly hidden) dynamical constraints, can make driven-dissipative phase transitions impossible.
• Based on work by Prosen and Seligman, we laid out a formalism for the solution of quasi-free and of quadratic open fermion and boson systems [10].
• In [11], we used this formalism to elucidate fundamental differences concerning criticality in open and closed systems, and found that, without symmetry constraints beyond invariance under single-particle basis and particle-hole transformations, all gapped quadratic Liouvillians belong to the same phase. We are currently employing Keldysh field theory and renormalization group analyses to investigate driven-dissipative phase transitions in interacting systems with a focus on fundamental differences to closed systems [5].

		[1] "Dephasing and the steady state in quantum many-particle systems" T. Barthel and U. Schollwöck arXiv:0711.4896, pdf, Phys. Rev. Lett. 100, 100601 (2008), also in Virtual Journal of Quantum Information 8, Issue 3 (2008) We discuss relaxation in bosonic and fermionic many-particle systems. For integrable systems, time evolution can cause a dephasing effect, leading for finite subsystems to steady states. We explicitly derive those steady subsystem states and devise sufficient prerequisites for the dephasing to occur. We also find simple scenarios, in which dephasing is ineffective and discuss the dependence on dimensionality and criticality. It follows further that, after a quench of system parameters, entanglement entropy will become extensive. This provides a way of creating strong entanglement in a controlled fashion.
		[2] "Quasiperiodic Bose-Hubbard model and localization in one-dimensional cold atomic gases" G. Roux, T. Barthel, I. P. McCulloch, C. Kollath, U. Schollwöck, and T. Giamarchi arXiv:0802.3774, pdf, Phys. Rev. A 78, 023628 (2008) We compute the phase diagram of the one-dimensional Bose-Hubbard model with a quasiperiodic potential by means of the density-matrix renormalization group technique. This model describes the physics of cold atoms loaded in an optical lattice in the presence of a superlattice potential whose wavelength is incommensurate with the main lattice wavelength. After discussing the conditions under which the model can be realized experimentally, the study of the density vs the chemical potential curves for a nontrapped system unveils the existence of gapped phases at incommensurate densities interpreted as incommensurate charge-density-wave phases. Furthermore, a localization transition is known to occur above a critical value of the potential depth V_2 in the case of free and hard-core bosons. We extend these results to soft-core bosons for which the phase diagrams at fixed densities display new features compared with the phase diagrams known for random box distribution disorder. In particular, a direct transition from the superfluid phase to the Mott-insulating phase is found at finite V_2. Evidence for reentrances of the superfluid phase upon increasing interactions is presented. We finally comment on different ways to probe the emergent quantum phases and most importantly, the existence of a critical value for the localization transition. The latter feature can be investigated by looking at the expansion of the cloud after releasing the trap.
		[3] "Magnetism, coherent many-particle dynamics, and relaxation with ultracold bosons in optical superlattices" T. Barthel, C. Kasztelan, I. P. McCulloch, and U. Schollwöck arXiv:0809.5141, pdf, Phys. Rev. A 79, 053627 (2009) We study how well magnetic models can be implemented with ultracold bosonic atoms of two different hyperfine states in an optical superlattice. The system is captured by a two-species Bose-Hubbard model, but realizes in a certain parameter regime actually the physics of a spin-1/2 Heisenberg magnet, describing the second-order hopping processes. Tuning of the superlattice allows for controlling the effect of fast first-order processes versus the slower second-order ones. Using the density-matrix renormalization-group method, we provide the evolution of typical experimentally available observables. The validity of the description via the Heisenberg model, depending on the parameters of the Hubbard model, is studied numerically and analytically. The analysis is also motivated by recent experiments [S. Foelling et al., Nature (London) 448, 1029 (2007); S. Trotzky et al., Science 319, 295 (2008)] where coherent two-particle dynamics with ultracold bosonic atoms in isolated double wells were realized. We provide theoretical background for the next step, the observation of coherent many-particle dynamics after coupling the double wells. Contrary to the case of isolated double wells, relaxation of local observables can be observed. The tunability between the Bose-Hubbard model and the Heisenberg model in this setup could be used to study experimentally the differences in equilibration processes for nonintegrable and Bethe ansatz integrable models. We show that the relaxation in the Heisenberg model is connected to a phase averaging effect, which is in contrast to the typical scattering driven thermalization in nonintegrable models. We discuss the preparation of magnetic ground states by adiabatic tuning of the superlattice parameters.
		[4] "Domain-wall melting in ultracold-boson systems with hole and spin-flip defects" J. C. Halimeh, A. Wöllert, I. P. McCulloch, U. Schollwöck, and T. Barthel arXiv:1307.0513, pdf, Phys. Rev. A 89, 063603 (2014) Quantum magnetism is a fundamental phenomenon of nature. As of late, it has garnered a lot of interest because experiments with ultracold atomic gases in optical lattices could be used as a simulator for phenomena of magnetic systems. A paradigmatic example is the time evolution of a domain-wall state of a spin-1/2 Heisenberg chain, the so-called domain-wall melting. The model can be implemented by having two species of bosonic atoms with unity filling and strong on-site repulsion U in an optical lattice. In this paper, we study the domain-wall melting in such a setup on the basis of the time-dependent density matrix renormalization group (tDMRG). We are particularly interested in the effects of defects that originate from an imperfect preparation of the initial state. Typical defects are holes (empty sites) and flipped spins. We show that the dominating effects of holes on observables like the spatially resolved magnetization can be taken account of by a linear combination of spatially shifted observables from the clean case. For sufficiently large U, further effects due to holes become negligible. In contrast, the effects of spin flips are more severe as their dynamics occur on the same time scale as that of the domain-wall melting itself. It is hence advisable to avoid preparation schemes that are based on spin-flips.
		[5] "Driven-dissipative Bose-Einstein condensation and the upper critical dimension" Y. Zhang and T. Barthel arXiv:2311.13561, pdf, Phys. Rev. A 109, L021301 (2024) Driving and dissipation can stabilize Bose-Einstein condensates. Using Keldysh field theory, we analyze this phenomenon for Markovian systems that can comprise on-site two-particle driving, on-site single-particle and two-particle loss, as well as edge-correlated pumping. Above the upper critical dimension, mean-field theory shows that pumping and two-particle driving induce condensation right at the boundary between the stable and unstable regions of the non-interacting theory. With nonzero two-particle driving, the condensate is gapped. This picture is consistent with the recent observation that, without symmetry constraints beyond invariance under single-particle basis transformations, all gapped quadratic bosonic Liouvillians belong to the same phase. For systems below the upper critical dimension, the edge-correlated pumping penalizes high-momentum fluctuations, rendering the theory renormalizable. We perform the one-loop renormalization group analysis, finding a condensation transition inside the unstable region of the non-interacting theory. Interestingly, its critical behavior is determined by a Wilson-Fisher-like fixed point with universal correlation-length exponent nu=0.6 in three dimensions.
		[6] "Dissipative quantum Church-Turing theorem" M. Kliesch, T. Barthel, C. Gogolin, M. Kastoryano, and J. Eisert arXiv:1105.3986, pdf, Phys. Rev. Lett. 107, 120501 (2011), see also D. Browne "Viewpoint: Quantum simulation hits the open road" Physics 4, 72 (2011) We show that the time evolution of an open quantum system, described by a possibly time dependent Liouvillian, can be simulated by a unitary quantum circuit of a size scaling polynomially in the simulation time and the size of the system. An immediate consequence is that dissipative quantum computing is no more powerful than the unitary circuit model. Our result can be seen as a dissipative Church-Turing theorem, since it implies that under natural assumptions, such as weak coupling to an environment, the dynamics of an open quantum system can be simulated efficiently on a quantum computer. Formally, we introduce a Trotter decomposition for Liouvillian dynamics and give explicit error bounds. This constitutes a practical tool for numerical simulations, e.g., using matrix-product operators. We also demonstrate that most quantum states cannot be prepared efficiently.
		[7] "Quasi-locality and efficient simulation of Markovian quantum dynamics" T. Barthel and M. Kliesch arXiv:1111.4210, pdf, Phys. Rev. Lett. 108, 230504 (2012) We consider open many-body systems governed by a time-dependent quantum master equation with short-range interactions. With a generalized Lieb-Robinson bound, we show that the evolution in this very generic framework is quasilocal; i.e., the evolution of observables can be approximated by implementing the dynamics only in a vicinity of the observables' support. The precision increases exponentially with the diameter of the considered subsystem. Hence, time evolution can be simulated on classical computers with a cost that is independent of the system size. Providing error bounds for Trotter decompositions, we conclude that the simulation on a quantum computer is additionally efficient in time. For experiments and simulations in the Schrödinger picture, our result can be used to rigorously bound finite-size effects.
		[8] "Algebraic versus exponential decoherence in dissipative many-particle systems" Z. Cai and T. Barthel arXiv:1304.6890, pdf, Phys. Rev. Lett. 111, 150403 (2013), covered by Phys.org "Quantum particles find safety in numbers" phys.org, 10/16/2013 The interplay between dissipation and internal interactions in quantum many-body systems gives rise to a wealth of novel phenomena. Here we investigate spin-1/2 chains with uniform local couplings to a Markovian environment using the time-dependent density matrix renormalization group (tDMRG). For the open XXZ model, we discover that the decoherence time diverges in the thermodynamic limit. The coherence decay is then algebraic instead of exponential. This is due to a vanishing gap in the spectrum of the corresponding Liouville superoperator and can be explained on the basis of a perturbative treatment. In contrast, decoherence in the open transverse-field Ising model is found to be always exponential. In this case, the internal interactions can both facilitate and impede the environment-induced decoherence.
		[9] "Superoperator structures and no-go theorems for dissipative quantum phase transitions" T. Barthel and Y. Zhang arXiv:2012.05505, pdf, Phys. Rev. A 105, 052224 (2022) In the thermodynamic limit, the steady states of open quantum many-body systems can undergo nonequilibrium phase transitions due to a competition between coherent and driven-dissipative dynamics. Here, we consider Markovian systems and elucidate structures of the Liouville superoperator that generates the time evolution. In many cases of interest, an operator-basis transformation can bring the Liouvillian into a block-triangular form, making it possible to assess its spectrum. The spectral gap sets the asymptotic decay rate. The superoperator structure can be used to bound gaps from below, showing that, in a large class of systems, dissipative phase transitions are actually impossible and that the convergence to steady states follows an exponential temporal decay. Furthermore, when the blocks on the diagonal are Hermitian, the Liouvillian spectra obey Weyl ordering relations. The results apply, for example, to Davies generators and quadratic systems and are also demonstrated for various spin models.
		[10] "Solving quasi-free and quadratic Lindblad master equations for open fermionic and bosonic systems" T. Barthel and Y. Zhang arXiv:2112.08344, pdf, J. Stat. Mech. 113101 (2022) The dynamics of Markovian open quantum systems are described by Lindblad master equations. For fermionic and bosonic systems that are quasi-free, i.e., with Hamiltonians that are quadratic in the ladder operators and Lindblad operators that are linear in the ladder operators, we derive the equation of motion for the covariance matrix. This determines the evolution of Gaussian initial states and the steady states, which are also Gaussian. Using ladder super-operators (a.k.a. third quantization), we show how the Liouvillian can be transformed to a many-body Jordan normal form which also reveals the full many-body spectrum. Extending previous work by Prosen and Seligman, we treat fermionic and bosonic systems on equal footing with Majorana operators, shorten and complete some derivations, also address the odd-parity sector for fermions, give a criterion for the existence of bosonic steady states, cover non-diagonalizable Liouvillians also for bosons, and include quadratic systems. In extension of the quasi-free open systems, quadratic open systems comprise additional Hermitian Lindblad operators that are quadratic in the ladder operators. While Gaussian states may then evolve into non-Gaussian states, the Liouvillian can still be transformed to a useful block-triangular form, and the equations of motion for k-point Green's functions form a closed hierarchy. Based on this formalism, results on criticality and dissipative phase transitions in such models are discussed in a companion paper.
		[11] "Criticality and phase classification for quadratic open quantum many-body systems" Y. Zhang and T. Barthel arXiv:2204.05346, pdf, Phys. Rev. Lett. 129, 120401 (2022) We study the steady states of translation-invariant open quantum many-body systems governed by Lindblad master equations, where the Hamiltonian is quadratic in the ladder operators, and the Lindblad operators are either linear or quadratic and Hermitian. These systems are called quasifree and quadratic, respectively. We find that steady states of one-dimensional systems with finite-range interactions necessarily have exponentially decaying Green's functions. For the quasifree case without quadratic Lindblad operators, we show that fermionic systems with finite-range interactions are noncritical for any number of spatial dimensions and provide bounds on the correlation lengths. Quasifree bosonic systems can be critical in D>1 dimensions. Last, we address the question of phase transitions in quadratic systems and find that, without symmetry constraints beyond invariance under single-particle basis and particle-hole transformations, all gapped Liouvillians belong to the same phase.

Somewhat surprisingly, we found that tensor network approximations and ideas from quantum Monte Carlo can be adapted to study stochastic nonequilibrium dynamics in classical networks. This is an exciting branch of nonequilibrium statistical physics with a wide range of applications such as spin glasses, technical networks, epidemic dynamics, and social networks. A common technique for the investigation of stochastic dynamics is Markov-chain Monte Carlo (MCMC), where one starts from an initial state and propagates in time according to a stochastic model. However, in this approach, observables converge very slowly when increasing the number of samples. Hence, MCMC does not allow for an efficient study of temporal correlations and rare (but consequential) events like big epidemic outbreaks or large-scale outages in technical networks.

• While tensor networks are used to encode entangled states in quantum many-body systems, we found that they can be employed to efficiently encode conditional probabilities for state trajectories in stochastic network dynamics [1,2], reducing the time complexity from exponential to polynomial. With simulations of Glauber-Ising dynamics on Boolean networks, we have established that the new tensor network technique has substantial advantages over MCMC methods. It gives access to precise temporal correlations, which are otherwise very difficult to assess in systems with bidirectional dynamics. There are still some limitations: The new tensor network technique is based on the dynamic cavity method which is an extension of Mézard's cavity method for equilibrium systems and, hence, restricted to locally tree-like networks. Also, computation costs increase strongly with increasing vertex degrees.
• Recently, I have developed an alternative Monte Carlo algorithm that is not restricted to tree-like graphs and allows us to sample stochastic dynamics conditioned on specific (rare) events. New efficient simulation techniques are needed to properly understand the preconditions of rare catastrophic events and to improve mitigation strategies. For example, large epidemic outbreaks often emerge via a confluence of rare, stochastic events. Potentially dangerous spillover infections from animals into humans are frequent but rarely take off: The introduced infectious disease typically dies quickly or spreads to a small number of local cases and fizzles out. Only a small percentage of these local outbreaks manage to pass bottlenecks in the heterogeneous contact network, to expand, and become widespread epidemics.

		[1] "A matrix product algorithm for stochastic dynamics on networks, applied to non-equilibrium Glauber dynamics" T. Barthel, C. De Bacco, and S. Franz arXiv:1508.03295, pdf, Phys. Rev. E 97, 010104(R) (2018) We introduce and apply a novel efficient method for the precise simulation of stochastic dynamical processes on locally tree-like graphs. Networks with cycles are treated in the framework of the cavity method. Such models correspond, for example, to spin-glass systems, Boolean networks, neural networks, or other technological, biological, and social networks. Building upon ideas from quantum many-body theory, the new approach is based on a matrix product approximation of the so-called edge messages -- conditional probabilities of vertex variable trajectories. Computation costs and accuracy can be tuned by controlling the matrix dimensions of the matrix product edge messages (MPEM) in truncations. In contrast to Monte Carlo simulations, the algorithm has a better error scaling and works for both, single instances as well as the thermodynamic limit. We employ it to examine prototypical non-equilibrium Glauber dynamics in the kinetic Ising model. Because of the absence of cancellation effects, observables with small expectation values can be evaluated accurately, allowing for the study of decay processes and temporal correlations.
		[2] "The matrix product approximation for the dynamic cavity method" T. Barthel arXiv:1904.03312, pdf, J. Stat. Mech. 013217 (2020) Stochastic dynamics of classical degrees of freedom, defined on vertices of locally tree-like graphs, can be studied in the framework of the dynamic cavity method which is exact for tree graphs. Such models correspond for example to spin-glass systems, Boolean networks, neural networks, and other technical, biological, and social networks. The central objects in the cavity method are edge messages - conditional probabilities of two vertex variable trajectories. In this paper, we discuss a rather pedagogical derivation for the dynamic cavity method, give a detailed account of the novel matrix product edge message (MPEM) algorithm for the solution of the dynamic cavity equation as introduced in Phys. Rev. E 97, 010104(R) (2018), and present optimizations and extensions. Matrix product approximations of the edge messages are constructed recursively in an iteration over time. Computation costs and precision can be tuned by controlling the matrix dimensions of the MPEM in truncations. Without truncations, the dynamics is exact. Data for Glauber-Ising dynamics shows a linear growth of computation costs in time. In contrast to Monte Carlo simulations, the approach has a much better error scaling. Hence, it gives for example access to low probability events and decaying observables like temporal correlations. We discuss optimized truncation schemes and an extension that allows to capture models which have a continuum time limit.