Simulation methods for open quantum many-body systems

TL;DR

This paper reviews various computational methods for simulating open quantum many-body systems, emphasizing recent progress in handling large systems and their steady states.

Contribution

It provides a comprehensive overview of recent simulation techniques tailored for open quantum systems, highlighting advancements over traditional methods for closed systems.

Findings

Recent methods enable simulation of larger systems.

Progress in understanding steady states of open systems.

Advances in relaxation dynamics modeling.

Abstract

Coupling a quantum many-body system to an external environment dramatically changes its dynamics and offers novel possibilities not found in closed systems. Of special interest are the properties of the steady state of such open quantum many-body systems, as well as the relaxation dynamics towards the steady state. However, new computational tools are required to simulate open quantum many-body systems, as methods developed for closed systems cannot be readily applied. We review several approaches to simulate open many-body systems and point out the advances made in recent years towards the simulation of large system sizes.

Tables

Table 1. Table 1: Ground state calculation in a closed quantum system (left) and Steady state calculation in an open quantum system (right). The former one requires an imaginary time evolution while the latter follows a real time evolution. Both the Hamiltonian H 𝐻 H and the vectorized Liouvillian ℒ ♯ subscript ℒ ♯ \mathcal{L}_{\sharp} can be decomposed as a sum of local terms. | e 0 ⟩ ket subscript 𝑒 0 |e_{0}\rangle is the ground state of the many-body Hamiltonian with e 0 subscript 𝑒 0 e_{0} as its ground state. | ρ s ⟩ ♯ subscript ket subscript 𝜌 𝑠 ♯ |\rho_{s}\rangle_{\sharp} is the non-equilibrium steady state of the Liouvillian in their vectorized forms.

Ground states	Steady states
$H={\sum }_{⟨i,j⟩}{h}^{[i,j]}$	${ℒ}_{♯}={\sum }_{⟨i,j⟩}{ℒ}_{♯}^{[i,j]}$
$e^{-Ht}$	$e^{{ℒ}_{♯}t}$
$|e_0⟩$	${|{\rho }_s⟩}_{♯}$
$⟨e_0\left|H\right|e_0⟩=e_0$	${}_{♯}{}^{}⟨|{\rho }_s⟩{ℒ}_{♯}|{\rho }_s⟩{}_{♯}=0$
Imaginary time	Real time

Table 2. Table 2: Table comparing the different simulation methods discussed in this review. We differentiate the methods by the system sizes that can be simulated, the spatial dimensions, contraints on the local Hilbert space dimension, whether fermionic systems can be treated, the simulation performance for inhomogeneous systems, and whether the correct critical exponents of phase transitions can be obtained.

	WFMC	TN	Variational Principle	VQMC	CMF	TWA
System size	20	TDL	TDL	16	TDL	400
Dimensions	1D,2D	1D,2D	any111Works better in higher dimensions	any	any222Works better in higher dimensions	any
Local Hilbert space	small	small	large	large	small	large
Fermionic systems	Yes	Yes	partially	No	partially	unknown
Inhomogeneous systems	good	good	bad	good	good	good
Critical exponents	good	good	good333For states with thermal statistics	unknown	bad	unknown

Equations

$$\rho =\sum _i{p}_i|{\psi }_i⟩⟨{\psi }_i|,$$

$$\begin{array}{rl}\frac{d}{dt}\rho & =\mathcal{L}[\rho ]\\ & =-i\left[H,\rho \right]+\sum _{\mu }\left(L_{\mu }\rho L_{\mu }^{†}-\frac{1}{2}{L}_{\mu }^{†}{L}_{\mu }\rho -\frac{1}{2}\rho L_{\mu }^{†}{L}_{\mu }\right),\end{array}$$

$$H=\frac{h}{2}\sum _i{\sigma }_x^{(i)}+\frac{V}{4}\sum _{⟨ij⟩}{\sigma }_z^{(i)}{\sigma }_z^{(j)},$$

$$H=-J\sum _{⟨i,j⟩}{b}_i^{†}{b}_j+\sum _i\left[\frac{U}{2}{n}_i^2-\Delta \omega n_i+F\left(b_i+b_i^{†}\right)\right].$$

$$\rho =\prod _i({\rho }_i^{MF}+\delta {\rho }_i),$$

$$⟨O⟩=\text{Tr}\left\{O\rho \right\}=\sum _i{p}_i⟨{\psi }_i|O|{\psi }_i⟩.$$

$$d{\psi }_i(t)=-i{H}_{eff}|{\psi }_i(t)⟩dt+\sum _j{M}_j|{\psi }_i(t)⟩d{W}_j,$$

$$H_{eff}=H+\sum _{j}2⟨c_j^{†}⟩c_j-c_j^{†}{c}_j-⟨c_j^{†}⟩⟨c_j⟩.$$

$$M_j=c_j-⟨c_j⟩.$$

$$H_{NH}=H-\frac{i}{2}\sum _j{c}_j^{†}{c}_j.$$

$$p_j=\mathcal{N}⟨{\psi }_i|c_j^{†}{c}_j|{\psi }_i⟩,$$

$$\chi =N^{\alpha }\tilde{\chi }(\lambda ),$$

$$|\psi ⟩=\sum _{s_1,\dots ,s_N=1}^d{A}_1^{s_1}\dots A_N^{s_N}|s_1,\dots ,s_N⟩$$

$$\begin{gathered} \rho =\sum _{s_1,s_1^{\prime },\dots ,s_N,s_N^{\prime }=1}^d(M_1^{s_1,s_1^{\prime }}\dots M_N^{s_N,s_N^{\prime }})|s_1,\dots ,s_N⟩ \\ \times ⟨s_1^{\prime },\dots ,s_N^{\prime }|, \end{gathered}$$

$$M_k^{s,s^{\prime }}=\sum _{a=1}^{d_k}{A}_k^{s,a}\otimes (A_k^{s^{\prime },a}{)}^{\ast }.$$

$$|\Psi ⟩=\sum _{s_1,\dots ,s_N}\sum _{a_1,\dots ,a_N}{A}_1^{s_1,a_1}\dots A_N^{s_N,a_N}|s_1{a}_1,\dots ,s_N{a}_N⟩$$

$$|\rho {⟩}_{♯}=\sum _{i_1=0}^{d^2-1}\cdots \sum _{i_N=0}^{d^2-1}{c}_{i_1\cdots i_N}|i_1{⟩}_{♯}\otimes \cdots \otimes |i_N{⟩}_{♯}.$$

$$[X{]}_{r_1,\dots ,r_N}^{s_1,\dots ,s_N}=\sum _{m_1,\dots ,m_{N-1}}{A}_{m_1}^{[1]s_1,r_1}{A}_{m_1,m_2}^{[2]s_2,r_2}\dots A_{m_{N-1}}^{[N]s_N,r_N}$$

$$\mathcal{H}=\sum _{r\in \mathbb{Z}}({\hat{\mathcal{L}}}_r^{†}\widehat{{\mathcal{L}}_r}{)}^{1/k}$$

$$|{\rho }_G⟩\approx \underset{\tau \to \infty }{\lim }\frac{e^{-\mathcal{H}\tau }|{\rho }_0⟩}{||e^{-\mathcal{H}\tau }|{\rho }_0⟩||}$$

$$|{\rho }_S⟩\approx \underset{T\to \infty }{\lim }\frac{e^{\mathcal{L}T}|{\rho }_G⟩}{||e^{\mathcal{L}T}|{\rho }_G⟩||}$$

$$\begin{array}{r}{\rho }^A=\sum _i{p}_i^A|{\varphi }_i^A⟩⟨{\varphi }_i^A|,\\ {\rho }^B=\sum _i{p}_i^B|{\varphi }_i^B⟩⟨{\varphi }_i^B|,\end{array}$$

$$\frac{d}{dt}|\rho {⟩}_{♯}={\mathcal{L}}_{♯}|\rho {⟩}_{♯}$$

$$\begin{array}{rl}{\mathcal{L}}_{♯}& \equiv -i\left(H\otimes \mathbb{I}-\mathbb{I}\otimes H^T\right)\\ & +\sum _{\mu }\left(L_{\mu }\otimes L_{\mu }^{\ast }-\frac{1}{2}{L}_{\mu }^{†}{L}_{\mu }\otimes \mathbb{I}-\frac{1}{2}\mathbb{I}\otimes L_{\mu }^{\ast }{L}_{\mu }^T\right).\end{array}$$

$$|\rho (t){⟩}_{♯}=e^{{\mathcal{L}}_{♯}t}|\rho (0){⟩}_{♯},$$

$${\mathcal{L}}_{♯}|{\rho }_s{⟩}_{♯}=0.$$

$$H=\sum _{⟨i,j⟩}\left(J_x{\sigma }_x^i{\sigma }_x^j+J_y{\sigma }_y^i{\sigma }_y^j+J_z{\sigma }_z^i{\sigma }_z^j\right)$$

$$D=\sum _{ij\in \mathcal{T}}\text{Tr}\left\{|{\dot{\rho }}_{ij}|\right\},$$

$$D[\varphi ]=\int dx\sum _m{v}_m[\nabla \varphi (x){]}^m+\sum _n{u}_n[\varphi (x){]}^n.$$

$$Z_{\text{eff}}=\int \mathcal{D}\varphi \exp \left(-{\beta }_{\text{eff}}D[\varphi ]\right).$$

Full text

1. I Introduction

2. I.1 The open quantum many-body problem

3. I.2 The Markovian quantum master equation

4. I.3 Steady state solution versus time evolution

5. I.4 Differences to equilibrium problems

6. I.5 Paradigmatic models

7. II Stochastic methods

8. III Tensor network methods

9. III.1 One spatial dimension

10. III.2 Extensions to higher dimensions

11. IV Variational methods

12. IV.1 The variational principle for open quantum systems

13. IV.2 Comparison with mean-field methods

14. IV.3 Variational tensor network methods

15. IV.4 Variational quantum Monte-Carlo methods

16. V Phase space and related methods

17. V.1 Truncated Wigner approximation

18. V.2 BBGKY hierarchy equations

19. VI Linked cluster expansion methods

20. VII Summary and outlook

I Introduction

I.1 The open quantum many-body problem

Open quantum many-body systems have witnessed a surge of interest in recent years, chiefly for two reasons. On the one hand, these systems offer the exciting possibility to use controlled dissipation channels to engineer interesting quantum many-body states as the stationary state of their dynamics Diehl et al. (2008); Verstraete et al. (2009); Weimer et al. (2010). On the other hand, open quantum many-body systems are attractive from a fundamental perspective, as their dynamics exhibits a wide range of features not found in equilibrium systems. As in the case of closed quantum systems, the complexity of the problem scales exponentially with the size of the system, requiring the use of sophisticated simulation methods to obtain useful results.

Interestingly, open quantum many-body systems are even harder to simulate on classical computers than closed systems, while at the same time the stationary state of an open quantum system is much easier to experimentally prepare than the ground state of a closed system. These properties make open quantum systems one of the prime candidates to show a quantum advantage of quantum simulators over classical methods within noisy intermediate-scale quantum devices Preskill (2018). However, this requires a thorough assessment of the capabilities of classical simulation methods, which we will provide in this review.

In our review, we first provide a general introduction to open quantum many-body systems, laying particular emphasis on the key differences compared to simulating closed quantum systems and on the paradigmatic models that have emerged to benchmark simulation methods for open systems. In the main part, we first review stochastic methods commonly known as wave-function Monte-Carlo techniques, which are based on a numerical exact treatment of the total Hilbert space of the problem. We then turn to tensor network simulation techniques aiming to describe the “physical corner” of the Hilbert space, i.e., the quantum states that are most relevant to describe the dynamical evolution and steady states of open quantum many-body systems. Subsequently, we review variational methods that employ very similar strategies, including variational methods that are based on a tensor network description. We also cover phase space methods and closely related counterparts. Finally, we have added a section on linked cluster expansion. Within our review, we will not cover methods derived from a field-theoretical description of open quantum systems within the Keldysh formalism, as this has already been extensively covered in a previous review article Sieberer et al. (2016). We will also not cover integrable models Prosen (2008); Medvedyeva et al. (2016); Foss-Feig et al. (2017); Prosen (2011b, a); Guo and Poletti (2018) for which analytical techniques such as the Bethe ansatz can be employed.

I.2 The Markovian quantum master equation

The state of an open system is described by its density operator $\rho$, which can be described as a statistical ensembles of pure states,

[TABLE]

where $p_{i}$ denotes the probability to find the system in the state $\mathopen{|}\psi_{i}\mathclose{\rangle}$. Note that the decomposition into pure states is not unique. In our review, we will limit ourselves to the discussion of Markovian systems, i.e., dynamical systems in which the generator of the dynamics $\mathcal{L}[\rho]$ (commonly called the Liouvillian) depends only on the state at the present time $t$ and not on the state at earlier times. Such Markovian systems form a dynamical semigroup and can be described by a quantum master equation in Lindblad form

[TABLE]

where $H$ is the Hamiltonian of the system and $\{L_{\mu},L_{\mu}^{\dagger}\}$ the Lindblad operators responsible for the incoherent dynamics arising from the coupling to an external environment, which are also known as the jump operators Gorini et al. (1976); Lindblad (1976).

The validity of the Lindblad master equation Eq. (2) for a concrete physical system crucially depends on the separation of several timescales. Considering a system of interest coupled to a larger environment, one first assumes a weak coupling between system and environment, such that the entanglement between system and environment remains low. Furthermore, the environment must not retain any memory of the system degrees of freedom. The approximations related to these conditions are commonly refered to as the Born-Markov approximation Breuer and Petruccione (2002) and require that the correlation time of the environment $\tau_{E}$ is much smaller than the relaxation time of the system $\tau_{R}$. Finally, the differences in eigenfrequencies in the system $\omega_{s}$ has to be large compared to the inverse relaxation time $\tau_{R}^{-1}$.

These approximations are well justified in quantum optical systems, in particular atoms coupled to electronically excited states Raitzsch et al. (2009); Baumann et al. (2010); Barreiro et al. (2011); Krauter et al. (2011); Malossi et al. (2014). There, the optical frequencies of the transition leave a large time scale to observe the complete relaxation to its equilibrium state. Additionally, the relaxation of the electronic excitation into the vacuum of the radiation field as the correlation time of the radiation field is related to the photon frequency Breuer and Petruccione (2002), which again is much larger than the relaxation rate $\tau_{R}^{-1}$. Artificial atomic systems such as the nitrogen-vacancy center in diamond Jelezko et al. (2004); Dutt et al. (2007); Robledo et al. (2011) offer similar benefits.

Another advantage of quantum optical systems for studying open quantum many-body systems is the possibility to drive them with time-dependent laser fields. Importantly, if all the jump operators in the master equation describe transitions between the eigenstates of the Hamiltonian, the resulting steady state of the system is guaranteed to be a thermal state Breuer and Petruccione (2002). However, if an oscillatory driving term is added to the system Hamiltonian, it is possible to observe non-equilibrium steady states in the rotating frame of the driving. Optically excited atoms can also exhibit strong interactions when excited to Rydberg states Saffman et al. (2010), which can be used to realize a rich variety of driven-dissipative quantum many-body systems Lee et al. (2011); Glaetzle et al. (2012); Ates et al. (2012); Lemeshko and Weimer (2013); Rao and Mølmer (2013); Carr and Saffman (2013).

There are also interesting solid state platforms to study strong interaction and dissipation. One example is that of semiconductor polaritonic systems see, e.g., Carusotto and Ciuti (2013)), where semiconductor microstructures are used to embed quantum wells or quantum dots, becoming a photonic resonator where strong interactions can be induced. Another example is that of circuit-QED systems Fitzpatrick et al. (2017); Ma et al. (2019), where superconducting circuits can be used to construct Bose-Hubbard lattices of microwave photons, and where dissipation can be engineered, so that one can have a tailored reservoir.

It is important to remark that the Lindblad operators are usually considered to be local, but this approximation holds only in the weak-coupling limit. To be more precise, a Markovian master equation with quasi-local Lindblad operators holds as long as the coupling between the system and the environment is weak, which in practice amounts to (1) a slow development of correlations between system and environment, (2) fast decay of excitations of the environment, and (3) neglect of fast-oscillating terms when compared to the typical system timescale. One should be careful, however, since when dealing with strongly correlated systems, strong interactions within the system of interest may lead to a breakdown of the local Lindblad dissipation Wichterich et al. (2007); Beaudoin et al. (2011). In these cases, it may be necessary to consider additional steps to derive the correct Lindblad operators Reiter and Sørensen (2012). For the purpose of our review, we assume that the correct Lindblad form has already been derived.

I.3 Steady state solution versus time evolution

Typically, there are two different aspects that are of interest when studying open quantum many-body systems. First, one wishes to understand the properties of one or several steady states that the system reaches in the long time limit. This is similar to understanding the ground state properties of a closed many-body system. Second, one is interested in the dynamical evolution of the system towards the steady state. The latter is particularly interesting when the system exhibits several steady states that can be reached depending on the initial condition of the system.

While the requirements for the appearance of a unique steady state are well understood for finite systems Spohn (1976), many-body systems add the additional complication that the long time limit and the thermodynamic limit do not necessarily commute. In some cases, even when chiefly interested in the steady state, it is more efficient to compute the full time evolution of the system. This is comparable to imaginary time evolution algorithms to find the ground state of a closed many-body system. In our review, we will contrast the two approaches and address this distinction when discussing individual simulation methods in the main part of our review.

Investigating the full time evolution also offers the possibility to investigate interesting many-body effects during the relaxation dynamics. For instance, it is possible for open many-body systems exhibiting a quite trivial steady state, while the relaxation behavior is dominated by complex glassy quantum dynamics Olmos et al. (2012).

I.4 Differences to equilibrium problems

To find the steady state of an open quantum many-body system, it might first be tempting to take well established methods for ground state calculations for closed systems and try to adapt it to the open case. Unfortunately, this approach fails in many cases. For example, quantum Monte-Carlo methods that are highly successful for ground state calculations, require to rewrite the partition function of the quantum system to a corresponding classical system. However, for the steady state of an open system it is unclear a priori (and often incorrect Sieberer et al. (2013)) whether the steady state of the system is a thermal state that can be described in terms of a partition function. The same argument holds for density functional theory approaches trying to minimize the ground state energy; usually, the steady state of an open system is completely different from the ground state of the Hamiltonian. This is even true in the limit of infinitely weak dissipation, as the strength of the dissipation will predominantly control the relaxation rate rather than the properties of the steady state.

Some methods from the study of closed quantum systems out of equilibrium can be adopted to open systems; we will discuss these cases in detail. In general, the simulation of an open quantum system is computationally much harder than for a closed system due to the statistical nature of the state.

Additionally, one can benefit to some extent from the vast body of works commited to the study of classical non-equilibrium dynamics. For example, the importance of the symmetries of the open quantum many-body dynamics is equally important as in the classical case Hohenberg and Halperin (1977) and allows for the classification of dissipative phase transitions in terms of their universality classes.

I.5 Paradigmatic models

Within the analysis of ground state many-body problems, there is a number of particular models that have found especially wide interest and are often used as a first example to benchmark a numerical method. These models include, e.g., the Ising model in a transverse field, the Heisenberg model, and the Hubbard model (both bosonic and fermionic). A similar observation can be made about open quantum many-body problems, where these paradigmatic models are often derived from the corresponding ground state counterparts, i.e., the Hamiltonian dynamics is the same. However, adding dissipation to a closed many-body model can be done in different ways and can lead to drastically different results. In the following, we present and briefly discuss the two most promiment dissipative many-body models; we will provide a more detailed discussion in later sections when referring to particular numerical strategies to tackle them.

One of the most widely studied open many-body models in recent years is the transverse field Ising model with longitudinal dissipation Lee et al. (2011). Its Hamiltonian is of the form of the conventional Ising model, given in terms of Pauli matrices $\sigma_{\alpha}$ by

[TABLE]

where $h$ is the strength of the transverse field and $V$ accounts for the Ising interaction. The dissipation is incorporated in terms of jump operators of the form $c_{i}=\sqrt{\gamma}\sigma_{-}$, with $\gamma$ being the rate of dissipative flips from the spin up to the spin down state. An important aspect is that the dissipation breaks the $Z_{2}$ Ising symmetry of the Hamiltonian, i.e., the quantum master equation does not exhibit such a symmetry. The model is also relevant to ongoing experiments in the field of interacting Rydberg atoms Carr et al. (2013); Malossi et al. (2014).

Within a mean-field calculation Lee et al. (2011), the model is predicted to support a large range of $h$ values for which the system exhibits two stable steady states. We will discuss in later sections of our review how different numerical approaches address the question on the existence of such a bistable thermodynamic phase. According to mean-field theory, the bistable region ends in a critical point that belongs to the Ising universality class Marcuzzi et al. (2014).

Another important dissipative model is the driven-dissipative Bose-Hubbard model. While there are different ways to generalize the famous Bose-Hubbard model Fisher et al. (1989) to the dissipative case, the most commonly studied one involves a dissipative particle loss that can be countered by a coherent driving term Le Boité et al. (2013); Carusotto and Ciuti (2013). Its Hamiltonian is given by

[TABLE]

In this model, $J$ describes the hopping of bosons between sites, while the on-site interaction $U$ involves the square of the density operator $n_{i}=b^{\dagger}_{i}b_{i}$. Furthermore, $\Delta\omega$ is the chemical potential for the bosons, and $F$ describes the aforementioned coherent driving. Finally, the quantum jump operators capturing the loss of a single boson are given by $c_{i}=\sqrt{\gamma}b_{i}$. While the dissipation term also breaks the $U(1)$ symmetry of the conventional Bose-Hubbard model, here, the symmetry is already broken on the level of the Hamiltonian by the inclusion of the driving term $F$.

As with the dissipative Ising model, the driven-dissipative Bose-Hubbard model has a very intriguing mean-field phase diagram, where several islands of multistability occur in a way that is somewhat reminiscent of Mott lobes Le Boité et al. (2013), see Fig. 1. The stability of the mean-field solutions has been evaluated by considering density matrices of the form

[TABLE]

with $\rho_{i}^{MF}$ being the mean-field solution for the steady state. Expanding the quantum master equation up to first order in $\delta\rho_{i}$ allows to evaluate the stability by checking whether none of eigenvalues of the Liouvillian has a positive real part.

II Stochastic methods

Upon first glance, the computational complexity of an open quantum system in terms of the Hilbert space dimension $d$ appears to be at least $O(d^{2})$, as there are $O(d^{2})$ independent entries in the density matrix $\rho$. However, the density matrix at an initial time $t_{0}$ can be written as a statistical ensemble of pure states, $\rho(t_{0})=\sum_{i}p_{i}\mathopen{|}\psi_{i}(t_{0})\mathclose{\rangle}\mathopen{\langle}\psi_{i}(t_{0})\mathclose{|}$. Instead of propagating the entire density matrix, the key strategy is to propagate the individual pure states $\mathopen{|}\psi_{i}\mathclose{\rangle}$ to the time $t$ and then calculate observables according to

[TABLE]

The probability distribution $p_{i}$ can then be sampled using standard Monte-Carlo techniques, which is why the approach is often called the wave function Monte-Carlo method. In practice, the most common strategy is to start from an initial pure state $\mathopen{|}\psi_{0}\mathclose{\rangle}$ and perform $M=1/p_{i}$ numerical simulations. Since the trajectories $\mathopen{|}\psi_{i}\mathclose{\rangle}$ are independent from each other, the statistical error associated with the observable will behave as $\Delta O\sim 1/\sqrt{M}$. The entire computational cost will be $O(Md)$, which is considerably lower than $d^{2}$ already for quite modest system sizes. Importantly, the requirement to repeat the simulation $M$ times results in the simulation time being significantly longer than for a comparable closed quantum system. Depending on the observable, $M\approx 1000$ is a reasonable choice to get the statistical error down to a few percent. For spin $1/2$ systems, this essentially means that the system sizes that can be studied in an open system consist of $\log_{2}M\approx 10$ particles less than in a closed system.

The central question is now how can a single trajectory $\mathopen{|}\psi_{i}\mathclose{\rangle}$ be propagated such that the ensemble of all trajectories satisfies $\rho(t)=\sum_{i}p_{i}\mathopen{|}\psi_{i}(t)\mathclose{\rangle}\mathopen{\langle}\psi_{i}(t)\mathclose{|}$. One possiblity is to describe the evolution of the density operator in terms of a quantum state diffusion approach Gisin and Percival (1992); Percival (1998), in which the incoherent dynamics from the Lindblad operators is captured in terms of a stochastic Schrödinger equation,

[TABLE]

where the $dW_{j}$ refer to Wiener increments. The effective Hamiltonian $H_{eff}$ describes the drift of the state vector in the Hilbert space,

[TABLE]

The diffusion operators $M_{j}$ describe the random fluctuations arising from each associated jump operator $c_{j}$,

[TABLE]

This stochastic Schrödinger equation conserves the norm of the state vector and can be solved by standard techniques for stochastic differential equations.

An alternative strategy to propagate a single trajectory is the quantum jump method Dalibard et al. (1992); Dum et al. (1992); Mølmer et al. (1993); Plenio and Knight (1998). This approach has been recently reviewed extensively in Daley (2014), so we will only cover the basic strategy. Within the quantum jump method, the dynamics is split into to parts. First, the state $\mathopen{|}\psi_{i}\mathclose{\rangle}$ is propagated under an effective non-Hermitian Hamiltonian $H_{NH}$,

[TABLE]

Once the norm of the state drops below a previously drawn random number $r$, a quantum jump occurs. Which quantum jump occurs is drawn from the probability distribution

[TABLE]

with $\mathcal{N}$ being a normalization factor. While the high order integration of $H_{NH}$ is straightforward, a high order simulation of the quantum jumps requires a more subtle identification of the time the jump operator needs to be applied. For instance, the popular QuTiP library Johansson et al. (2012, 2013) uses a logarithmic secant method to numerically solve the equation $\langle\psi_{i}(t)|\psi_{i}(t)\rangle=r$ for the time $t$.

No matter which approach is used to propagate a single trajectory, the computations can be highly parallelized since the trajectories are independent from each other by construction. Doing so, it is possible to simulate open many-body spin $1/2$ models with up to 20 spins Raghunandan et al. (2018). The relatively small system sizes when compared to equilibrium problems demand the development of new data analysis techniques, e.g., concerning finite size scaling methods. One possiblity is to use anisotropic system sizes to obtain more data points for a reliable finite size scaling extrapolation. Close to a phase transition, the susceptibility $\chi$ of a system may be expressed as

[TABLE]

where $N$ is the number of particles and $\alpha$ is an exponent associated with the underlying phase transition Binder and Wang (1989). The reduced susceptibility $\tilde{\chi}$ is only a function of the anisotropy $\lambda$ of the system and can be determined by symmetry considerations as well as numerical data Raghunandan et al. (2018).

The wave-function Monte-Carlo method has been used to analyze the one-dimensional dissipative Ising model of Eq. (3) Ates et al. (2012); Hu et al. (2013). While these works have not found a bistable phase as predicted by mean-field theory, a significant increase in the spin correlations has been reported in the same region Hu et al. (2013). Additionally, finite size scaling of a similar two-dimensional model believed to lie in the same universality class as the dissipative Ising model has found evidence for a first order transition Raghunandan et al. (2018).

III Tensor network methods

III.1 One spatial dimension

We would, first, like to describe the important numerical techniques that has been developed for studying open quantum many-body systems using Matrix Product States (MPS) which is the one-dimensional ansatz of the tensor network (TN) family. MPS, by far, is the most successful and widely used ansatz in comparison to other ansatz of the Tensor network family, thanks to the success of the density matrix renormalization group (DMRG) White (1992, 1993) and related techniques Vidal (2003, 2004). Not only are its properties very well understood, contraction of MPS tensors can be done efficiently and exactly unlike the case for its higher dimensional counterparts Schuch et al. (2007); Haferkamp et al. (2018). For these reasons, MPS have been used extensively producing extremely accurate results, however, mostly in the context of ground state calculations of many-body systems Schollwöck (2005). For a detailed review on MPS and other tensor networks in general, we ask the readers to refer to Orus (2014); Eisert (2013); Schollwöck (2011); Orus (2018); Verstraete et al. (2008); Cirac and Verstraete (2009); Biamonte and Bergholm (2017). The application to open quantum systems, meanwhile, is more rare and there are only a few known approaches one can take for such systems. Not only are open systems more computationally challenging (since we need to deal with matrices in place of vectors for the pure states), there are also several intrinsic bottlenecks such as the positivity, hermiticity in the numerical optimization of the density matrix. Nevertheless, many of the ideas in the pure state formalism have been successfully applied in the context of open systems using the concept of Matrix Product Operators (MPOs) or Matrix Product Density Operators (MPDOs) Pirvu et al. (2010); Cirac et al. (2017). These appraches have also been used to study thermodynamic properties of 1d systems. We discuss them below.

In 2004, Verstraete et al. (2004) introduced the concept of MPDO which extended the MPS formalism from pure to mixed states. Let us recall that an MPS can be written in the following form

[TABLE]

where the $A$’s are matrices whose dimension is bounded by some fixed number $D$ (also called the bond dimension $\chi$) and $d$ is the physical dimension of the Hilbert space. An MPDO $\rho$ of $N$ $d$-level particles with $(D_{1},D_{2},\dots,D_{N})$-dimensional bonds is then defined as

[TABLE]

where $M_{k}^{s_{k},s_{k}^{\prime}}$ are $D_{k}^{2}\times D_{k+1}^{2}$ matrices that can be decomposed as

[TABLE]

where $d_{k}$ is at most $dD_{k}D_{k+1}$ and the matrices $A_{k}^{s,k}$ are of size $D_{k}\times D_{k+1}$. Such a construction of MPDOs automatically ensures the positivity of the reduced density matrix $\rho$. This is shown in Fig. 2.

This MPDO can be expressed in terms of a pure state MPS by defining it over a larger Hilbert space and using the concept of purification Nielsen and Chuang (2000). This can be done by associating an ancilla with a Hilbert space of dimension $d_{k}$ with each physical system. One can then choose an orthonormal basis $|s_{k},a_{k}\rangle$ for these physical and ancilla indices. The corresponding MPS for this system can be written as

[TABLE]

The MPDO $\rho$ can be obtained by tracing over the ancillas i.e. $\rho=\text{Tr}_{a}\left(|\Psi\rangle\langle\Psi|\right)$. This process is illustrated in Fig. 3.

The original $A_{k}$ matrices can be recovered from $M_{k}$ by doing some eigenvalue decomposition. To determine the evolution of a Hamiltonian of a mixed state in real and imaginary time, they simply simulated the evolution of the purification by updating the $A_{k}$ matrices using an iterative procedure similar to the standard DMRG in this technique. The purification could then be used to reconstruct the density operator at any time and compute the expectation values of the observables. Such a purification scheme can be used for mixed state evolution under dissipation as well as for thermal equilibrium and can be implemented irrespective of periodic or open boundary conditions, finite or infinite systems. The main source of errors in this procedure, like most other TN techniques are (i) Trotter error and (ii) truncation error. Such an approach with ancillas was also applied to study the thermodynamic properties of several spin chains in Ref. Feiguin and White (2005). Although the MPDOs in Ref. Verstraete et al. (2004) are positive by construction, it was shown in Ref. las Cuevas et al. (2013) that such a MPDO descriptions of mixed states are not exactly equivalent to the one obtained using local purification schemes. In particular, it was shown that the bond dimension of the locally purified MPS $D^{\prime}$ is not upper bounded by the bond dimension of the MPDO $D$. In fact, the local purification techniques can be much more costly than the MPDO form itself. Thus, the authors concluded that a description of mixed states which is both efficient and locally positive semi definite does not exist and that one can only make approximations.

Around the same time, Zwolak and Vidal (2004) proposed another technique to study the mixed state dynamics in one dimensional lattice systems. Their technique, which is also based on MPS, used the Time Evolving Block Decimation (TEBD) to simulate the real time Markovian dynamics given by a master equation with nearest-neighbor couplings. At the heart of this algorithm, lies the concept of ‘Choi isomorphism’. It is more of a mathematical trick and it states that one can rewrite the coefficients of a matrix as those of a vector. In other words, this is simply turning a bra index into a ket index for a density matrix (understanding the coefficients of $\rho$ as those of a vectorized density matrix denoted by $|\rho\rangle_{\sharp}$). And in the language of TN diagrams, it can be regarded as reshaping one of the legs and gluing it with the other (Fig. 4).

Once vectorized, $|\rho\rangle_{\sharp}$ now lives in the $n$-fold tensor product of $\mathbb{C}_{d^{2}}$ and the master equation can be written in the vector form. The mixed state will now look like as follows

[TABLE]

where $|i_{l}\rangle_{\sharp}$ is an orthonormal basis of $\mathbb{C}_{d^{2}}$ for site $l$. Further assuming that the Liouvillian superoperator $\mathcal{L}$ can be decomposed into terms involving at most nearest-neighbors i.e. $\mathcal{L}[\rho]=\sum_{l}\mathcal{L}_{l,l+1}[\rho]$, one could in principle use the usual TEBD algorithm to solve Eq. 2 by starting from some initial Matrix Product Operator (MPO) (shown in in the left side of Fig. 4). This was the basic idea behind the technique in Ref. Zwolak and Vidal (2004). One of the first applications of this technique was the study of the driven-dissipative Bose-Hubbard model in the context of optical resonators Hartmann (2010). More detailed explanation of this vectorization process will be explained later when we discuss the case for higher dimensional systems. Although the technique proved to be extremely simple and efficient, the issue of positivity still remained at large. In fact, checking the positivity of a reduced density matrix is known to be a very hard problem in physics Kliesch et al. (2014).

Another approach was taken in Ref. Werner et al. (2016) to solve the problem of positivity. In this approach, instead of expressing $\rho$ directly as an MPO, at every stage of the algorithm, $\rho$ was kept in its locally purified $\rho=XX^{\dagger}$, where the purification operator $X$ is decomposed as a variational tensor network.

[TABLE]

where $1\leq s_{l}\leq d$ , $1\leq r_{l}\leq K$ and $1\leq m_{l}\leq D$. $A^{[l]}$ are rank-four tensors with physical dimension $d$, bond dimension $D$ and Kraus dimension $K$. Then, a technique similar to the usual TEBD was used to update the tensors. Such an approach never required to contract the two TN layers ($X$ and $X^{\dagger}$) together, thereby ensuring positivity at all times during the evolution. The technique also provided more control of the approximation error with respect to the trace norm.

In Ref. Cui et al. (2015), a very interesting and different approach based on MPO was taken for finding the steady states of dissipative 1D systems governed by the master equation of the Lindbladian form, $\frac{d\rho}{dt}=\mathcal{L}[\rho]$, where $\mathcal{L}$ is the Liouvillian superoperator. In this technique, instead of doing the full real time evolution of the Liouvillian, they proposed a variational method that searches for the null eigenvector of $\mathcal{L}$ which is, by definition, the steady state of the master equation in the Lindbladian form. Their results were based on the principle that if $\rho_{s}$ is the steady state of the Lindbladian master equation satisfying $\hat{\mathcal{L}}|\rho_{s}\rangle_{\sharp}=0$, then $|\rho_{s}\rangle_{\sharp}$ will also be the ground state of the non-local Hamiltonian $\hat{\mathcal{L}}^{\dagger}\hat{\mathcal{L}}$ (since it is Hermitian and positive semi-definite) where $|\rho_{s}\rangle_{\sharp}$ is the vectorized form of the steady state density matrix. Then using a variational algorithm, they directly targeted the ground state of $\hat{\mathcal{L}}^{\dagger}\hat{\mathcal{L}}$ to find the steady state of the Lindbladian master equation for a finite chain. One of the reasons why directly targeting the ground state of $\hat{\mathcal{L}}^{\dagger}\hat{\mathcal{L}}$ might be advantageous is that unlike imaginary time evolution, where the sequence of states visited by the algorithm is unimportant, the simulation of a master equation requires us to follow real time evolution. Therefore, if there are errors in the intermediate states visited by the algorithm, it may lead to problems in the convergence of our steady state. For example, some of the intermediate states may require large bond dimensions of the MPO although it is known that the final steady state can be well-represented by an MPO of small bond dimensions Cai and Barthel (2013); Bonnes et al. (2014). Also, one doesn’t need to worry about the large entanglement growth of real time evolution. A very similar approach was taken in Ref. Mascarenhas et al. (2015) where the algorithm instead of doing a time-evolution, searched for the null eigenvalue of the Liouvillian superoperator $\mathcal{L}$ by sweeping along the system. Their method claimed to work even in the weakly dissipative regime by slowly tuning the dissipation rates along the sweeps. However, it needs to be noted that such techniques, while advantageous numerically, cannot be used for obtaining the transient states.

In another paper Gangat et al. (2017), this idea was applied to infinite 1D systems (i.e. the thermodynamic limit) using a hybrid technique of both imaginary and real time evolution. They took a local auxiliary Hamiltonian $\mathcal{H}$ whose ground state is a good approximation to the ground state of the nonlocal Hamiltonian $\hat{\mathcal{L}}^{\dagger}\hat{\mathcal{L}}$ by taking its $k^{\text{th}}$ root as

[TABLE]

where $\hat{\mathcal{L}}=\sum_{r\in\mathbb{Z}}\hat{\mathcal{L}}_{r}$ since $\hat{\mathcal{L}}$ is a translationally invariant local operator. The $k^{\text{th}}$ root was taken in order to yield faster convergence. The idea is that if the gap between the two lowest eigenvalues of $\hat{\mathcal{L}}^{\dagger}_{r}\hat{\mathcal{L}_{r}}$ is less than one, then $k>1$ will increase the gap since $\hat{\mathcal{L}}^{\dagger}_{r}\hat{\mathcal{L}_{r}}$ is positive semi definite, thereby achieving faster convergence to the ground state. The authors then performed a real time evolution to obtain a more accurate steady state. In quick summary, the main steps of the algorithm is:

1. (i)

Imaginary time evolution of the auxiliary Hamiltonian $\mathcal{H}$ starting from some vectorized initial density matrx $|\rho_{0}\rangle$

[TABLE] 2. (ii)

Real time evolution of the Liouvillian superoperator starting from $|\rho_{G}\rangle$

[TABLE]

where $|\rho_{S}\rangle$ is the desired steady state of the Liouvillian master equation. Imaginary time evolution in step(i) ensures that one does not pass through highly entangled transient regime. Step(ii) increases the accuracy of the stationary state since $|\rho_{G}\rangle$ is the ground state of $\mathcal{H}$ which is a truncated approximation of the non local Hamiltonian $\hat{\mathcal{L}}^{\dagger}\hat{\mathcal{L}}$.

It is worth mentioning that in one spatial dimension, many of the above techniques and their combinations have been used for studying not only other important disspative models Pižorn (2013); Höning et al. (2012); Mascarenhas et al. (2015); Carollo et al. (2019), including the dissipative Ising model of Eq. (3) Höning et al. (2013); Mendoza-Arenas et al. (2016), but also in dissipative preparation of topologically ordered materials Iemini et al. (2016) as well as in the energy transport Guo et al. (2015). Very recently, MPO based techniques have been applied to study vibronic states which extends the application to quantum biology and organic photo voltaics Somoza et al. (2019). and also to study the dynamics of photonic circuits with time delays and quantum feedback Pichler and Zoller (2016). We do not discuss the later two works due to the non-Markovian nature of the problem which is beyond the scope of this review. Similar MPS based techniques that go beyond the Lindblad master equation Xu et al. (2019) or the markovian approximation Guo et al. (2018) are also not discussed here.

III.2 Extensions to higher dimensions

Unlike the case for 1d, the generalization of MPS in higher dimensions, also known as Projected Entangled Pair States (PEPS) or Tensor Product States (TPS) comes with some serious limitations and there are still many open problems Cirac et al. (2019). Not only do the PEPS algorithm require serious programming effort, exact contratcion of PEPS is known to be a mathematically hard problem Schuch et al. (2007); Haferkamp et al. (2018). To achieve this, one requires additional PEPS contraction algorithms Orús (2012); Orús and Vidal (2009); Jordan et al. (2008) that are nevertheless known to give very accurate results, in particular, for gapped systems. Even for critical systems with algebraically decaying correlations, the PEPS contraction schemes are known to provide reasonably accurate results with sufficiently high bond dimension of the environment Orús and Vidal (2009). In fact, recently, techniques have been introduced to capture the infinite correlation length of 2D critical systems using iPEPS based on finite correlation length scaling Corboz et al. (2018); Rader and Läuchli (2018). Thus, despite the higher requirement of numerical dedications and limitations, PEPS algorithms are becoming state of the art numerical tools for strongly correlated two dimensional systems. Recently, PEPS provided the best variational energy for the 2D Hubbard model Corboz (2016a), have offered several new insights on paradigmatic models and real materials in the lab Liao et al. (2017); Corboz and Mila (2014); Matsuda et al. (2013); Kshetrimayum et al. (2019). The successes of PEPS so far, however, is mostly confined to ground state calculations and partially to thermal states Czarnik et al. (2012); Czarnik and Dziarmaga (2015); Czarnik et al. (2016); Kshetrimayum et al. (2019); Dai et al. (2017) using the concept of Projected Entangled Pair Operators (PEPOs) or Tensor Product Operators (TPOs), which we will discuss in more detail later, and, more recently, to time evolution Hubig and Cirac (2019); Czarnik et al. (2019); Kshetrimayum et al. (2019). For the context of open dissipative quantum system, so far there is only one known approach using PEPS Kshetrimayum et al. (2017) and another one using a Corner Space Renormalization method Finazzi et al. (2015). We describe them below. We will also discuss briefly other potential implementation techniques and possible issues while using PEPS formalism in particular for such open systems.

The Corner Space Renormalization method Finazzi et al. (2015) solves the master equation in a corner of the Hilbert space through an iterative procedure. It starts by finding the steady state density matrix for small lattice systems (say $\rho^{A}$ and $\rho^{B}$ for systems $A$ and $B$ respectively). This can be done by a brute force integration of the master equation since the system size is very small. The steady state density matrices can be diagonalized and written as

[TABLE]

where the states $|\phi_{i}^{A}\rangle$ form an orthonormal basis for $\mathcal{H}_{A}$ (the Hilbert space corresponding to system $A$) and $p_{i}^{A}$ are the corresponding probabilities. Similar notation follows for system $B$. The two systems are then merged and the $\chi$ most probable product states spanning the so-called corner space are selected i.e. we only keep the subspace generated by the orthonormal basis $\{|\phi_{i1}^{A}\rangle|\phi_{i^{\prime}1}^{B}\rangle,|\phi_{i2}^{A}\rangle|\phi_{i^{\prime}2}^{B}\rangle,\dots,|\phi_{i{\chi}}^{A}\rangle|\phi_{i^{\prime}{\chi}}^{B}\rangle\}$ where the product of the probabilities of the two systems are arranged in decreasing order of magnitude. In this way, we only keep the $\chi$ most probable pair of states. The steady state of the density matrix in this corner space can be determined either by direct numerical integration in time (for small $\chi$) or by using a stochastic wave function Monte Carlo algorithms for large $\chi$, see Sec. II. One can then increase the size of the corner $\chi$ until convergence in some observables is reached. Larger systems can be simulated by merging more systems as we discussed in the initial steps. A simplified summary of the steps involved is shown in Fig. 5.

The proposed CSR method was used to study the driven-dissipative Bose-Hubbard model in 2D in both periodic and open boundary conditions for system sizes up to $16\times 16$ lattice sites. The technique has also been used to study the critical Heisenberg model Rota et al. (2017) for system size up to $6\times 6$ lattice sites and more recently the critical regime in the Bose Hubbard model Rota et al. (2019) for up to $8\times 8$ lattices. The size of the lattice that can be simulated using this technique depends on the entanglement of the steady state. Even if not obvious at first sight, the structure of the density operator generated by the Corner Space Renormalization method amounts to that of a Tree Tensor Network Shi et al. (2006). As such, this particular method, even if understood in terms of TNs, is tailored to driven-dissipative systems of finite size. For generalizing it to the thermodynamic limit or for non-driven non-dissipative systems, one needs to use more general TN techniques. We will discuss one such technique which we developed recently below.

In Kshetrimayum et al. (2017), we make use of the concept of PEPO by vectorizing them. PEPOs are simply the operator version of PEPS, in the same way that an MPO is the operator version of MPS for the 1d case. Hence, PEPOs are used to represent mixed states $\rho$ in 2D, even beyond dissipative systems, e.g., for thermal states Czarnik et al. (2012); Kshetrimayum et al. (2019); Czarnik and Dziarmaga (2015); Czarnik et al. (2016). As mentioned before, such a construction of density matrices using PEPOs does not automatically guarantee the positivity of the density matrix. However, for simulations targeting the steady states, this lack of exact positivity is not a bottleneck if the fixed point is not very highly entangled. For the moment, we will restrict our discussion to this case. Once we have our PEPO, we vectorize it i.e. rewrite the coefficients of the PEPO as a PEPS (also called Choi’s isomorphism). Once vectorized, the PEPO $\rho$ can be understood as a PEPS of physical dimension $d^{2}$ and bond dimension $D$ (now called $|\rho\rangle_{\sharp}$), as shown also in Fig. 6.

The vectorized form of the Lindblad master equation Eq. (2) can be written as

[TABLE]

where the vectorized Liouvillian operator is given by

[TABLE]

$H$ is the Hamiltonian of the system and $\mathbb{I}$ corresond to the identity operator. $L_{\mu}$ and $L_{\mu}^{\dagger}$ correspond to the on-site Lindblad/jump operators, responsible for dissipation. The tensor product $\otimes$ separates the operator acting on the ket and bra index of $\rho$ before the vectorization. When the vectorized Liouvillian superoperator $\mathcal{L}_{\sharp}$ is independent of time, Eq. 23 can be integrated as

[TABLE]

where $|\rho(0)\rangle_{\sharp}$ is some vectorized initial density matrix, written as a PEPS. In the limit of $t\rightarrow\infty$, we obtain the non-equilibirum steady state (NESS) as the fixed point of the master equation which we denote by $|\rho_{s}\rangle_{\sharp}$. From Eq. 23, it is also obvious that $|\rho_{s}\rangle_{\sharp}$ is the right eigen vector of $\mathcal{L}_{\sharp}$ with zero eigen value so that

[TABLE]

For a Liouvillian $\mathcal{L}$ consisting of local terms say $\mathcal{L}[\rho]=\sum_{\langle i,j\rangle}\mathcal{L}^{[i,j]}[\rho]$, the vectorized form of the Lindblad equation Eq. (2) yields a parallelism with the calculation of ground states of local Hamiltonians by imaginary-time evolution, which we detail in Table 1.

Given the parallelism above, it is clear that one can adapt, at least in principle, the methods to compute imaginary time evolution of a pure state as generated by local Hamiltonians, to compute also the real time evolution of a mixed state as generated by local Liouvillians. This was, in fact, the approach taken in Ref.Zwolak and Vidal (2004) for finite-size 1d systems, using Matrix Product Operators (MPO) to describe the 1d reduced density matrix, and proceeding as in the Time-Evolving Block Decimation (TEBD) algorithm for ground states of 1d local Hamiltonians Vidal (2003, 2004, 2007); Orús and Vidal (2008) as we have discussed previously. In Kshetrimayum et al. (2017), we extended this implementation for the case of 2D systems using the concept of PEPO with physical dimension $d$ and bond dimension $D$, see Fig.6. For the case of an infinite-size 2d system, this setting is actually equivalent to that of the infinite-PEPS algorithm (iPEPS) to compute ground states of local Hamiltonians in 2d in the thermodynamic limit. Thus, in principle, one can use the full machinery of iPEPS to tackle as well the problem of 2d dissipation and steady states.

There seems to be, however, one problem with this idea: unlike in imaginary-time evolution, we are now dealing with real time. In the master equation, part of the evolution is generated by a Hamiltonian $H$, and part by the Lindblad operators $L_{\mu}$. The Hamiltonian part corresponds actually to a unitary “Schrödinger-like” evolution in real time, which typically increases the “operator-entanglement” in $|\rho\rangle_{\sharp}$, up to a point where it may be too large to handle for a TN representation (e.g., 1d MPO or 2d PEPO) with a reasonable bond dimension. In 1d this is the reason why the simulations of master equations are only valid for a finite amount of time. In 2d, simple numerical experiments indicate that in a typical simulation the growth of entanglement is even faster than in 1d. Luckily, this is not a dead-end: if the dissipation is strong compared to the rate of entanglement growth, then the evolution drives the system into the steady state before hitting a large-entanglement region. In fact, even if there is too much entanglement for the TN at intermediate times, the dissipation may still drive the evolution towards a good approximation of the correct steady state. In short, dissipation limits the growth of entanglement if the fixed point attractor is strong enough. This can be verified numerically by plotting the operator entanglement entropy for different dissipation strengths as it flows into the NESS. This is shown in Fig. 7. Details on how to compute this quantity can be found in Kshetrimayum et al. (2017); Kshetrimayum (2017).

Hence, one can apply the iPEPS machinery to compute the time evolution in 2d with a local Liouvillian $\mathcal{L}$ and some initial state. This procedure was used to investigate the dissipative Ising and the XYZ model, confirming and offering several insights that were inaccesible before using mean field and other techniques. For example, for the dissipative Ising model of Eq. (3), given by the Hamiltonian $H=\frac{V}{4}\sum_{\langle i,j\rangle}\sigma_{z}^{[i]}\sigma_{z}^{[j]}+\frac{h_{x}}{2}\sum_{i}\sigma_{x}^{[i]}$ and Lindblad operators $L_{\mu}=\sqrt{\gamma}\sigma_{-}^{[\mu]}$. The phase diagram is controversial with some papers suggesting the existence of a bistable steady state Lee et al. (2011); Marcuzzi et al. (2014) and others supporting a first order transition Weimer (2015a, b); Maghrebi and Gorshkov (2016). Our technique found bistability for low bond dimensions of the PEPO ($D=1,2$) which was replaced by a first order transition for higher $D$s, thus confirming that the bistability is an artifact of mean field. This is shown in Figure 8.

Furthermore, some studies suggested the existence of an antiferromagnetic region in the presence of the transverse field $h_{x}$ Lee et al. (2011); Weimer (2015a). Once again, while our technique found evidence of such an AF region, it eventually shrank with increasing bond dimension until it finally disappears for large enough $D$. Results of the dissipative Ising model have been reproduced independently using a different update scheme Czarnik et al. (2019) compared to the one used in Kshetrimayum et al. (2017). While the technique also employed vectorization along with iPEPS, the update scheme is based on maximizing the fidelity between two consecutive steps of the update of the iPEPS tensors. For the case of the dissipative Heisenberg model Lee et al. (2013) with the Hamiltonian

[TABLE]

and the same Lindblad operators as before, our studies found no phenomenon of re-emergence in the phase diagram, confirming a prediction by studies using cluster mean-field approaches Jin et al. (2016).

To the best of our knowledge, we have discussed most of the state-of-the-art numerical techniques based on TN for the study of open quantum many-body systems in both one and two spatial dimensions. We would, now, like to discuss some of the possible ideas that could be helpful in improving the existing algorithm and possible new implementation techniques specially in 2d. First of all, we remark that the 2d algorithm suggested above does not guarantee the positivity of the density matrices. This problem can be solved by starting from an initial state that is positive by construction, for example taking the product of two PEPOs which are the conjugate of each other ($A$ and $A^{*}$). One can then think about using a positivity preserving algorithm such as the one in Werner et al. (2016). Such an algorithm will ensure the positivity of the density matrix at all times of the evolution. We can call this initial density matrix as Projected Entangled Pair Density Operator (PEPDO) as shown in Figure 9.

While such an approach may avoid the problem of negative eigenvalues of the density matrix, in practice, it may require very high bond dimension of the PEPDO and one, therefore, needs to consider the practical aspect of the implementation.

The other possibility would be to target the ground state of the Hermitian and positive semidefinite operator $\mathcal{L}_{\sharp}^{\dagger}\mathcal{L}_{\sharp}$. This ground state could be computed, e.g., by an imaginary time evolution. However, there are two major hurdles associated with this approach. First, the crossed products in $\mathcal{L}_{\sharp}^{\dagger}\mathcal{L}_{\sharp}$ are non-local, and therefore the usual algorithms for time evolution are difficult to implement unless one introduces extra approximations in the range of the crossed terms. Another option is to approximate the ground state variationally, e.g., via the Density Matrix Renormalization Group White (1992, 1993) in 1d, or variational PEPS in 2d Corboz (2016b). In the thermodynamic limit, however, this approach does not look very promising because of the non-locality of $\mathcal{L}_{\sharp}^{\dagger}\mathcal{L}_{\sharp}$ mentioned before. In any case, one could always represent this operator as a PEPO (in 2d), which would simplify some of the calculations, but at the cost of introducing a very large bond dimension in the representation of $\mathcal{L}_{\sharp}^{\dagger}\mathcal{L}_{\sharp}$. For instance, if a typical PEPO bond dimension for $\mathcal{L}_{\sharp}$ is $\sim 4$, then for $\mathcal{L}_{\sharp}^{\dagger}\mathcal{L}_{\sharp}$ it is $\sim 16$, which in 2d implies extremely slow calculations. Another option would be to target the variational minimization of the real part for the expectation value of $\mathcal{L}$.

IV Variational methods

Variational techniques are often very powerful tools to analyze quantum many-body systems, as demonstrated by the successes of density functional theory Kohn (1999) and matrix product state approaches Schollwöck (2011) for ground state problems. As we will discuss in this section, variational methods can also be successfully applied to open quantum many-body systems.

IV.1 The variational principle for open quantum systems

Variational methods generically consist of two steps. The first step is a paramterization of the state of the system in terms of a set of variational parameters $\{\alpha_{i}\}$. For open quantum system, it is convenient to parameterize the density matrix, i.e., $\rho=\rho(\{\alpha_{i}\})$, although parameterizations based on statistical ensembles of pure states are also possible Transchel et al. (2014). The second step is to identify a suitable functional that can be optimized by tuning the variational parameters. For open quantum systems, it is very natural to apply a variational principle to find the steady state of the quantum master equation, which can be found by solving the equation $\dot{\rho}=0$. Importantly, the exact steady state can no longer be determined after the variational parametrization, as the steady state will generically lie outside the variational manifold. Hence, the best possible option is to find the variational parameters that will minimize the functional $||\mathcal{L}\rho||$ for a suitable norm Weimer (2015b).

The correct norm for the variational optimization can be identified as the trace norm $||\mathcal{L}\rho||=\text{Tr}\{|\dot{\rho}|\}$, i.e., the sum of the absolute values of the eigenvalues of $\dot{\rho}$ Weimer (2015b). This choice can be motivated on different grounds. First, the trace distance being the natural distance measure for density matrices Nielsen and Chuang (2000) is highly suggestive of the trace norm being the natural norm for the tangent space $\dot{\rho}$. This can be formalized in the sense that the trace norm describes an optimal measurement to distinguish $\dot{\rho}$ from the zero matrix Gilchrist et al. (2005). A second way to motivate the trace norm is to consider classes of possible alternatives. It can be shown that all Schatten $p$-norms of the form $(|\dot{\rho}|^{p})^{1/p}$ are inherently biased towards the maximally mixed state for all values of $p>1$ Weimer (2015b). Since functionals with $p<1$ do not constitute proper norms, this leaves the trace norm as the only valid choice. One can also understand the variational principle as a direct solution of the overdetermined steady state equation $\mathcal{L}\rho=0$ in terms of a trace norm minimization.

In general, the evaluation of the variational functional is still an exponentially hard problem, as the computation of the trace norm requires the diagonalization of the matrix $\dot{\rho}$. However, it is possible to construct upper bounds to the variational norm that retain the variational character Weimer (2015b) and appear to introduce only small quantitative deviations even close to phase transitions Weimer (2015a). The upper bound depends on the variational manifold and its tangent space, i.e., the degree of additional correlations that can be build up by applying the Liouvillian to states within the variational manifold. For example, for a variational class of product states of the form $\rho=\prod_{i}\rho_{i}$, the upper bound $D$ can be given as

[TABLE]

where $\mathcal{T}$ contains pairs of sites that are connected to each other by the Liouvillian Weimer (2015b).

The variational principle has been applied to find the steady states of the dissipative Ising and Bose-Hubbard models introduced in Eqs. (3) and (4), respectively Weimer (2015b, a), as well as dissipative Ising models including a $\mathbb{Z}_{2}$ symmetry Overbeck et al. (2017), purely dissipative Heisenberg models Weimer (2017), dissipative Rydberg gases Weimer (2015a), dissipative ensembles of nitrogen-vacancy centers Raghunandan et al. (2018), entanglement generation in cavity QED arrays Lammers et al. (2016), and dissipative Fermi-Hubbard models Kaczmarczyk et al. (2016). In the latter case, the study of fermionic models was realized by employing a two-dimensional Jordan-Wigner transformation, where the appearance of nonlocal Wigner strings was ruled out by the choice of the variational manifold.

In the case where the steady state of the system is close to criticality, it is possible to construct a dissipative Ginzburg-Landau theory based on the variational principle Overbeck et al. (2017). The essential step is to perform a series expansion of the variational norm of Eq. 28 in terms of an order parameter field $\phi(x)$ and its spatial gradient $\nabla\phi(x)$, leading to

[TABLE]

All the coefficients $v_{n}$ and $u_{n}$ can be calculated from the microscopic quantum master equation. The series can be truncated at low orders of $m$ and $n$, as higher order terms are irrelevant close to criticality. In the case of steady states with thermal statistics due to the presence of a dynamical symmetry Sieberer et al. (2013), it is possible to construct a Ginzburg-Landau-Wilson functional integral for an effective partition function Hohenberg and Krekhov (2015), given by

[TABLE]

Here, the effective inverse temperature $\beta_{\text{eff}}$ can be derived from the $u_{0}$ coeffecient, as this coefficient captures the strength of fluctuations beyond a spatially homogeneous order parameter field Overbeck et al. (2017). The subsequent statistical field theory of Eq. (30) can then be analyzed using standard techniques such as the perturbative renormalization group.

Finally, the variational principle can also be extended towards the full time evolution of open quantum systems Overbeck and Weimer (2016), following very similar ideas discussed in the context of the time-dependent variational principle Kraus and Osborne (2012); Transchel et al. (2014). There, the variational functional is replaced by a variational integration of the quantum master equation for small time steps $\tau$. For example, in the lowest order Euler approximation, it is given by

[TABLE]

where $\rho(t+\tau)$ is the density matrix containing the variational parameters. Higher-order schemes exist as well, but constructing an upper bound similar to Eq. (28) requires to consider higher-order correlations due to multiple applications of the Liouvillian to the density matrix. A good compromise is the implicit midpoint method, which is exact up to second order in $\tau$ while only requiring a single application of the Liouvillian Overbeck and Weimer (2016).

IV.2 Comparison with mean-field methods

For equilibrium problems, the variational method based on product states is exactly equivalent to a mean-field decoupling of the interaction terms. Remarkably, this is not the case for open quantum systems. Within the mean-field approach to open systems Tomadin et al. (2010); Diehl et al. (2010a), a set of effective single site master equations is considered that is obtained by tracing out the rest of the system. For the $i$th site, the mean-field master equation reads

[TABLE]

where $H_{i}^{MF}$ and $\mathcal{D}_{i}$ are the mean-field Hamiltonian and the mean-field dissipators, respectively. This set of equation is then solved self-consistently, while for translationally invariant systems it is often sufficient to consider an effective single site problem.

Due to the nonlinear structure of the mean-field equations of motion, it is possible to have two or more independent solutions for the steady state Lee et al. (2011), see Fig. 10. This also occurs within mean-field theory for equilibrium systems close to first order transitions. However, there one can always resort to the free energy, which has to be minimal in thermal equilbrium. Unless one invokes the variational principle, one cannot decide which of the solutions of mean-field theory are stable and which ones are not. Interestingly, the solution according to the variational principle and mean-field theory are only identical in the limit of infinite dimensions, where both approaches become exact Weimer (2015b).

Mean-field theory predicts bistability for a wide range of models, including the dissipative Ising model Lee et al. (2011); Lee and Cross (2012); Marcuzzi et al. (2014) or extended spin models Parmee and Cooper (2018), as well as driven-dissipative Bose-Hubbard models Le Boité et al. (2013); Jin et al. (2013); Le Boité et al. (2014); Mertz et al. (2016). So far, mean-field bistability has been found in the absence of symmetries in the underlying master equation, i.e., the two solutions are not connected by a symmetry transformation. These properties have led to speculation that bistability could be a genuine nonequilibrium phase, which has stimulated several investigations whether this could indeed be the case. However, the results of these investigation have all been negative so far. Specifically, the variational principle predicts that bistability is replaced by a first order transition both in the dissipative Ising and in the driven-dissipative Bose-Hubbard model Weimer (2015b). For the dissipative Ising model, the existence of the first order transition has been confirmed in tensor network simulations, where bistability is found for low bond dimensions, but a first order transition appears for higher bond dimensions Kshetrimayum et al. (2017), see Sec. III.2. In the case of the driven-dissipative Bose-Hubbard model, the first order transition has also been found in a field-theoretic treatment based on the Keldysh formalism Maghrebi and Gorshkov (2016), again confirming the variational prediction. These results underscore that the conventional argument of mean field theory becoming qualitiatively correct if the spatial dimension becomes large enough appears to be incorrect for open quantum systems. On the other hand, this argument seems to be much more justified when applied to the variational principle (especially, when considering the connection to equilibrium statistical physics through the existence of the dissipative Ginzburg-Landau theory of Eq. (29)), however, even there one may have possible counterexamples Mesterházy and Hebenstreit (2017), which are not yet fully understood.

Nevertheless, these findings do not rule out genuine bistability in open quantum systems per se, but only that such mean-field results need to be taken with caution. Classical models exhibiting extended coexistence regions Muñoz et al. (2005) might still exhibit bistability after including quantum fluctuations. The situation is similar when it comes to limit cycles of open quantum many-body systems Chan et al. (2015), which have been predicted to exist in sufficiently high-dimensional systems Owen et al. (2018).

One systematic extension of mean-field theory is cluster mean-field theory, where the trace in Eq. (32) is not carried out over all but one site but results in a larger cluster that has again to be solved self-consistently Jin et al. (2016). This strategy is in close analogy to the cluster mean-field theory for statistical mechanics and ground state problems Bethe (1935); Oguchi (1955). Essentially, cluster mean-field approaches treat the short-range physics more accurately than bare mean-field theory, leading to better quantitative estimates for phase transitions. However, the qualitative limitations of bare mean-field theory remains, as these are resulting of long-range fluctuations in the system. For open quantum many-body models, cluster mean-field theory has been used to calculate the phase diagram of the dissipative Heisenberg model given by Eq. 27 Jin et al. (2016) and dissipative Ising models with and without a $\mathbb{Z}_{2}$ symmetry Jin et al. (2018).

Finally, it is also possible to systematically go beyond the mean-field approximation using open system dynamical mean-field theory (DMFT). DMFT is a mapping of a many-body lattice model onto a single impurity problem that has to be solved in a self-consistent way Georges et al. (1996). Within DMFT, the approach is to start with an effective dynamical Green’s function $\mathcal{G}_{0}$, which serves as a time-dependent version of a mean-field coupling. Considering the Fermi-Hubbard model as an example, $\mathcal{G}_{0}$ can used to express the effective action of a single site as

[TABLE]

where $f_{\sigma}$ annihilates a fermion with spin $\sigma$, $\beta$ is the inverse temperature, and $U$ is the on-site interaction. The central idea of DMFT is to consider a self-consistent solution that repoduces the dynamical Green’s function $\mathcal{G}_{0}$. This constraint is satisfied by the solution to the DMFT equations for the local Green’s function $G_{0}$, the dynamical Green’s function $\mathcal{G}_{0}$, and the self-energy $\Sigma$ evaluated at the Matsubara frequencies $\omega_{n}=(2n+1)\pi/\beta$,

[TABLE]

where $\mu$ is the chemical potential and $N(\epsilon)$ is the density of states Kollar (2011). The first step in bringing DMFT to open systems has been to use effective Lindblad master equations have to describe quantum transport in closed quantum systems using DMFT Arrigoni et al. (2013); Titvinidze et al. (2015, 2016). Recently, this approach has been extended to the case where already the initial many-body problem describes an open quantum system Panas et al. (2019).

A different method to systematically extend mean-field theory is to use projection operator methods. The central idea is to consider a single site of the many-body problem, with the rest of the system forming a non-Markovian environment. This non-Markovian master equation is then solved using standard projection operator techniques such as the Nakajima-Zwanzig method or the time-convolutionless master equation Breuer and Petruccione (2002). Initially, this approach has been used to describe the relaxation dynamics of local observables in a closed quantum system Weimer et al. (2008), which has later been extended to the Lindblad dynamics of open systems Degenfeld-Schonburg and Hartmann (2014). There, the initial step is to introduce corrections $\Delta\mathcal{L}$ to the mean-field Liouvillian $\mathcal{L}_{\text{MF}}$ are introduced according to

[TABLE]

The projection $\mathcal{P}$ removes all correlations and projects the system onto a product state, i.e.,

[TABLE]

If the initial state at time $t_{0}$ is also a product state, the projected Lindblad master equation may be formally written as

[TABLE]

where the generator $\mathcal{K}$ has been introduced Degenfeld-Schonburg and Hartmann (2014). The generator $\mathcal{K}$ may then be expanded in terms of a power series of the beyond mean-field corrections $\Delta\mathcal{L}$. This projection operator approach has been used to investigate both dissipative XY models Degenfeld-Schonburg and Hartmann (2014) and the dissipative Heisenberg model Owen et al. (2018). Remarkably, in the latter case a limit cycle behavior has been reported, which for sufficiently large spatial dimensions also survives under inclusion of the terms beyond mean-field. Consequently, it would be interesting to learn whether the projection operator approach is also capable to correctly identify the replacement of mean-field bistability by a first order transition in the dissipative Ising model.

IV.3 Variational tensor network methods

Given the successes of tensor network methods discussed in Sec. III, it appears natural to combine them with variational methods for the study of open quantum many-body systems. However, the main challenge is that the natural trace norm for constructing the variational principle cannot be calculated efficiently in a tensor network respresentation. This has led to the use of different norms as possible alternatives Cui et al. (2015); Mascarenhas et al. (2015), see Sec. III.1.

On the one hand, the choice of the norm is not really relevant if the value of the norm is very low (i.e., comparable to the machine precision of the numerical simulation), as then the solution is almost exact from any point of view. On the other hand, choosing a non-natural norm is a potential source of errors that is not under control of the variational algorithm. In practice, this difficulty will mostly manifest itself for higher-dimensional problems, as there the bond dimensions that can be reached are severly constrained by the compuational resources Kshetrimayum et al. (2017). But even for one-dimensional systems, there are computationally challenging problems involving long relaxation times Carollo et al. (2019), where an arbitrarily low variational norm might not be reachable.

A way out of this problem can be realized by representing the density matrix in terms of an ensemble of pure states and use a variational tensor network formulation for these pure states Transchel et al. (2014). In this case, the density matrix is parametrized according to

[TABLE]

where $\mathopen{|}\psi(\alpha)\mathclose{\rangle}$ is a variational wave function with variational parameters $\alpha$ and $p(\alpha,\bar{\alpha})$ is the associated probability distribution. Crucially, the variational norm associated with the effective Hamiltonian of the master equation $H_{\text{eff}}=H-i/2\sum_{i}c_{i}^{\dagger}c_{i}$ can now be calculated as

[TABLE]

This expression can both be computed efficiently using tensor network methods and corresponds to the natural trace norm when evaluated over the full ensemble. The quantum jump terms of the master equation can be treated in a similar fashion Transchel et al. (2014).

IV.4 Variational quantum Monte-Carlo methods

The central idea behind quantum Monte-Carlo methods is to rewrite a quantum many-body problem in terms of a sampling over a classical probability distribution Batrouni and Scalettar (2011). However, the existence of destructive interference in quantum mechanics can lead to corresponding classical probabilities that are negative, which is the root of the famous sign problem. One common workaround is to sample over the absolute value of the probability distribution instead, but this comes at the price of the complexity of the computation increasing exponentially with the system size Troyer and Wiese (2005). Open quantum many-body systems are especially prone to the sign problem since the eigenvalues of the Liouvillian can even be complex Nagy and Savona (2018). Nevertheless, Monte-Carlo sampling can be useful even in the presence of the sign problem, if the required resources for the Monte-Carlo sampling are lower than for a full solution of the problem.

The first quantum Monte-Carlo simulation of an open quantum many-body problem has been based on a non-variational full-configuration-interaction Monte Carlo algorithm Nagy and Savona (2018), which is better equipped to deal with the sign problem without completely negating it. For the magnetization of a dissipative XYZ model on small lattices, the quantum Monte-Carlo simulation is in excellent agreement with wave-function Monte-Carlo results.

Recently, variational Monte-Carlo methods have been applied to open quantum systems Hartmann and Carleo (2019); Yoshioka and Hamazaki (2019); Nagy and Savona (2019); Vicentini et al. (2019a). These approaches are inspired by using variational wave function corresponding to restricted Boltzmann machines (RBMs) Carleo and Troyer (2017), which were first introduced in the context of neural network simulations. The main idea behind RBM wave functions is shown in Fig. 11, where an additional hidden layer introduces variational parameters associated with the quantum correlations of the many-body state. The entries of the vectorized density matrix are then given by

[TABLE]

where the $W_{ij}$, $a_{i}$, and $h_{j}$ are variational parameters. Interestingly, there is a close connection between RBM wave functions and matrix product states Deng et al. (2017); Chen et al. (2018), however, RBMs are potentially also capable to describe long-range entangled quantum states.

For the variational Monte-Carlo samplings, different norms have been put forward. One possibility is to consider the Hilbert-Schmidt norm of the time evolution Hartmann and Carleo (2019) or the steady state Vicentini et al. (2019a). Interestingly, in the latter case, the variational norm $D$ has been normalized according to the purity $\text{Tr}\{\rho^{2}\}$, i.e.,

[TABLE]

This norm is not biased towards the maximally mixed state as mentioned above. An alternative approach to construct a suitable norm is to mimize the Hermitian $\mathcal{L}^{\dagger}\mathcal{L}$ in close analogy to a ground state problem Yoshioka and Hamazaki (2019). Finally, it is possible to consider the equivalent of an expectation value for vectorized density matrices according to ${}_{\sharp}\langle\rho|\mathcal{L}_{\sharp}|\rho\rangle_{\sharp}/_{\sharp}\langle\rho|\rho\rangle_{\sharp}$ Nagy and Savona (2019). With the respect to the more natural trace norm for density matrices, the RBM approaches behave similarly to the tensor network simulations discussed in Sec. IV.3. However, since RBMs can be applied to two-dimensional models in a straightforward way, it will be very interesting to see how these methods perform for the investigation of dissipative phase transitions, in particular in critical systems.

V Phase space and related methods

Other methods have also been used with relative success in the study of open quantum systems, such as phase space methods, as well as methods based on hierarchy equations. In this section we explain two of such examples, namely, truncated Wigner approximations and Bogoliubov-Born-Green-Kirkwood-Yvon (BBGKY) hierarchies. The resulting methods are very general in purpose, and can be applied to a wide variety of systems, yet in what follow we discuss concrete examples.

V.1 Truncated Wigner approximation

In the context of phase-space related methods, truncated Wigner approximations were first used in Carusotto and Ciuti (2005) to driven-dissipative microcavity polariton system coherently driven into the optical parametric oscillator regime, also reviewed in Carusotto and Ciuti (2013) and revisited in Dagvadorj et al. (2015) as an example of a 2d driven-dissipative non-equilibrium phase transition. The Hamiltonian for the system is given by

[TABLE]

with cavity and photon field operators $\psi_{X,C}(\vec{r},t)$, spatial coordinate $\vec{r}=(x,y)$, $m_{X,C}$ the exciton and photon masses, $g_{X}$ the exciton-exciton interaction strength, and $\Omega_{R}$ the Rabi splitting. One introduces the effect of an external drive (pump) as well as incoherent decay by adding a system-bath Hamiltonian given by

[TABLE]

with $\psi_{l,\vec{k}}(t)$ the Fourier transform of the field operators in real space, $B_{l,\vec{k}}$ and $B_{l,\vec{k}}^{\dagger}$ the bath’s bosonic anihilation and creation operators with energy $\omega_{l,\vec{k}}$, which describes the decay for both excitons and cavity photons. The decay is compensated by an external homogeneous coherent pump $F(\vec{r},t)=f_{p}e^{i(\vec{k}_{p}\cdot\vec{r}-\omega_{p}t)}$, injecting polaritons with momentum $\vec{k}_{p}$ and energy $\omega_{p}$.

By using standard quantum optical methods, one can trace out the bath within the Markovian approximation and obtain a Master equation for the system. There is, however, an alternative approach by means of phase-space techniques. In particular, one can represent the quantum fields as quasiprobability distribution functions. The Fokker-Planck partial differential equation that governs the dynamics of such distributions can be mapped to a stochastic differential equation, which can be solved using different techniques. For the example that we are discussing, one solves the equation on a finite grid with lattice spacing $a$. The most suitable quasiprobability distribution for this example is the Wigner representation, which is also the most suitable one for numerical implementation. By truncating the corresponding Fokker-Planck equation in the limit $(g_{X}/\kappa_{X,C}a^{2})\ll 1$, with $\kappa_{X,C}$ the exciton and photon decay rates, and keeping up to second-order derivatives only, one obtains the stochastic differential equation

[TABLE]

In this equation, $dW_{l=X,C}$ are Wiener noise terms, and $H^{\prime}_{MF}$ is given by

[TABLE]

The resulting stochastic differential equation can then be solved using standard methods and softeare packages for this purpose.

Let us add that recently, this method has also been used to study critical slowing down in photonic lattices Vicentini et al. (2018), as well as extended to disordered quantum many-body system (the so-called optical stochastic unraveling for disordered systems) Vicentini et al. (2019b).

V.2 BBGKY hierarchy equations

It is also possible to study open quantum systems via the so-called Bogoliubov-Born-Kirkwood-Yvon hierarchy Liboff (2003). In a nutshell, this is a hierarchy of equations aimed to describe a system of a large number of interacting particles. As such, the idea is very generic. But as shown in Navez and Schützhold (2010), it can also be applied directly in the context of open dissipative systems in order to obtain a hierarchy of equations for the different reduced density matrices. The way this approach works is quite intuitive. Consider the reduced density matrices for one lattice site $\rho_{\mu}$, for two lattice sites $\rho_{\mu\nu}$, and so on. We separate the correlated parts as $\rho_{\mu\nu}=\rho_{\mu\nu}^{c}+\rho_{\mu}\rho_{\nu}$, as well as $\rho_{\mu\nu\lambda}=\rho_{\mu\nu}^{c}\rho_{\lambda}+\rho_{\mu\lambda}^{c}\rho_{\nu}+\rho_{\nu\lambda}^{c}\rho_{\mu}+\rho_{\mu}\rho_{\nu}\rho_{\lambda}$, and so on. The method that will be discussed in what follows in based on the scaling hierarchy of correlations

[TABLE]

with $|\mathcal{S}|$ the number of lattice sites in set $\mathcal{S}$. The different reduced density matrices can also be computed using the generating functional $\mathcal{F}(\alpha_{\mu})=\log{\rm Tr}(\rho\prod_{\mu}(\mathbb{I}_{\mu}+\alpha_{\mu}))$, with $\alpha_{\mu}$ an arbitrary operator acting on site $\mu$. Using such a functional one has $\rho_{\mu}=\partial\mathcal{F}/\partial\alpha_{\mu}|_{\alpha=0}$, as well as $\rho_{\mu\nu}^{c}=\partial^{2}\mathcal{F}/\partial\alpha_{\mu}\partial\alpha_{\nu}|_{\alpha=0}$, and so on. Next, the Liouville operators $\mathcal{L}_{\mu}$ and $\mathcal{L}_{\mu\nu}$ acting on one and two sites are introduce via the dissipation equation $i\partial_{t}\rho=[H,\rho]+\sum_{\mu}\mathcal{L}_{\mu}\rho+\sum_{\mu\nu}\mathcal{L}_{\mu\nu}\rho/Z$, with $Z$ the coordination number of the Hamiltonian (e.g., the number of tunneling neighbours at any given site for a Hubbard-like Hamiltonian). Following these equations, the time evolution of $\mathcal{F}$ is given by

[TABLE]

Using this equation, one can take derivatives and obtain a set of equations for the correlated density matrices,

[TABLE]

with $\mathcal{L}_{\mu\nu}^{S}=\mathcal{L}_{\mu\nu}+\mathcal{L}_{\nu\mu}$. This hierarchy of equations for the reduced density matrices is preserved in time. Moreover, it allows us to write explicit equations for the one- and two-site density matrices. For the one-site matrix one gets

[TABLE]

and for the two-site matrix one has

[TABLE]

By combining the above expressions with Eq.(48), one can expand in powers of $1/Z$ and obtain different approximations for the one- and two-particle behaviour.

This approach can be implemented for a variety of systems (spins, bosons, fermions…) and has the advantage of being independent of the dimensionality of the system. For instance, in Navez and Schützhold (2010) it was applied to a lattice Bose-Hubbard model. The method can be used to obtain analytical expansions, as well as to facilitate efficient numerical simulations.

VI Linked cluster expansion methods

Methods based on linked-cluster expansions have also been recently put forward in the study of open quantum many-body systems, so far focusing on the study of two-dimensional spin systems with incoherent spin relaxation Biella et al. (2018). The method numerically targets expectation values of observables in the steady state (at long times) of the master equation.

Mathematically, the procedure is as follows: let us assume (without loss of generality) that the Liouvillian can be expanded as a sum of two-body terms, i.e.,

[TABLE]

with $\alpha_{ij}$ some local coupling strength. For the sake of simplicity let us define $k\equiv(i,j)$ as a combined index. The expectation value $O$ of an observable $\hat{O}$ can be expanded in terms of powers of $\alpha_{k}$, i.e.,

[TABLE]

with $n_{k}$ running over all non-negative integers for all $k$. It is clear that all possible polynomials in $\alpha_{k}$ are included in the above expression, which can be reorganized in clusters as follows:

[TABLE]

with $c$ a non-empty set of $k$-indexes identifying the sites belonging to the cluster. The cluster weight $W_{[O]}(c)$ contains all the terms in the expansion with at least one power of $\alpha_{k}$, for all $k$ in $c$, and no powers of $\alpha_{k}$ of $k$ does not belong to $c$. These terms obey the recurrence relation

[TABLE]

with

[TABLE]

being the expectation value of the observable in the steady-state $\rho_{s}(c)$ for the finite cluster $c$. Taking into account symmetries in the system, the expectation value per site in the thermodynamic limit can be written as

[TABLE]

with $L\rightarrow\infty$ the size of the system, the outer sum running over all possible cluster sizes $n$, and the inner sum over all topologically different clusters $c_{n}$ of size $n$, with $l(c_{n})$ their multiplicity. This series expansion can be truncated up to a cluster size $R$, thus giving rise to a plausible approximation method also valid for open systems.

The linked cluster expansion works very well for the dissipative Heisenberg model Biella et al. (2018), where an exact product state solution can be used as a starting point of the expansion. In this case, it is even possible to calculate phase boundaries and critical exponents of a dissipative phase transition between a paramagnet and a ferromagnet. The situation is quite different for the dissipative Ising model, where the expansion series failed to converge even for a 10${}^{\text{th}}$ order expansion Jin et al. (2018).

VII Summary and outlook

The enormous effort to develop novel simulation methods to investigate open quantum many-body systems has enabled us to review a large variety of numerical methods. To be more specific, in this review we considered methods for the Markovian quantum master equation (assuming a weak-coupling limit), including mean-field, stochastic methods, tensor networks, variational methods, quantum Monte-Carlo, truncated Wigner approximation, BBGKY hierarchy equations, and linked cluster expansions. While so far, no method has emerged that is universally optimal for all cases, there have been several very promising developments with different methods for different regimes. Even with such major technical advances discussed in this review, there are still many open problems which are inacessible with these state-of-the-art numerical techniques. To give concrete examples of actual physical problems, one may consider a very common setting in the context of Rydberg atoms where the interaction is often long-ranged and cannot be approximated with just a nearest-neighbour Hamiltonian Schachenmayer et al. (2015); Labuhn et al. (2016); Browaeys et al. (2016). Even TN techniques will face a difficult challenge specially in 2D while encountering such problems although there has been promising developments even in this direction recently O’Rourke and Kin-Lic Chan (2019). Other challenging problems include the existence of AF order in 3D dissipative Ising model which is an open question that appears hard to answer. This is again relevant to ongoing experiments with Rydberg atoms, which one cannot reliably simulate at the moment Carr et al. (2013); Malossi et al. (2014); Helmrich et al. (2018). Phase transitions and universality classes of dissipative models is another class of problem which has proven to be quite difficult for numerical techniques Diehl et al. (2010b); Carmichael (2015); Fink et al. (2017); Biondi et al. (2017).

Certainly, the largest confidence in a simulation result can be achieved if it is reproducible using a complementary simulation approach. Despite these caveats, one can make several key observations about the particular methods covered in this review. The first observation is that mean-field methods are considerably less reliable for open system than their counterparts for closed systems, although the reason for this discrepancy is still an open question. Furthermore, tensor network methods have demonstrated their ability to successfully tackle many hard problems surrounding open many-body systems and resolve long-standing open questions. A particularly interesting and promising case is that of open 2d systems, which is unexplored territory to a great extent. As for the variational methods discussed in this review, there appears to be a tradeoff between the formal suitability of the norm and its efficient computability. It will be interesting to see if and how this tradeoff will be resolved in future work. We provide a summary in Tab. 2 comparing the different techniques we have discussed above.

The progress in recent years in simulating open quantum systems has brought the field to a level where one has a wide range of tools at hand to systematically compare to experimental results, in particular in the context of quantum simulations. Combined with the experimental ease of preparing the steady state of an open quantum system, these are good reasons to believe that the study of strongly-correlated open quantum many-body systems will become a research topic with impact in other areas of science, such as material design and quantum computation.

Acknowledgements.

We thank C. Ciuti for valuable feedback on our manuscript. This work was funded by the Volkswagen Foundation, by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) within SFB 1227 (DQ-mat, project A04), SPP 1929 (GiRyd), and under Germany’s Excellence Strategy – EXC-2123 QuantumFrontiers – 390837967.