operator one may introduce a projector on to the space orthogonal to that spanned by the GGE solutions.

Other applications of perturbation theory in the many-body context have included the thermalisation of open many-body systems [230] and many-body localisation in the presence of an external driving [334].

## 4.4 Mean-field approximations

Mean-field theory can be a powerful approximation scheme to simplify many-body problems. In some cases it is a crude approximation, but can nonetheless help identify possible phases and the structure of the phase diagram (even if not the correct location of phase boundaries). In other cases, mean-field theory can prove to be a good, even exact, approximation. As explored further below, mean-field theory can also be the starting point for a variety of expansions that systematically increase the extent of correlations which are taken into account.

At its core, the key idea of mean-field theory is to replace an interacting many-body system with one where one part of a system evolves in a “mean field” arising due to the effects of other parts of the system. The key approximation here is to neglect correlations between these parts of the system. Mean-field approximations can be formulated in a variety of contexts, including pure quantum dynamics (where it becomes an approximation for the wavefunction), classical problems, and as we discuss here, for open quantum systems. In this last case, the essential structure of a mean-field theory is to assume some form of factorisable density matrix:

$$ \hat{\rho} = \bigotimes_{\mu} \hat{\rho}_{\mu}, \qquad (94) $$

where $\hat{\rho}_{\mu}$ is the density matrix of the $\mu$-th part of the system—for example, factorisation in real space (for which part $\mu$ means different sites), or momentum space (for which part $\mu$ means different momentum states). We discuss many different examples of this in the rest of this section. We note that this formal statement of the Ansatz will hold for most, but not all, of the cases we discuss below; the conceptual structure of a factorised Ansatz will however hold in all cases. Such a description is clearly abstract, in that there can be many ways to divide up the system. As such, there can often be more than one possible mean-field decoupling, as we discuss in the sections below.

In many of the cases we discuss below, there can be alternative statements of how one proceeds to derive a mean-field theory—for example, replacing operators by classical variables, factorising expectations of products of operators by products of expectations, or assuming only pairwise correlations among particles exist. While these rules may appear simpler than the statement of factorisation, they obscure the nature of the approximation that is actually being made, and can give the impression that different types of mean-field theory are distinct and unrelated. Our aim in this section is to make clear the links between different forms of mean-field theory (as well as their relation to beyond-mean-field approaches). These links come precisely from considering that mean-field theories are all related to assumptions of factorisation. We will nonetheless make clear in each of the following sections where an alternate approach can lead to the same end point.

When considering how to use mean-field approaches, two significant challenges arise. The first challenge is to decide which form of mean-field factorisation is most appropriate. We will discuss below some general considerations on when a given mean-field approach is likely to be accurate, which can be used to guide such a decision. Beyond this, one can consider known experimental behaviour of a system, or similarity to other known systems to propose a particular form of factorisation. The second challenge is how to then use the Ansatz to find either the steady state or the dynamics. For an isolated system in its ground state, one can use a variational approach to minimise the ground state energy; i.e. for some variational
Ansatz $|\Psi_v\rangle$ one may vary parameters to minimise the energy $E_v = \langle\Psi_v|\hat{H}|\Psi_v\rangle$. For a system in thermal equilibrium at non-zero temperature, one may similarly consider a variational Ansatz for the density matrix $\hat{\rho}_v$ and minimise the free energy, $F = \text{Tr}[\hat{H}\hat{\rho}_v + k_B T \hat{\rho}_v \ln \hat{\rho}_v]$. For an open quantum system the equivalent procedure is less obvious since the steady state is not determined by optimising a thermodynamic potential. Ultimately one seeks a solution of the equation $\hat{\mathcal{L}}[\rho] = 0$, however any restricted Ansatz may be unable to exactly solve this equation. One can convert this condition into a minimisation problem in various ways; a popular choice involves taking the modulus squared of the vectorised form, i.e. $\|\hat{\mathcal{L}}[\hat{\rho}]\|^2 = \langle\rho|\hat{\mathcal{L}}^\dagger\hat{\mathcal{L}}|\rho\rangle$.

In many cases, one is interested in the time evolution rather than just in the steady state. Here one can use an alternative variational approach, the time-dependent variational principle [335–337] (TDVP). This consists of assuming that one has a state $\hat{\rho}_t$ within the variational manifold at time $t$, and then finding the state at time $t + \delta t$ by minimising some distance measure $D[\hat{\rho}_{t+\delta t}, \hat{\rho}_t + \delta t \hat{\mathcal{L}}[\hat{\rho}_t]]$, denoting the distance between the best variational state at the new time and the time evolution of the previous state. In the analogous approach for pure quantum states it is clear that the distance measure to be used is the Euclidean norm $D[|\psi_1\rangle, |\psi_2\rangle] = \sqrt{\langle\psi_1 - \psi_2|\psi_1 - \psi_2\rangle}$; this will maximise overlap between the variational state at the next time step and its “correct” value. As with finding the steady state, the choice of correct distance function to minimise for density matrices is less clear, but approaches minimising the Euclidean norm [190] are often used.

### 4.4.1 Real space factorisation and Gutzwiller mean-field theories

In bosonic lattice models (such as the Bose-Hubbard or related models) or quantum spin chains, one natural mean-field factorisation is between different lattice sites. In this case $\hat{\rho} = \bigotimes_i \hat{\rho}_i$ holds with $i$ running over the site labels. Such an Ansatz is closely related to the Gutzwiller Ansatz first introduced in the context of the Fermi-Hubbard model [338]. However, as we will discuss later, for bosonic models this factorised Ansatz becomes exact in the limit of infinite lattice connectivity in contrast to fermionic analogues. As such, we will refer to mean-field theories of this form as Gutzwiller mean-field theories [339].

An alternative context where similar factorisation holds comes for problems involving many emitters coupled to a single cavity [231]. Here, a common factorisation is between the state of the cavity modes and the state of the emitters, and so in this case Eq. (94) again holds with $i$ running over the cavity and the emitters.

To discuss this spatial factorisation more explicitly, suppose one considers a Lindblad master equation where the Hamiltonian takes the form $\hat{H} = \sum_i \hat{H}_i + \sum_{i \neq j} \hat{H}_{ij}$, describing terms that act on a single part of the system, $i$ and terms that couple two parts, $i$ and $j$. If we further assume that each jump operator involves operators on only one part of the system, $\hat{L}_\mu = \hat{L}_{\mu,i}$, then one may find a simplified equation of motion within the space of factorised density matrices:

$$ \frac{d\hat{\rho}_i}{dt} = -i[\hat{H}_i + \hat{\Theta}_i, \hat{\rho}_i] + \sum_{\mu} \hat{D}[\hat{L}_{\mu,i}, \hat{\rho}_i], \qquad (95) $$

where $\hat{\Theta}_i = \sum_{j \neq i} \text{Tr}_j(\hat{\rho}_j \hat{H}_{ij})$ describes the mean field of the other parts of the system, i.e. parts $j$, acting back on part $i$. In cases where the jump operators are not local, one can derive a more complex equation involving partial traces of these terms. The structure of this equation means that the evolution of different $\hat{\rho}_i$ are coupled, but only through a mean-field expectation. This coupling means the evolution of each $\hat{\rho}_i$ depends on expectations of terms on other sites.

To give concrete examples of such an approach, we may consider the XYZ spin model and dissipative Bose-Hubbard model that were introduced in Sec. 3.1.5. For the XYZ model, the Hamiltonian in Eq. (28) contained only pairwise interaction terms involving nearest-neighbour
sites, while the dissipation, Eq. (29) acted on a single site. The mean-field decoupling then leads to

$$ \hat{\Theta}_i = \sum_{j \in \partial_i} \mathrm{Tr} \left[ \hat{\rho}_j \sum_{\alpha} J_{\alpha} \hat{\sigma}_i^{\alpha} \hat{\sigma}_j^{\alpha} \right], $$

where $j \in \partial_i$ indicates a sum over the sites $j$ that are nearest neighbours to site $i$. We can thus write the resulting mean-field equation of motion as:

$$ \frac{d\hat{\rho}_i}{dt} = -i \left[ \sum_{\alpha} h_i^{\alpha} \hat{\sigma}_i^{\alpha}, \hat{\rho}_i \right] + \hat{\mathcal{D}}[\sqrt{\gamma} \hat{\sigma}_i^{-}, \hat{\rho}_i], \quad h_i^{\alpha} = J_{\alpha} \sum_{j \in \partial_i} \mathrm{Tr} \left[ \hat{\rho}_j \hat{\sigma}_j^{\alpha} \right]. \quad (96) $$

In the limit of weak coupling $J^\alpha$, the steady state will be dominated by the dissipation, leading to a spin down state, with $h_i^x = h_i^y = 0$. When the couplings $J^x, J^y$ become large enough, one has a transition to a state where there is a self-consistent solution with non-zero $h_i^x, h_i^y$.

Applying the same logic to the Bose-Hubbard Hamiltonian given by Eq. (36) we find

$$ \frac{d\hat{\rho}_i}{dt} = -i \left[ \phi_i \hat{a}_i^\dagger + \phi_i^* \hat{a}_i + \frac{U}{2} \hat{n}_i (\hat{n}_i + 1), \hat{\rho}_i \right] + \sum_{\mu} \hat{\mathcal{D}}[\hat{L}_{i,\mu}, \hat{\rho}_i], \quad \phi_i = -t_H \sum_{j \in \partial_i} \mathrm{Tr} [\hat{\rho}_j \hat{a}_j]. \quad (97) $$

Here we have allowed for a general set of on-site dissipative processes, as described when introducing this model. For this mean-field factorisation, there is a transition between a superfluid state with non-zero $\phi_i$ and a normal state where $\phi_i = 0$. In the absence of dissipation the ground state with $\phi_i = 0$ is a Mott insulating state with integer occupation per site, however this is modified by any pumping or dissipation (or non-zero temperature).

We may note that this form of real-space mean-field factorisation is not appropriate for the Fermi-Hubbard model. This is because physical Hamiltonians always conserve fermion parity, whereas decoupling the hopping term for the Fermi-Hubbard model would give a model with single fermion operators, and would require a state with a non-zero expectation of a single fermion operator. We discuss an alternate approach that resolves this issue—dynamical mean-field theory—in section 4.5.

When considering Gutzwiller mean-field theories, it is worth distinguishing the mean-field approximation from two other additional approximations that are sometimes associated with mean-field factorisations in real space, but which are in fact additional approximations.

**Spatial homogeneity.** The first is the assumption of spatial homogeneity. If one considers Eq. (95) for lattice model, one may note that in general the equation allows for $\hat{\rho}_i$ to evolve independently on each site. If one is focused on steady states, and one expects translational invariance, one can make an additional assumption that $\hat{\rho}_j = \hat{\rho}_i$ for all $i, j$, reducing to a single self-consistent equation of motion. One may though note that this is an additional assumption on top of mean-field theory. Note further that even when the model is translationally invariant, the steady state need not be. For example, antiferromagnetic coupling between sites might lead to a state with staggered order, and more complex forms of coupling can lead to other commensurate or incommensurate states, e.g. [234, 340–342].

**Semiclassical approximation.** The second approximation is the semiclassical approximation, when considering bosonic modes. This is the assumption that bosonic operators $\hat{a}$ can be replaced by classical fields $\psi$ throughout the equations of motion. There are some cases where mean-field theory in real space does imply a semiclassical approximation. One such example is the Dicke model, defined by Eqs. (26) and (27). If the mean-field factorisation introduced above is made between the cavity, $\hat{\rho}_c$, and the emitters, $\hat{\rho}_{e_i}$, i.e. $\hat{\rho} = (\otimes_i \hat{\rho}_{e_i}) \otimes \hat{\rho}_c$, then the
mean-field equation for the cavity mode becomes:

$$ \frac{d\hat{\rho}_c}{dt} = -i[\omega_c \hat{a}^\dagger \hat{a} + \phi(\hat{a}^\dagger + \hat{a}), \hat{\rho}_c] + \kappa \hat{D}[\hat{a}], \quad (98) $$

where the term $\phi$ comes from a trace over the emitters, i.e. $\phi = g \sum_i \text{Tr}[\hat{\sigma}_i^x \hat{\rho}_{e_i}]$ for the model in Eq. (26). One may easily see that the stationary solution of this equation is a coherent state, $\hat{\rho}_c = |z\rangle\langle z|$ with $|z\rangle = \exp(z\hat{a}^\dagger - z^*\hat{a})|0\rangle$, $z = -\phi/(\omega_c - i\kappa/2)$. Thus, in this specific case the mean-field factorisation implies that the state of the cavity bosonic mode is one for which the semiclassical approximation holds.

Despite examples such as that above, Gutzwiller mean-field factorisation is not in general equivalent to assuming a coherent state or that replacing operators with fields holds. For example, in the Bose-Hubbard model discussed above, the on-site Hamiltonian involves interacting terms, and as a result, the steady state is not a coherent state.

### 4.4.2 Momentum space factorisation

A second form of mean-field factorisation that often arises is that in momentum space. This corresponds to considering a density matrix of the form $\hat{\rho} = \otimes_k \hat{\rho}_k$, with each term involving only the pair of momenta $k, -k$. This structure is chosen to be compatible with states with overall zero momentum while allowing only correlations between excitations with equal and opposite momenta. Such approaches have been used to consider the dynamics of pairing in the BCS model (see e.g. [343, 344]), as well as for driven dissipative problems [244, 245] or for lossy gases [345].

Momentum factorisation is often accompanied by other special properties, that are sometimes thought of as intrinsic properties of mean-field approaches. As an example, we will next show that the standard BCS ground state Ansatz can be recovered from two different starting points: one is from considering factorisation in momentum space, another is from considering translationally invariant Gaussian states. From the perspective introduced at the start of this section, the key feature of the BCS state is its factorised nature, rather than its Gaussian form, although (as we will show), both assumptions lead to the same end result. We consider the BCS ground state because this is the most famous example of such momentum state factorisation, and a particularly simple case to discuss. These considerations can be extended to open quantum systems with some additional work, but for simplicity and clarity we restrict this discussion to a closed-system example.

One may write down the momentum-factorised Ansatz for the ground-state wavefunction of the BCS model by three considerations: assuming factorisation of the state in momentum space, $|\Psi\rangle = \prod_k |\psi_k\rangle$, requiring a state with definite fixed fermion number parity (i.e., not mixing odd and even fermion number states), and requiring a state with zero total momentum. These requirements yield the state:

$$ |\Psi\rangle = \prod_{k \ge 0} \left( \psi_k^{(0)} + \sum_{\sigma, \sigma' \in \uparrow, \downarrow} \psi_{k, \sigma \sigma'}^{(1)} \hat{c}_{-k, \sigma}^{\dagger} \hat{c}_{k, \sigma'}^{\dagger} + \psi_{k}^{(2)} \hat{c}_{-k, \downarrow}^{\dagger} \hat{c}_{-k, \uparrow}^{\dagger} \hat{c}_{k, \uparrow}^{\dagger} \hat{c}_{k, \downarrow}^{\dagger} \right) |0\rangle. \quad (99) $$

Note that the notation $k \ge 0$ in this product indicates a product over the half-space of momentum $k$, since each term contains all contributions from both $\pm k$. The corresponding density matrix $\hat{\rho} = |\Psi\rangle\langle\Psi|$ constructed from this takes exactly the form of momentum-factorised Ansatz discussed above.

If one further considers using this Ansatz for the ground state of the BCS Hamiltonian (39) one finds that one can consider a restricted version:

$$ |\Psi\rangle = \prod_k \left( \cos \theta_k + \sin \theta_k e^{i\phi_k} \hat{c}_{-k, \downarrow}^{\dagger} \hat{c}_{k, \uparrow}^{\dagger} \right) |0\rangle. \quad (100) $$
Note that in this expression the product over $k$ includes all terms. The simplification to this form can be understood by noting that after factorising over momentum, the non-vanishing terms in Eq. (39) are all quadratic in fermions, and the anomalous terms take the forms $\hat{c}_{-k,\downarrow}^{\dagger}\hat{c}_{k,\uparrow}^{\dagger}$ or $\hat{c}_{-k,\downarrow}\hat{c}_{k,\uparrow}$, leading to the form given. The use of $\theta_k, \phi_k$ to write the expression ensures normalisation. This quadratic form also means, as discussed next, that the BCS mean-field theory can alternatively be found from a Gaussian Ansatz.

**Gaussian Ansatz.** As in our discussion of real-space factorisation, there are features that arise in the momentum space factorisation of some problems (including the BCS Ansatz) that are not generic properties of mean-field theory, but which frequently accompany with mean-field theory. In particular, one may note that the BCS Ansatz written above can equivalently be recovered by assuming a Gaussian state. Such a state can be written either in real space or momentum space. In real space this means writing:

$$|\Psi\rangle = \exp\left(\iint d\mathbf{r}d\mathbf{r}'\phi(\mathbf{r}-\mathbf{r}')\hat{c}_{\uparrow}^{\dagger}(\mathbf{r})\hat{c}_{\downarrow}^{\dagger}(\mathbf{r}')\right)|0\rangle. \qquad (101)$$

In writing this we have imposed translational invariance of this Gaussian Ansatz. As a result, if we switch to momentum space operators by using

$$\hat{c}_{\sigma}^{\dagger}(\mathbf{r}) = \frac{1}{L^{d/2}} \sum_{\mathbf{k}} \hat{c}_{\mathbf{k},\sigma}^{\dagger} e^{-i\mathbf{k}\cdot\mathbf{r}}, \qquad (102)$$

one finds the equivalent momentum-space representation:

$$|\Psi\rangle = \exp\left(\sum_{\mathbf{k},\mathbf{k}'} \tilde{\phi}_{\mathbf{k},\mathbf{k}'} \hat{c}_{\mathbf{k},\uparrow}^{\dagger} \hat{c}_{\mathbf{k}',\downarrow}^{\dagger}\right)|0\rangle, \qquad (103)$$

with

$$\tilde{\phi}_{\mathbf{k},\mathbf{k}'} = \frac{1}{L^d} \iint d\mathbf{r} d\mathbf{r}' \phi(\mathbf{r}-\mathbf{r}') e^{i(\mathbf{k}\cdot\mathbf{r}+\mathbf{k}'\cdot\mathbf{r}')} = \chi_{\mathbf{k}} \delta_{\mathbf{k},-\mathbf{k}'}, \quad \chi_{\mathbf{k}} \equiv \int d\mathbf{s} \phi(\mathbf{s}) e^{i\mathbf{k}\cdot\mathbf{s}}, \qquad (104)$$

which then yields the state:

$$|\Psi\rangle = \exp\left(\sum_{\mathbf{k}} \chi_{-\mathbf{k}} \hat{c}_{-\mathbf{k},\uparrow}^{\dagger} \hat{c}_{\mathbf{k},\downarrow}^{\dagger}\right)|0\rangle = \prod_{\mathbf{k}} \left(1 + \chi_{-\mathbf{k}} \hat{c}_{-\mathbf{k},\uparrow}^{\dagger} \hat{c}_{\mathbf{k},\downarrow}^{\dagger}\right)|0\rangle. \qquad (105)$$

This recovers the BCS Ansatz (up to normalisation) if we identify $\chi_k = \tan(\theta_k)e^{i\phi_k}$. Note that the combination of translational invariance and Gaussianity implies the state must factorise in momentum space. In contrast the state in Eq. (101) does not factorise in real space—i.e. it has correlations between the state at different positions in space.

As with our discussion of the relation between real-space factorisation and the semiclassical approximation, not all momentum space factorisations necessarily imply a Gaussian state. Indeed, the most general momentum space factorisation given in Eq. (99) is clearly not Gaussian. Further examples of non-Gaussian but momentum-factorised states are given by the Jastrow wave function for fermionic or bosonic Mott insulators [346, 347], which can be factorised in momentum space, when the system is translational invariant, yet does not satisfy Wick's theorem because of the Jastrow projector.
### 4.4.3 Particle-based factorisation: Hartree-Fock theory and Gross-Pitaevskii equation

The above two approaches consider factorisation between different single-particle orbitals, i.e. different modes of the system. One can also consider factorisation between the states of different particles. If particles are distinguishable, this would have the straightforward meaning that, as a wavefunction, $\Psi(\mathbf{r}_1, \mathbf{r}_2, \dots, \mathbf{r}_N) = \prod_{i=1}^N \psi_i(\mathbf{r}_i)$. For indistinguishable particles, this simple product must be replaced by a Slater determinant of a matrix of wavefunctions for fermions, or a matrix permanent for bosons. These constructs are required so that the wavefunction obeys the required anti-symmetry or symmetry. Such an approach corresponds to Hartree-Fock theory [348]. Similar expressions can be written for density matrices as we discuss further below.

For bosons (but not for fermions) one can also consider the special case in which all particles are in the same state. In this case the matrix permanent takes a simple form: $\Psi(\mathbf{r}_1, \mathbf{r}_2, \dots, \mathbf{r}_N) = \prod_{i=1}^N \psi(\mathbf{r}_i)$. Such an approximation is the basis for the widely used Gross-Pitaevskii equation that is the mean-field theory for the weakly interacting dilute Bose gas (WIDBG) model, introduced in Sec. 3.1.5

In second-quantised notation corresponding to Eq. (37), the product state Ansatz takes the form:

$$|\Psi(t)\rangle = \frac{1}{\sqrt{N!}} \left( \sum_{\mathbf{k}} \tilde{\psi}_{\mathbf{k}}(t) \hat{a}_{\mathbf{k}}^{\dagger} \right)^N |0\rangle, \quad (106)$$

creating $N$ particles in the single-particle state defined by the Fourier transform, $\tilde{\psi}_{\mathbf{k}}$ of the wavefunction $\psi(r)$. The equation of motion corresponding to evolution of the function $\psi(r)$ under purely Hamiltonian evolution is then (in the large $N$ limit) the Gross-Pitaevskii equation:

$$i \frac{d}{dt} \psi(\mathbf{r}, t) = \left( \frac{-\nabla^2}{2m} + N U |\psi(\mathbf{r}, t)|^2 \right) \psi(\mathbf{r}, t). \quad (107)$$

It is worth noting that the same Gross-Pitaevskii equation can also be derived from an alternate picture, based on coherent states. Here one considers

$$|\Psi(t)\rangle = \exp\left(\sqrt{N}\sum_{\mathbf{k}} \tilde{\psi}_{\mathbf{k}}(t) \hat{a}_{\mathbf{k}}^{\dagger}\right) |0\rangle = \exp\left(\sqrt{N}\int d^d\mathbf{r} \psi(\mathbf{r}, t) \hat{a}(\mathbf{r})\right) |0\rangle. \quad (108)$$

As such, this limit also corresponds to the semiclassical approximation for the WIDBG model.

One may extend this Ansatz to consider a density matrix $\hat{\rho} = |\Psi\rangle\langle\Psi|$, which allows one to consider evolution under the Lindblad master equation written in Eq. (38). Using the Ansatz $\hat{\rho} = |\Psi\rangle\langle\Psi|$ with $|\Psi\rangle$ as above one then finds the complex Gross-Pitaevskii equation:

$$i \frac{d}{dt} \psi(\mathbf{r}, t) = \left( \frac{-\nabla^2}{2m} + N U |\psi(\mathbf{r}, t)|^2 + \frac{i}{2} [\kappa_+ - \kappa_- - N \kappa_-^{(2)} |\psi(\mathbf{r}, t)|^2] \right) \psi(\mathbf{r}, t). \quad (109)$$

Such equations have been widely used to study the physics of nonequilibrium polariton condensates [17, 163]. They also can be derived as the saddle point of a Schwinger-Keldysh path integral, as we will discuss in Sec. 4.7

## 4.5 Large connectivity limit and dynamical mean field theory

As noted above, mean-field theory is an approximation. Since mean-field approximations involve neglecting correlations between different parts of the system, such approximations should best hold when such correlations are weak. A key property controlling this is the connectivity of sites, $z$, i.e. how many other sites affect a given site. When the connectivity is large, $z \gg 1$, the field seen by a given site comes from summing the contributions of many
sites. That is, for large connectivity, a central limit theorem may apply [349], so that the sum over contributions can be replaced by its mean value. In the case of quantum systems one may show that the large connectivity limit for spins and bosons recovers mean-field theory in real space, as discussed above. For other systems—e.g. fermions or disordered systems—the large connectivity limit is more subtle, for example, for fermions the relevant limit is Dynamical Mean-Field Theory (DMFT) [350–352]. Below we first discuss different ways of reaching the large connectivity limit where mean-field theory becomes valid, and then discuss DMFT for both bosons and fermions.

### 4.5.1 Large connectivity and validity of mean-field theories

As discussed above, when the number of neighbours *z* of a given lattice site is large, statistical and quantum fluctuations induced by the neighbouring sites become small and can be treated in an approximate way, while the local, on-site physics must always be accounted for exactly. From this perspective, mean-field theories can be formulated by studying the large connectivity limit *z* → ∞. As we discuss below, this large connectivity limit can be reached in different ways in different models.

**Bosonic and spin lattice models.** Let us consider a model on a regular lattice with connectivity *z*, such as the Bose-Hubbard model or its hard-core limit given by quantum spin chains. In such lattice models, high connectivity generally requires high dimensionality (e.g. for one-dimensional lattices, *z* = 2, while for cubic three-dimensional lattices *z* = 8). As we will see in the next section in more detail, the limit *z* → ∞ reduces to Gutzwiller mean-field theory, i.e. to a factorised Ansatz for the system density matrix [240, 339]. This is ultimately due to the fact that bosonic and spin operators which can take a non-zero expectation value (describing spin or bosonic coherent states that break a symmetry). Corrections to mean-field theory can be included by treating the leading 1/*z* corrections within Dynamical Mean-Field Theory [353, 354], as we will discuss in the next section for driven dissipative bosons.

**Fermionic lattice models.** For fermionic models on lattices with connectivity *z*, the limit *z* → ∞ does not reduce to a static mean-field theory. This is due to the fact that fermionic parity cannot be broken, i.e. a single fermionic mode cannot condense [350], as we mentioned already in Sec. 4.4.1. Rather in the *z* → ∞ limit a fermionic lattice model is mapped onto a self-consistent quantum impurity model describing a single site of the lattice embedded in a fermionic bath describing self-consistently the rest of the system. The problem therefore retains local quantum fluctuations even in the *z* → ∞ limit. DMFT has been applied to the nonequilibrium dynamics of correlated electrons [351, 355]. The effect of Markovian dissipation has been studied for example in Ref. [356].

**Fully connected spin models and models with one-to-many coupling.** The mean-field approximation can be also formulated as an exact solution of a model on a fully connected lattice of *N* sites, where each lattice site is coupled to *N*−1 neighbours, in the thermodynamic limit *N* → ∞. In these fully connected models the thermodynamic limit therefore coincides with the limit in which mean-field becomes exact, as opposed to lattice models discussed above in which fluctuations beyond mean-field survive in the thermodynamic limit.

Fully connected models can arise either in presence of all to all couplings or, effectively, in cases of models with one-to-many coupling such as the Dicke model, in which a single photon mode acts as a central site, coupling to many emitters as satellite sites. In such models, there are cases where a static mean-field theory can be shown to be exact in the
limit that the number of emitters goes to infinity [357–360]. This holds as the model has a high connectivity (set by the number of emitters), and so the arguments noted above should apply. One may note however that counterexamples exist, particularly in cases where the central site has a weak effect back on the other sites. In such cases, the mean-field effect of the satellite sites vanishes, and so the otherwise-subleading fluctuation terms are dominant [361, 362].

**Fully connected models with quenched disorder.** Here again the fully connected nature of the model makes mean-field theory exact in the thermodynamic limit. However the presence of quenched disorder makes mean-field theory non-trivial as for the fermionic case. Specifically, the large connectivity limit does not coincide with a factorised Ansatz, as mentioned above, but is again described by a dynamical mean field [352]. In the context of open quantum systems examples have been recently studied of random Lindbladian with all to all couplings, for which DMFT provides the exact solution in the thermodynamic limit [363, 364].

### 4.5.2 DMFT for driven-dissipative bosonic lattices

To make concrete the ideas in this section we discuss DMFT for driven-dissipative bosons on a lattice with coordination number z and follow the discussion presented in Ref. 284. In the following, since we are interested in the large connectivity limit we define the hopping as $t_H \equiv \tilde{t}_H/z$, and $\tilde{t}_H$ will be finite as $z \to \infty$, which is the correct bosonic scaling to have a well defined infinite connectivity limit [240, 353, 354]. For concreteness we consider the driven-dissipative Bose-Hubbard model introduced in Eq. (36). The Hamiltonian can be written by separating explicitly terms which are local (on-site) from those that couple neighbouring sites

$$ \hat{H} = -\frac{\tilde{t}_H}{z} \sum_{\langle i,j \rangle} (\hat{a}_i^\dagger \hat{a}_j + \text{H.c.}) + \sum_i \hat{h}_{\text{loc}}[\hat{a}_i^\dagger, \hat{a}_i], \quad (110) $$

where the first term describes the hopping term while the second one accounts for generic local interactions. We will also consider a generic set of local jump operators, $\hat{L}_i$, defined on each site i.

The starting point to formulate DMFT for driven-dissipative bosons is the Keldysh action formulation of the Lindblad lattice problem, which following Sec. 3.2 we can write as

$$ S = \int dt \sum_{\zeta=f,b} s_{\zeta} \left[ \sum_{i} \psi_{i\zeta}^{*} i \frac{d}{dt} \psi_{i\zeta} - \sum_{i} h_{\text{loc}}(\psi_{i\zeta}^{*}, \psi_{i\zeta}) \right] \\ -i \int dt \sum_{i} \left[ L_{ib} L_{if}^{*} - \frac{1}{2} (L_{if}^{*} L_{if} + L_{ib}^{*} L_{ib}) \right] \\ + \frac{\tilde{t}_{H}}{z} \int dt \sum_{\zeta=f,b} s_{\zeta} \sum_{\langle i,j \rangle} (\psi_{i\zeta}^{*} \psi_{j\zeta} + \text{H.c.}), \quad (111) $$

where $s_{f,b} = \pm 1$ and $L_{i\zeta}, L_{i\zeta}^*$ are the jump operators written in terms of the bosonic fields $\psi_{i\zeta}^*, \psi_{i\zeta}$.

By taking the large connectivity limit one can formally map the effective action of the problem, obtained by integrating out all but one site of the lattice, onto a quantum impurity model. This describes an interacting Markovian single site, characterised by the same local Hamiltonian $H_i$ and local jump operators $L_{i\zeta}$ entering Eq. (95), coupled to a time-dependent field $\Theta_{i,\text{eff}}(t)$ acting as a coherent drive and a non-Markovian quantum bath (Fig. 14, top panel). These take into account the effect of the neighbouring sites and have to be determined
self-consistently through the calculation of impurity properties. The DMFT effective Keldysh action reads (assuming for simplicity a normal phase where the bosons remain incoherent):

$$ S_{\text{eff}}[\psi_{\zeta}^{*}, \psi_{\zeta}] = S_{\text{loc}}[\psi_{\zeta}^{*}, \psi_{\zeta}] + \int dt \sum_{\zeta=f,b} s_{\zeta} \Phi_{\text{eff}\zeta}^{*}(t) \psi_{\zeta}(t) \\ - \frac{1}{2} \iint dt dt' \sum_{\zeta, \zeta'=f,b} s_{\zeta} s_{\zeta'} \psi_{\zeta}^{*}(t) \Lambda^{\zeta\zeta'}(t, t') \psi_{\zeta'}(t'). \quad (112) $$

The first term in Eq. (112), $S_{\text{loc}}[\psi_{\zeta}^{*}, \psi_{\zeta}]$ is the local, on-site, contribution of the original lattice problem which can be read off from the first two terms in Eq. (111) and therefore includes interactions, as well as Markovian incoherent drive and dissipation leading to off-diagonal terms in Keldysh space. The second and third terms describe the feedback of the rest of lattice onto the site through its neighbours, in terms of an effective coherent drive $\Phi_{\text{eff}\zeta}^{*}(t)$ and an effective non-Markovian bath with hybridisation function $\Lambda^{\zeta\zeta'}(t, t')$. Both these quantities have to be determined self-consistently, in particular the effective coherent drive reads:

$$ \Phi_{\text{eff}\zeta}^{*}(t) = \tilde{t}_{H} \Phi_{\zeta}^{*}(t) + \int dt' \sum_{\zeta'=f,b} s_{\zeta}^{\prime} \Phi_{\zeta'}^{*}(t') \Lambda^{\zeta\zeta'}(t', t), \quad (113) $$

and has two contributions, the first coming from the average of the bosonic field as in Gutzwiller mean field theory $\Phi_{\zeta}^{*} = \langle \psi_{\zeta}^{*} \rangle_{\text{eff}}$ and the second one coming from the non-Markovian bath, a non-trivial finite z correction accounting for the feedback of neighbouring sites on the local effective field [354,365,366]. The notation $\langle \dots \rangle_{\text{eff}}$ appearing above indicates an integral weighted by the effective action, i.e.

$$ \langle F(\psi_f, \psi_b) \rangle_{S_{\text{eff}}} = \int \mathcal{D}[\psi_f, \psi_b] F(\psi_f, \psi_b) \exp(iS_{\text{eff}}). \quad (114) $$

The DMFT self-consistency condition depends on the specific choice of the lattice and provides in general a functional relation between the hybridisation function of the non-Markovian bath to the impurity connected Green's function

$$ G^{\zeta\zeta'}(t, t') = -i \langle \psi_{\zeta}(t) \psi_{\zeta'}^{*}(t') \rangle_{S_{\text{eff}}}. \quad (115) $$

For example, on a Bethe lattice with connectivity z the DMFT self-consistency condition takes the form [366]

$$ \Lambda^{\zeta\zeta'}(t, t') = \frac{\tilde{t}_{H}^2}{z} G^{\zeta\zeta'}(t, t'; \Lambda). \quad (116) $$

Note here we have indicated that G depends self-consistently on $\Lambda$ via the form of $S_{\text{eff}}$. The DMFT solution of the original Markovian lattice problem thus requires one to solve the Keldysh action (112), computing in particular the impurity Green's function (115) and the average of the bosonic field (113), for given values of the non-Markovian bath $\Lambda$ and effective field $\Phi$, and to iterate (113-116) until self-consistency.

It is instructive at this point to take explicitly the limit of infinite coordination number $z \to \infty$. In this limit, the DMFT effective action (112) becomes completely local in time

$$ S_{\text{eff}}[\psi_{\zeta}^{*}, \psi_{\zeta}] \xrightarrow{z=\infty} S_{\text{loc}}[\psi_{\zeta}^{*}, \psi_{\zeta}] + \int dt \sum_{\zeta} s_{\zeta} \Phi_{\text{eff}\zeta}^{*}(t) \psi_{\zeta}(t), \quad (117) $$

since the contribution from the non-Markovian bath scales as $1/z$ (see Eq. (116)), and as such can be unfolded back into a master equation for a single-site density matrix $\rho$, which
<figure><img src="image_17.png" /><figcaption>Figure 14: Dynamical Mean-Field Theory (DMFT) for bosonic open quantum systems. In the limit of infinite lattice connectivity, a driven-dissipative bosonic array is mapped via Gutzwiller mean-field theory onto a single-site in a self-consistent field. DMFT includes $1/z$ corrections from the neighbouring sites via the solution of a dissipative quantum impurity problem, where the single site is coupled to a non-Markovian (i.e. frequency dependent) bath, whose spectrum and distribution function, encoded in the retarded and Keldysh components of the bath hybridisation function $\Lambda^R(t, t')$, $\Lambda^K(t, t')$, are determined self-consistently by the local Green's function. Adapted from Ref. [284].</figcaption></figure>

satisfies Eq. (95) where $\mathcal{L}$ is the local part of the Lindbladian and the feedback from the neighbouring sites is carried by $\Phi(t) = \text{Tr}[\hat{a}\hat{\rho}(t)]$. This corresponds to a factorised Gutzwiller-like Ansatz for the lattice many body density matrix. In other words, as for equilibrium or closed systems [354, 366] also for driven-dissipative lattice systems the infinite connectivity limit of bosons coincides with Gutzwiller mean-field theory. We note that when $\Phi = 0$, this mean-field describes completely uncoupled sites, while DMFT ($z < \infty$) captures the feedback from neighbouring sites through the self-consistent bath $\Lambda$.

To gain insights on the physics behind the DMFT approximation it is useful to compare it with the Gutzwiller mean-field theory discussed in Sec. 4.4.1. This approach corresponds to a factorized density matrix Ansatz over different lattice sites and it is therefore equivalent to the solution of a single site problem in a self-consistent field (Fig. 14, bottom panel). By construction such a mean-field theory cannot describe short-range correlations and processes arising from coupling to neighboring sites, either of coherent or dissipative origin. This is particularly true for interacting normal phases such as Mott insulators of bosons or quantum Zeno phases which are featureless within Gutzwiller mean-field theory. The corrections due to DMFT in the large connectivity limit reintroduce part of these short-range fluctuations via the coupling to an effective self-consistent bath as in Eq. (112), and lead to a richer descriptions of isolated or dissipative Mott insulating phases [367, 368]. The price to pay is the solution of a dissipative quantum impurity model, and in particular the calculation of its Green's functions, which although simplified with respect to the full master equation, still poses a major computational challenges, particularly in presence of non-linear dissipative processes. Different approaches have been developed to tackle these types of dissipative impurities, including exact diagonalisation of the impurity Lindbladian [356, 369, 370] or strong coupling self-consistent hybridisation expansions [284, 371, 372].
We may note that there have been other approaches to bosonic open quantum systems— such as the self-consistent projection operator approach [373]—which have close connections to DMFT.

### 4.5.3 DMFT for driven-dissipative fermions

We conclude our discussion of DMFT by briefly discussing its application to driven-dissipative fermionic problems, and the main differences from the bosonic version.

A first key point is that the large connectivity limit of fermions involves a different rescaling from the bosonic problem discussed above. For fermions one should instead rescale the hopping as $t_H \equiv \tilde{t}_H/\sqrt{z}$, with $\tilde{t}_H$ finite as $z \to \infty$. This different scaling can be traced back to the fact individual fermionic operators cannot have a non-zero expectation value; i.e. $\langle \hat{c} \rangle = 0$ due to conservation of fermionic parity in physical states. As such the trivial mean-field decoupling of the hopping—discussed above for bosons and which naturally led to the $1/z$ scaling—does not apply. As discussed in Ref. [350], this key difference changes the structure of the large connectivity limit.

To give a concrete example we consider a Fermi–Hubbard-like mode with next-nearest neighbour hopping and a local on-site interaction, with Hamiltonian:

$$ \hat{H} = -\frac{\tilde{t}_H}{\sqrt{z}} \sum_{\langle i,j \rangle} \sum_{\sigma} (\hat{c}_{i\sigma}^{\dagger} \hat{c}_{j\sigma} + \text{H.c.}) + \sum_{i\sigma} \hat{h}_{\text{loc}}[\hat{c}_{i\sigma}^{\dagger}, \hat{c}_{i\sigma}]. \quad (118) $$

In addition we can consider a set of local jump operators $\hat{L}_i$, describing losses, dephasing or other local dissipative processes. Writing the Keldysh action associated to this Lindblad master equation one can follow the steps discussed above and obtain the effective quantum impurity model. The main difference with respect to the bosonic case is that now the coherent field— related to the expectation value of the impurity field—vanishes, and only the effective non- Markovian bath $\Lambda^{\zeta\zeta'}(t, t')$ survives. As previously, this bath must be found self-consistently by computing the Green’s function of the fermionic impurity. A broad array of methods to solve these fermionic impurity models in and out of equilibrium have been developed, includ- ing for example exact diagonalisation [356], diagrammatic Monte Carlo [374–378], tensor- networks methods [379–381], auxiliary Master Equation methods [208] or hybridisation ex- pansion methods [382, 383]. We note that for a specific type of dissipation, corresponding to linear coupling to a quadratic fermionic bath, one can proceed exactly by integrating out the environment and retain a full non-Markovian description of the dissipation. This approach has been used in connection with DMFT in several works, see for example Ref. [384–387]. A different approach to dissipation within DMFT has been explored via coupling to a bosonic (phonon) bath, see Refs. [388–390].

## 4.6 Systematic expansions beyond mean-field theory.

There are a number of approaches by which one may systematically improve on the results of mean-field approximations. By “systematic”, we refer here to methods in which one can vary some convergence parameter, $P$, so as to interpolate between mean-field theories (as described above) in the limit $P = 1$, and numerically exact approaches in the limit $P \to \infty$ (or above some large but finite $P$ for a finite system). Each of these approaches described below corresponds to a systematic expansion in this sense and they all differ in the way they organise what kinds of correlations are kept in the calculation.

Within such expansions there are two further questions that determine the usefulness of an approach: first, is there uniform convergence on the exact solution? and second, is a sufficient level of accuracy achievable in practicable calculations? Regarding uniformity of convergence,
<figure><img src="image_18.png" /><figcaption>Figure 15: Illustration of cluster mean-field theory. The four sites illustrated are considered as a single cluster, while the effect of other sites (other clusters) are represented by their mean fields. From [283].</figcaption></figure>

this is defined by asking whether, for some error $\epsilon$, one can find a $P_\epsilon$ such that the error $E_p$ associated with calculations with parameter $P \ge P_\epsilon$ always satisfies $|E_p| < \epsilon$. It is worth noting that even when uniform convergence exists, there can be occasions where $E_p$ increases from one value of $P$ to the next—i.e. $E_p$ need not be monotonic. Uniform convergence though implies that for some sufficiently large $P_\epsilon$, all subsequent errors $E_p$ are sufficiently small. The fact that convergence need not be monotonic does imply the need for care when determining whether numerical results do or do not imply a converged result.

Regarding the practicability of calculations, we note that in some cases (such as matrix-product-state approaches) it has often been possible to reach convergence to reasonable error thresholds $\epsilon$. In such cases, one essentially has an exact result. In other cases, such as cluster mean field theory or cumulant expansions, reaching convergence has not generally been possible: the computational resources required grow too fast with $P$, or the convergence of error with $P$ is slow or non-uniform. Nonetheless, each of these approaches shares the idea that in principle, tuning a control parameter of the method (as opposed to changing the problem) can lead to exact results. We next discuss in turn various examples of such expansions.

In most of the cases discussed below, one can formulate the approximation both in terms of a recipe of how to perform the calculation, and also in terms of a variational Ansatz. The variational Ansatz formulation can be useful in understanding exactly what approximation is being considered: what forms of correlations are kept and what are lost. In principle, once one has formulated such an approach, one can then use the Dirac–Frenkel time-dependent variational principle [335,336] as introduced at the start of Sec. 4.4 to obtain a set of equations of motion for time evolution. In some cases there are however more direct approaches of how to perform time evolution, or to determine the steady state.

### 4.6.1 Cluster mean field theory.

In this approach, starting from the real-space mean-field theory described above, one systematically expands the size of what is considered as a “site” in the mean-field factorisation. The convergence parameter $P$ here is the size of the cluster. This means that one considers a cluster of several bare sites as one effective site, as shown in Fig. 15(a), and performs a factorisation of correlations between these different sites [283]. In the limit of an infinitely large cluster, this approach corresponds to including all correlations and so becomes exact. However, the cost of calculations increases exponentially with cluster size due to the growth of the size of Hilbert space. Such an approach would work well in cases where correlations are short ranged, so that including correlations between neighbouring sites captures these effects.

### 4.6.2 Tensor-network methods.

Tensor network methods provide a systematic expansion in terms of the degree of correlation between different “sites” in a system, using a structure which is efficient if the degree of cor-
relation remains small. Such methods are particularly effective in one dimension, in which case the tensor network simplifies to a structure known as a matrix product state (MPS). In the following we will use this one-dimensional matrix-product-state approach to illustrate the key ideas, before expanding to discuss cases beyond one dimension. There are many reviews covering both matrix product states and more general tensor networks, including [391–394]. Even within the one-dimensional MPS approach, there are multiple ways the concept of an MPS can be used to model a many-body open quantum system. As such we introduce the general concept first before discussing the different ways it can be applied to open systems.

**Introduction to matrix product states.** The essential point of all matrix-product-state methods is that a tensor with many indices (i.e. of high order), $T_{v_1, v_2, \dots, v_{N_s}}$ can always be written as a product of lower-order tensors:

$$ T_{v_1, v_2, \dots, v_{N_s}} = \sum_{\mu_1, \dots, \mu_{N_s-1}} A_{(1)}^{\mu_1} A_{(2)}^{\mu_1, \mu_2} \dots A_{(N_s-1)}^{\mu_{N_s-2}, \mu_{N_s-1}} A_{(N_s)}^{\mu_{N_s-1}} A_{v_{N_s}}^{\mu_{N_s}}, \quad (119) $$

as illustrated in Fig. 16(a). Let us assume the indices $v_i$ run over the values $v_i = 1, \dots, D$. The key benefit of such a rewriting can be seen by considering what happens if the indices $\mu_1, \mu_2, \dots, \mu_{N-1}$ can be restricted to run over a small number of values. If one writes $\mu_i = 1, \dots, M_i$ then $M_i$ is known as the bond dimension of the $i$th bond. One may then note that while specifying the original tensor requires $D^{N_s}$ values, specifying the matrix product state form requires $\sum_{i=1}^{N_s} DM_{i-1}M_i$ values (where we take $M_0, M_{N_s} \equiv 1$). If one can truncate $M_i$ to a small value, this replaces exponential scaling with $N_s$ with polynomial scaling. The maximum bond dimension $M_i$ is the convergence parameter for this approach.

In principle one can construct the MPS decomposition of a high-order tensor directly by performing sequential singular value decompositions (SVDs). For example, one can combine indices in the tensor $T$ to write $T_{v_1, \bar{v}_1}$ where $\bar{v}_1$ is the combination of all remaining indices, $v_2, \dots, v_{N_s}$. The singular value decomposition then decomposes this matrix as:

$$ T_{v_1, \bar{v}_1} = \sum_{\mu_1} [U]_{v_1, \mu_1} \sigma_{\mu_1} [V^{\dagger}]_{\mu_1, \bar{v}_1} = \sum_{\mu_1} A_{(1)}^{\mu_1} T_{(1)}^{\mu_1} \bar{v}_1. \quad (120) $$

Here $U, V$ are unitary matrices and $\sigma$ are (positive) singular values of $T$. The truncation of the sum over $\mu_1$ is made by choosing the largest $M_1$ values $\sigma_{\mu_1}$; these describe the most significant forms of correlation between the first index and the other indices. Depending on the application, there are different ways to identify the tensors $A_{(1)}, T_{(1)}$ from the SVD. These include the symmetric form that identifies $A_{(1)}^{\mu_1}, v_1 = [U]_{v_1, \mu_1} \sqrt{\sigma_{\mu_1}}, T_{(1)}^{\mu_1} = \sqrt{\sigma_{\mu_1}} [V^{\dagger}]_{\mu_1, \bar{v}_1}$, or left-canonical form where $A_{(1)}^{\mu_1}, v_1 = [U]_{v_1, \mu_1}$ and $\sigma$ is included entirely in $T$. The combined index $\bar{v}_1$ can then be split into $v_2, \bar{v}_2$, and the process repeated to split $T_{(1)}$ into tensors $A_{(2)}, T_{(2)}$, and so on until the entire tensor has been split up. We note however that such a direct decomposition of a state is rarely useful: it requires that one already has calculated the exponentially large object $T$. Instead practical algorithms discussed below work directly with the compressed representation, i.e. the tensors $A_{(i)}$.

While the approximation of a tensor by an MPS cannot hold in general, it turns out that in problems where correlations are bounded and relatively local, such an approximation can be good. That is, in many problems, the set of physically relevant states are all states that can be represented by matrix product states with small bond dimension. One may also note that in the limit $M = 1$, where all indices $\mu_i$ are restricted to take a single value, this state becomes a product state, corresponding to the mean-field factorisation Ansatz.

The more-general term “tensor network states” captures any Ansatz based on contracting networks of tensors to represent an object; we discuss some examples of this below.
<figure><img src="image_19.png" /><figcaption>Figure 16: (a) Matrix product state decomposition of a high-order tensor T. The six-index tensor (top) can be decomposed into multiple three-index tensors (bottom) as discussed in the text. (b) Matrix product operator form of a density matrix. (c) Equivalent matrix product state form after vectorisation. (d) Purification Ansatz tensor network for a density matrix.</figcaption></figure>

**Application to open quantum systems.** The matrix product state decomposition of a tensor can apply to many objects: wavefunctions, classical probability distributions, or density matrices [281,395]. There are many algorithms that have been developed to find or manipulate MPS for wavefunctions: we mention here two key algorithms. Firstly, the density matrix renormalisation group (DMRG) [396]—which actually predates the MPS formalism [397], but which is now understood in terms of MPS [392]—is an algorithm to perform variational minimisation of an MPS to find the ground state of a Hamiltonian. Secondly, time-evolving block decimation (TEBD) is an algorithm to time-evolve an MPS state with nearest-neighbour interactions.

To apply these ideas to density matrices, there are different ways one may proceed. One approach, illustrated in Fig. 16(b), is to represent these as a matrix product operator (MPO). An MPO differs from an MPS in that the object on each site has two indices, corresponding to an operator, rather than one operator, corresponding to a state. The MPO approach is often used to represent operators that act on quantum states, and one can in principle use the same concept to represent a density matrix. An alternative approach is to vectorise the density matrix as discussed in Sec. 3.1.2. This corresponds to writing the many-site density matrix as

$$ \hat{\rho} = \sum_{v_1, v_2, \dots, v_{N_s}} \rho_{v_1, v_2, \dots, v_{N_s}} \hat{\tau}_{v_1} \otimes \hat{\tau}_{v_2} \otimes \dots \otimes \hat{\tau}_{v_{N_s}}, \quad (121) $$

where $\hat{\tau}_{v_i}$ are a basis for the density matrix of site i. One then writes the tensor $\rho_{v_1, v_2, \dots, v_{N_s}}$ in the form of Eq. (119), where the label $v_i = 1, \dots, d^2$ runs over the space of local density matrices corresponding to a d-dimensional Hilbert space on a given site. This vectorised approach is illustrated in Fig. 16(c). Even within this vectorised approach, there remain different choices as to the density matrix basis: One can form this density-matrix basis from outer products of pure basis states, $\tau_v = |\psi_{l_v} \rangle \langle \psi_{r_v}|$, or one can use a basis such as the Pauli basis for two-level systems and their n-level generalisation through generalised Gell-Mann matrices.

A third approach to use MPS (or more precisely tensor-network) methods for open quantum systems is the purification approach [398]—this constructs an explicit Ansatz for how the density matrix can be written as a sum over pure matrix product states, and is illustrated in Fig. 16(d).

As noted above for closed-system wavefunctions, there are a variety of algorithms that one can apply to find or manipulate MPS representations of states; the same is true for MPS or MPO representations of density matrices. The TEBD algorithm for time evolution can be directly
adapted to time evolution of a density matrix, including both Hamiltonian and Lindblad terms. If all sites are distinct, then application of this approach has a computational effort that grows (albeit polynomially) with the number of sites. However, for translationally invariant problems (with periodic boundary conditions) one can assume all site tensors are equivalent, and directly access the infinite system limit—this approach is often referred to as infinite time-evolving block decimation (iTEBD) [399].

There also exist approaches that seek to directly find the steady state of the open quantum system—analogous to DMRG approaches for finding the ground state of a closed system. For the open system after vectorisation, rather than minimising the expectation of the Hamiltonian, one must seek a solution of the steady-state condition $\hat{L}|\rho\rangle = 0$. This can be turned into an optimisation problem by considering the minimisation of either $|\langle\rho|\hat{L}|\rho\rangle|$ or of the Hermitian version $\langle\rho|\hat{L}^\dagger\hat{L}|\rho\rangle$. One can then apply standard methods to minimise these expressions. For example, Ref. [400] performs such minimisation on the first of these forms, and Ref. [401] on the Hermitian version.

When considering the computational cost of MPS calculations for open quantum systems compared to closed systems a first consideration is that the larger state space—i.e. the $d^2$ basis states for a density matrix vs the $d$ states for a wavefunction—would make MPS calculations for open systems more costly. However, the effects of dissipation and dephasing in an open quantum system mean that entanglement (and correlations more generally) can be suppressed, so one may potentially have a lower bond dimension in the open system, despite the larger state space. This has been discussed in various contexts, including recently for the question of how qubit errors suppress the bond dimension needed to simulate quantum circuits [402,403]. It is though worth noting here another distinction between the closed system and open system MPS. For a closed system, the bond dimension relates directly to the entanglement entropy, which is given by the von-Neumann entropy of the singular values (i.e. Schmidt coefficients [190]) of the quantum state. For the open system, there can be both classical and quantum correlations [404], and both kinds of correlation contribute to the required bond dimension. It is also worth considering how the bond dimension changes over time. In the closed system, non-local interactions in the Hamiltonian lead the bond dimension to grow, so that there is typically an “entanglement barrier” to evolving to late times. For the open system, dissipation competes against this as discussed above. In some cases this can lead to non-monotonic dependence of entanglement on time when simulating dynamics following a quench of parameters or evolution from a trivial state: entanglement initially increases due to interactions before being suppressed due to dissipation. This means that for some open quantum systems it can be preferable—i.e. required bond dimension can be lower—to find steady states either by variational approaches, or by adiabatically varying parameters to explore the steady state phase diagram.

**Combining quantum trajectories and matrix product states.** As an alternative to studying density matrices, one may also proceed by unravelling the Lindblad dynamics and then combine MPS with quantum jumps and study the evolution of the wavefunction along the different quantum trajectories [269]. Following the discussion above, one may note that in the unravelled case, the bond dimension depends only on the entanglement within a given trajectory: classical correlations do not contribute to the bond dimension. As compared to closed system dynamics, bond dimensions however can be smaller, as jump operators describing dephasing processes will reduce entanglement. As noted recently [262, 263, 405], different unravellings can lead to quite different bond dimensions in the resulting matrix product state. In some cases, the bond dimension for the unravelled approach can significantly exceed that for the density matrix approach [406]. Another issue worth noting with the unravelled system is that when considering dissipation acting on each site individually, unravelling breaks translational
<figure><img src="image_20.png" /><figcaption>Figure 17: (a) Matrix product state representation imposed on a 2D lattice, following a serpentine path. The diagonal lines indicated the “exposed” indices describing the state on a given site. (b) Tensor network (e.g. projected entangled pair state) for a 2D lattice.</figcaption></figure>

invariance as the measurement records for each site can be distinct. This prevents one using iTEBD approaches as mentioned above.

**Beyond one dimension.** The discussion in this section so far has focussed on one dimensional lattice systems. However, the concepts discussed above can be applied or extended to consider other geometries. Two broad approaches exist: application of MPS forms to other geometries, or extension to consider tensor networks. In the former case, one may consider a lattice of arbitrary geometry and then apply the same MPS decomposition as written in Eq. (119). For example, this can be done as shown in Fig. 17(a), following a serpentine path on a two-dimensional lattice. This will however lead to complications in both time evolution and minimisation, as the structure of the MPS no longer matches the connectivity of the lattice. This causes a problem for time evolution as the interaction will not be nearest-neighbour in terms of the MPS representation. The alternative approach is to generalise from the one-dimensional matrix product state to a network of tensors—a tensor network state. An example of this is shown in Fig. 17(b), corresponding to projected entangled pair states (PEPS) in two dimensions [391, 393, 394]. These involve decomposing the many-body density matrix into a form that has a fifth order tensor on each site, with four internal legs connecting each site to its neighbours. For applications to open systems, see e.g. [331, 407, 408].

In addition to lattice problems, there has also been work on models with one-to-many coupling such as the Dicke model. In cases where the “many” sites—i.e. the satellite sites—are bosonic modes, one can use “chain mapping” techniques [409–412] to produce an effective one-dimensional chain model with nearest-neighbour interactions. Such an approach is based on the idea that for any spectral density arising from the set of *N* satellite modes can be recovered by engineering the on-site energies and couplings of a one-dimensional chain [410]. After such a mapping, one may then use MPS/MPO methods to time evolve the system.

In another set of works [413, 414], a more direct mapping to a one-dimensional chain was made for a model of *N* spins coupled to a bosonic cavity mode. Here the sites on the chain correspond directly to the *N* spins and the bosonic mode. This would lead to non-local (i.e. not nearest neighbour) interactions between each spin and the bosonic mode—as such this is similar using an MPS state to represent a two-dimensional lattice as discussed above. To overcome this non-locality, swap operations were repeatedly used to move the cavity site to be adjacent to each spin site in turn. Since the matrix-product-state geometry here does not
<figure><img src="image_21.png" /><figcaption>Figure 18: Restricted Boltzmann Machine Ansatz for a wavefunction, i.e. Eq. (122). Plotted for $N_s = 6$ input states (red circles) and $M = 4$ hidden nodes (blue squares).</figcaption></figure>

reflect the physical geometry, there can be significant correlations between neighbouring sites required to capture the correlations between each satellite site and the central cavity, leading to large bond dimensions $M$.

### 4.6.3 Neural network Ansatz

The tensor network approach described above makes clear the nature of the density matrix as a map from a pattern of indices to a complex number. The structure of the tensor network reflects how correlated the dependence on the indices should be. Recognising this structure suggests the possibility of alternative approaches, where a neural network is used to represent the map from input indices to the output value.

Neural-network approaches were introduced for ground-state wavefunctions in Ref. [415], specifically considering a Restricted Boltzmann Machine (RBM) Ansatz. An RBM is a two-layer network, with links only between (not within) the layers. One layer, the “visible” layer consists of the input indices, $v_1, ..., v_{N_s}$. For a wavefunction, these input indices correspond to the states of a set of a lattice of sites, i.e. $|\Psi\rangle = \sum_{v_1, v_2, ..., v_{N_s}} \Psi_{v_1, v_2, ..., v_{N_s}} |v_1\rangle |v_2\rangle \cdots |v_{N_s}\rangle$. We assume the input indices should take the form $v_i = \pm 1$, corresponding to the states of a spin 1/2 on each site. The other “hidden” layer, $h_1, ..., h_M$, with $M$ nodes, is traced over. For a wavefunction one thus writes:

$$ \begin{aligned} \Psi_{v_1, v_2, \dots, v_{N_s}} &= \sum_{h_1, \dots, h_M = \pm 1} \exp \left( \sum_{i=1}^{N_s} a_i v_i + \sum_{j=1}^{M} b_j h_j + \sum_{i,j=1}^{N_s M} W_{ij} v_i h_j \right) \\ &= \exp \left( \sum_{i=1}^{N_s} a_i v_i \right) \prod_{j=1}^{M} 2 \cosh \left( b_j + \sum_{i=1}^{N_s} W_{ij} v_i \right). \end{aligned} \quad (122) $$

See Fig. 18 for a sketch of this with $N_s = 6$, $M = 4$ showing the network before tracing out the hidden layer. The parameters $a_i, b_j, W_{ij}$ are all variational parameters, and various standard methods (e.g. stochastic gradient descent [416] or stochastic reconfiguration [417]) can be used to minimise the expectation of the Hamiltonian for such a state.

This RBM approach was later extended to density matrices. In Refs. [285–287, 418] an extension of the above Ansatz was made to represent a purification of the density matrix (thus guaranteeing its positivity). Such an Ansatz is written in terms of left indices $v_1^{(l)}, ..., v_{N_s}^{(l)}$ and right indices $v_1^{(r)}, ..., v_{N_s}^{(r)}$ for the density matrix. There are three hidden layers: one connected to just the left indices, one to the right indices, and one connecting the two indices together. To ensure Hermiticity the weight matrices connecting the left and right layers are related by complex conjugation, see Fig. 19. We thus have variational parameters $b_j, W_{ij}$ and their complex conjugates for the left and right layers respectively, with $i = 1, ..., N_s$ and
<figure><img src="image_22.png" /><figcaption>Figure 19: Restricted Boltzmann Machine Ansatz for a density matrix, Eq. (123). Plotted for $N_s = 6$ input states and $M_1 = M_2 = 4$, i.e. four nodes for both types of hidden layer. The hidden layers (blue squares) are not explicitly labelled as Eq. (123) presents the result after summing over these layers. Note the structure here has a set of hidden layers (on the left) coupled only to the left (ket) indices, $v_i^{(l)}$ of the density matrix, a set of hidden layers (on the right) coupled only to the right (bra) indices $v_i^{(r)}$, and a set of hidden layers in the middle that couple the two sets of indices.</figcaption></figure>

$j = 1, \dots, M_1$, and parameters $c_k, X_{ik}$ for the connecting layers, where $k = 1, \dots, M_2$. After tracing out the hidden layers this Ansatz takes the form:

$$ \begin{aligned} \rho_{(v_1^{(l)}, \dots, v_{N_s}^{(l)}), (v_1^{(r)}, \dots, v_{N_s}^{(r)})} = & \exp \left( \sum_{i=1}^{N_s} a_i v_i^{(l)} + a_i^* v_i^{(r)} \right) \\ & \times \prod_{j=1}^{M_1} 4 \cosh \left( b_j + \sum_{i=1}^{N_s} W_{ij} v_i^{(l)} \right) \cosh \left( b_j^* + \sum_{i=1}^{N_s} W_{ij}^* v_i^{(r)} \right) \\ & \times \prod_{k=1}^{M_2} 2 \cosh \left( c_k + \sum_{i=1}^{N_s} (X_{ik} v_i^{(l)} + X_{ik}^* v_i^{(r)}) \right). \end{aligned} \quad (123) $$

Such a variational Ansatz describes a restricted subset of the Liouville space. As discussed in the previous section, in contrast to minimising the ground state energy for a closed system, variational minimisation for open quantum systems presents various choices. Just as with the MPS Ansatz discussed above, the steady-state condition $\hat{\mathcal{L}}[\hat{\rho}] = 0$ can be turned into an optimisation problem minimisation of either $|\langle \rho | \hat{\mathcal{L}} | \rho \rangle|$ or the Hermitian version $\langle \rho | \hat{\mathcal{L}}^\dagger \hat{\mathcal{L}} | \rho \rangle$. These then allow the application of standard variational methods to optimise the RBM Ansatz, such as stochastic gradient descent and stochastic reconfiguration.

One alternative way to write an RBM Ansatz is to consider vectorisation of the indices on each site individually, and construct an RBM as a function of these indices [419, 420]. Because these vectorised indices now have more than two states, one cannot directly use the Ansatz as written previously. However, as discussed in [420], if one labels the four states $\{v_i^{(l)}, v_i^{(r)}\} = (++, +-, +-, --)$ via a combined index $v_i = (+2, +1, -1, -2)$ one can write:

$$ \rho_{v_1, \dots, v_{N_s}} = \exp \left( \sum_{i=1}^{N_s} \sum_{\xi=1}^{3} a_i^{(\xi)} (v_i)^{\xi} \right) \prod_{j=1}^{M} 2 \cosh \left( b_j + \sum_{i=1}^{N_s} \sum_{\xi=1}^{3} W_{ij}^{(\xi)} (v_i)^{\xi} \right), \quad (124) $$

involving three different powers of the indices. This form of Ansatz is better able to efficiently describe products of mixed states than the form in Eq. (124).

As defined so far, this Ansatz does not have an obvious control parameter to systematically improve the approximation. The role of the convergence parameter here corresponds to the number of hidden neurons (or more specifically, the ratio of hidden neurons to visible neurons
$M/N_s$). Increasing $M/N_s$ increases the expressivity of the Ansatz, and convergence to the known exact result on increasing $M/N_s$ has been seen [415]. Another possibility one may also consider is to add further hidden layers to the neural network, i.e. going beyond the restricted Boltzmann machine.

Another way to describe open quantum systems using an more general neural networks has been to first make use of an alternate representation of the density matrix via the outcomes of positive operator-valued measures (POVMs) [190]. That is, rather than time evolving the system density matrix, one considers quantities $P^a = \text{Tr}[\hat{\rho}\hat{M}^a]$ where $\hat{M}^a = \hat{M}^{a_1} \otimes \hat{M}^{a_2} \cdots \otimes \hat{M}^{a_{N_s}}$, and each $\hat{M}^{a_i}$ is a POVM on site $i$. In such approaches, the key idea is that $P^a$ is a probability distribution which can then be represented by any form of neural network suitable to describe probabilities [421–423]. An alternate recent approach has been to use a form of purified density matrix via a Gram–Haddamard matrix [424].

### 4.6.4 Corner-space renormalisation

One further approach to truncating the full Hilbert space was introduced in Ref. [282] as “Corner-space renormalisation”. This approach makes use of an iterative construction, first solving a problem on a small lattice, and then using that to determine an effective basis for a larger lattice.

This method can be understood by considering one step of the process. Suppose one has found a good approximation to the steady state of the master equation $\hat{\rho}_A$ on subsystem A, and $\hat{\rho}_B$ on subsystem B.<fnref n="3" /> The aim is then to find an approximation on the system $A \cup B$. The density matrices can be written as

$$ \hat{\rho}_A = \sum_{\mu} p_{\mu}^{(A)} |\phi_{\mu}^{(A)}\rangle\langle\phi_{\mu}^{(A)}|, \quad \hat{\rho}_B = \sum_{\nu} p_{\nu}^{(B)} |\phi_{\nu}^{(B)}\rangle\langle\phi_{\nu}^{(B)}| \quad (125) $$

The corner-space approach is to construct a basis for the Hilbert space $\mathcal{H}_{A \cup B}$ by considering the $M$ states $|\phi_{\mu}^{(A)}\rangle|\phi_{\nu}^{(B)}\rangle$ which have the largest values of $p_{\mu}^{(A)}p_{\nu}^{(B)}$. That is, one determines a set of indices $\mu_m, \nu_m$ for $m = 1, \dots, M$ such that $p_{\mu_m}^{(A)}p_{\nu_m}^{(B)} > p_{\mu_{m+1}}^{(A)}p_{\nu_{m+1}}^{(B)}$, and then uses these to construct the basis:

$$ C(M) = \{ |\phi_{\mu_1}^{(A)}\rangle|\phi_{\nu_1}^{(B)}\rangle, |\phi_{\mu_2}^{(A)}\rangle|\phi_{\nu_2}^{(B)}\rangle, \dots, |\phi_{\mu_M}^{(A)}\rangle|\phi_{\nu_M}^{(B)}\rangle \}. \quad (126) $$

One then solves the master equation within this reduced “Corner-space” basis $C(M)$. The parameter $M$ serves as a convergence parameter, increasing $M$ until the solution is converged. This then yields a good approximation for the system $A \cup B$; the whole procedure can then be repeated to consider a larger system. We note that for this approach, one cannot straightforwardly formulate the approximation via a variational Ansatz for the density matrix. One may note that the basic idea—of truncating a basis formed by combining the bases of two subparts of the system—is closely related to the original formulation of the density matrix renormalisation group (DMRG) [396, 397]. We note that DMRG approaches can be reformulated in terms of matrix product states [392], no similar reformulation is yet known for corner-space renormalisation.

### 4.6.5 Cumulant expansion.

Cumulant expansions provide an alternate approach to systematically include correlations between different parts of the system, by systematically including correlations between $M$ sites,

<footnote n="3">³For generality this allows cases where $\hat{\rho}_A \neq \hat{\rho}_B$, however in translationally invariant problems, one expects $\hat{\rho}_A = \hat{\rho}_B$, since these are solutions to the same problem.</footnote>
$\langle \hat{O}_1 \hat{O}_2 \cdots \hat{O}_M \rangle$, where $\hat{O}_i$ is some operator on site $i$. One typically finds that the equations of motion for such correlations will involve terms arising from $M+1$ or more sites. These can be related to correlations of $M$ or fewer sites by assuming that all cumulants beyond $M$th order vanish and using this to write the expectations in terms of lower order terms [425]. Thus, $M$ serves here as the convergence parameter. The case where $M=1$ is equivalent to the mean-field approximation described above. For example, for $M=2$, the 3rd order cumulant $\langle \hat{O}_1 \hat{O}_2 \hat{O}_3 \rangle_c$ is defined as

$$ \langle \hat{O}_1 \hat{O}_2 \hat{O}_3 \rangle_c = \langle \hat{O}_1 \hat{O}_2 \hat{O}_3 \rangle - \langle \hat{O}_1 \hat{O}_2 \rangle \langle \hat{O}_3 \rangle - \langle \hat{O}_1 \hat{O}_3 \rangle \langle \hat{O}_2 \rangle - \langle \hat{O}_2 \hat{O}_3 \rangle \langle \hat{O}_1 \rangle + 2 \langle \hat{O}_1 \rangle \langle \hat{O}_2 \rangle \langle \hat{O}_3 \rangle, \quad (127) $$

and by setting it to zero, $\langle \hat{O}_1 \hat{O}_2 \hat{O}_3 \rangle_c = 0$, allows to obtain the correlation of 3 sites in terms of those of fewer sites. In the limit of $M \to \infty$ all correlations are kept, so there is no approximation.

For open systems, the Ansatz underpinning such an expansion is to assume a generalisation of density matrix factorisation, allowing two point correlations, e.g. for equivalent sites one may consider:

$$ \hat{\rho} = \bigotimes_{i} \hat{\rho}_{i}^{(1)} + \sum_{i \neq j} \rho_{ij}^{(2)} \bigotimes_{k} \hat{\rho}_{k}^{(1)} + \sum_{i \neq j \neq k \neq l} \rho_{ij}^{(2)} \otimes \rho_{kl}^{(2)} \bigotimes_{m \neq i,j,k,l} \hat{\rho}_{m}^{(1)} + \dots, \quad (128) $$

where $\rho_{ij}^{(2)}$ are the traceless part of two-site density matrices. The second and third terms written here include correlations within a single pair of sites and within two pairs of sites respectively. The standard second-order cumulant approach (as discussed above) is equivalent to including all possible pairwise correlations. This would mean that e.g. a fourth-order term $\langle \hat{O}_1 \hat{O}_2 \hat{O}_3 \hat{O}_4 \rangle$ would have an expansion that includes terms with two pairwise correlations, consistent with what would be found by setting all 3rd and 4th order cumulants to zero.

This approximation has been considered in the context of problems with one-to-many and many-to-many couplings. For such problems, as noted above, the matrix product state approach is challenging. In such cases, where all satellite sites are equivalent, one may note that the computational effort to model $N$ sites at a given order of cumulant expansion is independent of $N$: the value of $N$ appears only as a parameter in the equations, and does not affect the number of equations to be solved. The cumulant expansion has been used also in the context of driven-dissipative spin chains, see for example Ref. [426] where it was benchmarked with matrix product state calculations. An open source implementation of the cumulant expansion to open quantum systems exists in Ref. [427]; see references therein for other applications.

As discussed so far, cumulant methods have been applied to simplify the master equation for the ensemble-averaged density matrix. Cumulant methods have also been combined with stochastic unravellings, to capture correlations within a single quantum trajectory [428].

### 4.6.6 Multi-configuration time-dependent Hartree (MCTDH).

This approach builds on the particle-based—i.e. Hartree–Fock—version of mean-field theory. As noted above, such a mean-field theory corresponds to considering a single configuration of single-particle states. to be occupied by the particles.

The multi-configurational approach is based on extending to consider a number $M > 1$ of different configurations, and increasing the parameter $M$ until numerical convergence is reached. This was originally introduced for describing wavefunctions of distinguishable particles [429,430], and later extended to indistinguishable fermions [431–433] and bosons [434, 435]. In the context of bosons the MCTDH-B Ansatz takes the form:

$$ |\Psi(t)\rangle = \sum_{n_1, n_2, \dots, n_M} C_n(t) \prod_{i=1}^{M} \frac{[\hat{b}_i^{\dagger}(t)]^{n_i}}{\sqrt{n_i!}} |0\rangle, \quad \hat{b}_i^{\dagger}(t) = \int d^d r \psi_i(r, t) \hat{a}(r). \quad (129) $$
Here the vector **n** defines the (time-dependent) occupation of the different configurations, while $\psi_i(\mathbf{r}, t)$ describes the spatial profile of these time-dependent configurations. These different configurations are constrained to remain orthogonal, and then time evolution of $C_n(t)$, $\psi_i(\mathbf{r}, t)$ is found by using the Dirac–Frenkel time-dependent variational principle [335, 336]. In the limit $M = 1$ this recovers the Ansatz that leads to the Gross–Pitaevskii equation.

There are various approaches to extending MCTDH to density matrices [436, 437]. An open-source implementation exists [438] and has been used to study the dynamics of cold atoms in optical cavities [127].

## 4.7 Schwinger–Keldysh path integral methods.

In this section we discuss various approaches to solving problems based on the Schwinger–Keldysh formalism, as introduced in Sec. 3.2. For “non-interacting” problems, such as Eq. (54), where the action is quadratic in the fields, one can directly use this equation to evaluate any correlation functions of interest. For example, for the non-interacting single bosonic mode problem discussed above, one can calculate two point functions by inverting the the matrix in Eq. (54) to give

$$ \begin{pmatrix} D_K & D_R \\ D_A & 0 \end{pmatrix}, \qquad D_K = -D_R [D^{-1}]_K D_A, \qquad (130) $$

where $D_{R,A}$ are the simple inverses of the matrices above. In frequency space this means $D_{R,A}(\nu) = [\nu - \omega \pm i\frac{\omega}{2}(\kappa_- - \kappa_+)]^{-1}$. By Fourier transforming back to the time domain these Green's functions give (for bosons):

$$ D_R(t) = -i\theta(t) \langle [\hat{a}(t), \hat{a}^\dagger(0)] \rangle, \qquad D_K(t) = -i \langle \{\hat{a}(t), \hat{a}^\dagger(0)\} \rangle, \qquad (131) $$

where $\theta(t)$ is the Heaviside step function. The two forms given here correspond to the commutator and anti-commutator respectively, thus allowing any particular operator ordering to be determined. For Gaussian problems (i.e. where the action has only quadratic terms), all higher-order correlation functions can be factorised into products of two-point correlation functions.

Physically the retarded Green's function corresponds to a response function—hence its causal structure, with non-zero value only for $t > 0$ in the time domain, and all poles in the lower half-plane in the frequency domain (assuming the stable configuration, $\kappa_- > \kappa_+$). In distinction, the Schwinger–Keldysh Green's function describes the noise spectrum. In a thermal equilibrium system, these noise and response functions are related by

$$ D^K(\nu) = D^R(\nu)F(\nu) - F(\nu)D^A(\nu), \qquad (132) $$

where $F(\nu)$ is a distribution function, where $F(\nu) = 1 \pm 2n_{B,F}(\nu)$ for bosons or fermions respectively, with $n_{B,F}$ being the Bose or Fermi occupation function. This identity is a form of the fluctuation dissipation theorem [55, 56](FDT), relating $D^K(\nu)$; it is also the frequency domain version of the Kubo–Martin–Schwinger (KMS) [439, 440] condition on equilibrium Green's functions. In a general nonequilibrium context $F(\nu)$ will not adopt an equilibrium form, however one can often extract a low-frequency effective temperature [231, 441] by fitting the behaviour as $\nu \to 0$; e.g. for a thermal bosonic system, $F(\nu) = \coth(\beta \nu/2) \approx 2k_B T/\nu$, so one may define $k_B T_{eff} \equiv \lim_{\nu \to 0} \nu F(\nu)/2$.

For interacting problems, reformulating the problem in terms of such a path integral is not itself a method of solution, but provides a framework to then consider a number of theoretical approaches. In the remainder of this section we provide a summary of some of these approaches.
### 4.7.1 Saddle point and mean field theory.

One can often determine the value of a path integral approximately by considering the saddle point configuration—the path $\psi_{cl,q}(t)$ which makes the action stationary—and then expanding around this. One may show [54] that a “classical saddle point” will generally exist—i.e. a solution with $\psi_q(t) = 0$—at least in the absence of quantum sources, i.e. if $J_q = (J_f - J_b)/\sqrt{2} = 0$. Quantum saddle points may also exist when considering boundary conditions with different initial and final states: these can describe tunnelling as well as thermal activation connecting different initial and final states [54].

As with any path integral, this saddle-point approach coincides with a form of mean-field theory. For example, for the action for the weakly interacting dilute Bose gas given in Eq. (56), taking the saddle point with respect to the standard bosonic field recovers the complex Gross-Pitaevskii equation in Eq. (109). Specifically to recover this one should evaluate the variation $\delta S_{WIDBG}/\delta \psi_q^* = 0$ and then take $\psi_q \to 0$, corresponding to the classical saddle point. This yields:

$$i \frac{d}{dt} \psi_{cl} = \left( -\frac{\nabla^2}{2m} + \frac{U}{2} |\psi_{cl}|^2 + \frac{i}{2} \left[ \kappa_+ - \kappa_- - \frac{\kappa_-^{(2)}}{2} |\psi_{cl}|^2 \right] \right) \psi_{cl}. \quad (133)$$

To recover Eq. (109) one should make the identification $\psi_{cl} \to \sqrt{2N}\psi$, corresponding to difference in normalisation between the single particle wavefunction $\psi$ used to derive Eq. (109) and the Keldysh-rotated many-particle field $\psi_{cl}$ appearing in the path integral.

A similar approach can be formulated to recover the Gutzwiller mean-field theory of Sec. 4.4.1 as a saddle point. Starting from the Keldysh action of the lattice and performing the Hubbard-Stratonovich decoupling [254] of the hopping in terms of a field $\phi_i, \phi_i^*$ one can obtain the effective action

$$S_{\text{eff}} = \int_C dt \sum_{ij} \phi_i^* t_{ij}^{-1} \phi_j + \sum_i \Gamma[\phi_i^*, \phi_i]. \quad (134)$$

The second term here represents the generating functional of the Green functions of decoupled sites:

$$\Gamma[\phi_i^*, \phi_i] = -i \log \left\langle \mathcal{T}_C \exp \left( i \int_C dt (\phi_i^* \hat{a}_i + \hat{a}_i^* \phi_i) \right) \right\rangle_0, \quad (135)$$

where $\mathcal{T}_C$ indicates time ordering along the Keldysh contour $C$ shown in Fig. 12. We note that this average is taken with respect to a density matrix where different sites are decoupled, therefore formally reducing the many-body problem to solving a driven-dissipative single site [368, 442]. By taking the saddle point of the effective action $\delta S_{\text{eff}}/\delta \phi_i = 0$ one obtains $\phi_j = \sum_i t_{ij} \langle \hat{a}_i \rangle$ thus recovering the Gutzwiller mean-field theory discussed in Sec. 4.4.1.

For problems involving fermions, one cannot take a saddle point with respect to Grassmann variables. In such cases one must first perform whatever Hubbard-Stratonovich decouplings are required to yield a non-interacting fermion problem, and then integrate out the fermions. This will yield an action in terms of only complex fields, for which one can then take a saddle point. As a general principle, to find observable properties using the saddle point approach requires several steps. One should first add source currents $J_{cl,q}$ as noted in Sec. 3, evaluate the saddle point in the presence of those currents, and then take derivatives of the action with respect to the currents to evaluate physical observables. In many cases however a simpler approach can be realised, of identifying the saddle point value of a field $\psi$ with the expectation of the corresponding operator. For cases where operator correlation functions are required however, the full approach of evaluating the saddle point in the presence of the classical fields.
4.7.2 Quadratic fluctuations and ladder series expansions.

After having evaluated a saddle point, one can then systematically expand around this. Con-
sidering fluctuations up to quadratic order yields an effective action for which one can still
exactly evaluate the path integral. As such, this allows one to find the effect of leading order
fluctuation corrections to the mean-field result.

When such an approach is used to calculate correlations or response functions, the result
of the saddle point and fluctuation approach is equivalent to re-summation of a diagrammatic
expansion, typically corresponding to a ladder series in terms of the bare Green's functions.
For example, this approach was discussed in Ref. [163] as a route to calculate the superfluid
stiffness of a driven-dissipative WIDBG model.

Similarly, by expanding the effective action $S_{\text{eff}}$ in Eq. (134) around the mean-field con-
figuration one can obtain the fluctuation contribution to the saddle point, controlled in this
case by the local Green's functions of the driven-dissipative single-site problem. These objects
can be computed exactly for example via a spectral decomposition of the single site Lindbla-
dian [443,444].

4.7.3 Renormalisation group.

The Schwinger–Keldysh action provides a useful starting point for many renormalisation group approaches [31, 32, 242, 445–447]. The renormalisation group is a very broad area with many books and reviews covering it; introductions to this topic include Refs. [254,448,449]. The key idea is that by integrating over “fast” fields (i.e. the components of fields at high frequencies or momentum) one can find the effective action for slow modes, i.e. those at long wavelength and low frequency. After eliminating such modes, one then rescales the remaining momenta and frequencies to restore the original structure of the action, but with modified coefficients for each term. In almost all cases, nonlinearity of the action will mean that the fast modes couple to the slow modes, so that integrating out fast modes modifies the coefficients in the action that determine the slow modes. The goal of such analysis is to determine what form the low energy theory takes. This approach can also be particularly useful to identify phase transitions and study their critical behaviour, by identifying different fixed points of the renormalisation group flow that define either phases or critical boundaries between different phases.

The simplest approach one may consider to explore the likely consequences of renormali-
sation is to use “power counting” to determine which terms in the Keldysh action are relevant
at long wavelengths and low energies [31, 32, 450]. This means looking at each term in the
action, and counting powers of momentum, *k*, associated with each term. For example, consid-
ering the WIDBG action in Eq. (56), one has the following considerations: Spatial derivatives
naturally scale as *k*¹, while the dispersion of the WIDBG means time derivatives scale as *k*².
As such, the integral over both time and space gives a scaling *k*^-(d+2) in *d* dimensions. The
overall action must be invariant, i.e. *k*⁰, as must the constant factors in the quadratic terms κ±.
From these considerations we see the quantum field must scale as ψ_q ~ *k*^((d+2)/2). Considering
the other quadratic terms in the action we see that ψ_cl ~ *k*^((d-2)/2). These considerations then
lead to a useful observation regarding the quartic terms: the dominant terms there are those
involving only one quantum field (and three classical fields). All other quartic terms will be
smaller by a factor of *k*², and thus smaller in the long wavelength *k* → 0 limit. This yields the
effective long wavelength action:

$$S_{WIDBG,SC} = \int_{t_0}^{t_1} dt \int d^d \mathbf{r} \left\{ (\psi_{cl}^* \quad \psi_q^*) \begin{pmatrix} 0 & [D_A^{(0)}]^{-1} \\ [D_R^{(0)}]^{-1} & [D^{(0)-1}]_K \end{pmatrix} \begin{pmatrix} \psi_{cl} \\ \psi_q \end{pmatrix} - \left[ \left( \frac{U}{2} + i \frac{\kappa}{4} \right)^2 \psi_{cl}^* \psi_{cl} \psi_q + \text{c.c.} \right] \right\}. \quad (136)$$

This simplified action can also be considered as a semiclassical limit.

The power counting approach will allow one to determine which terms can survive in a low energy theory, but does not determine how the renormalised coefficients relate to the original coefficients.

To perform more complete calculations, one must use an explicit scheme to integrate out degrees of freedom and see how coefficients change. One fruitful approach there is the functional renormalisation group [449], which has been applied to Keldysh actions [31, 32, 445]. This works by introducing a “regulator” into the effective action which suppresses long wavelength modes with $k < K$; one then studies an effective action as a function of $K$, and derives an equation for how the action evolves as $K$ changes. This ultimately yields equations for how the form of effective action changes with the momentum scale at which one works.

### 4.7.4 Stochastic interpretations.

As discussed in section 8, there are various methods by which one can transform a master equation for the ensemble averaged density matrix into equations for a set of individual stochastic trajectories. The Schwinger-Keldysh path integral—and in particular the semiclassical limit noted above—can be used as one such route to derive a stochastic equation of motion that captures the ensemble behaviour.

Following the notation of Ref. [450], one can consider deriving the stochastic equation by starting from a form of action that frequently arises in the semiclassical limit,

$$S = \int dt \int d^d r \left[ \psi_q^* \left( i \frac{d}{dt} \psi_{cl} + \frac{\delta f[\psi_{cl}]}{\delta \psi_{cl}} \right) + \text{c.c.} + i \Gamma |\psi_q|^2 \right]. \quad (137)$$

This “semiclassical” form allows a general dependence on the field $\psi_{cl}$, but restricts how $\psi_q$ appears. The $\psi_q$ dependence that occurs in this expression is in line with the terms that are relevant at long lengths scales according to the power-counting argument discussed above.

Starting from this action, the key idea is to perform a Hubbard-Stratonovich decoupling of the term $i\Gamma|\psi_q|^2$, by using the identity:

$$\exp\left[-\int dt \int d^d r \Gamma |\psi_q|^2\right] = \int \mathcal{D}(\xi) \exp\left[-\int dt \int d^d r \left(\frac{|\xi|^2}{\Gamma} + i\xi^* \psi_q + i\xi \psi_q^*\right)\right]. \quad (138)$$

This puts the action into a form that is linear in $\psi_q$, $\psi_q^*$. As such the path integral over $\psi_q$, $\psi_q^*$ introduces a Dirac delta function requiring that $\psi_{cl}$ obey the Langevin equation:

$$i \frac{d}{dt} \psi_{cl} = - \frac{\delta f[\psi_{cl}]}{\delta \psi_{cl}} + \xi, \quad (139)$$

The dependence of the action of $\xi$ implies it has a white noise spectrum:

$$\langle \xi(\mathbf{r}, t) \rangle = 0, \quad \langle \xi(\mathbf{r}, t) \xi^*(\mathbf{r}', t') \rangle = \Gamma \delta(t-t') \delta(\mathbf{r}-\mathbf{r}'). \quad (140)$$

The procedure here, extracting a Langevin equation from a Keldysh action is the inverse of the Martin-Sigia-Rose approach [54], which derives a Keldysh action starting from a Langevin
equation. The interpretation of the above equation is that for such semiclassical actions, the resulting correlation functions can be equivalently found by simulating a set of trajectories of $\psi_{cl}$ from the Langevin equation in Eq. 139, and then averaging over noise.

As discussed in Refs. [445,450], a further simplification occurs if the functional $f$ appearing above is purely imaginary. In such a case, one can give an explicit expression for the probability distribution of $\psi_{cl}$ that results from averaging over such trajectories:

$$P[\psi_{cl}] = \exp\left(\frac{if[\psi_{cl}]}{\Gamma/2}\right). \qquad (141)$$

(Note that the exponent here is real, since $f$ is imaginary.) In this case, despite the derivation from a nonequilibrium formalism, the behaviour can be understood as corresponding to an equilibrium Gibbs ensemble, with an energy controlled by the functional $-if[\psi]$, and a temperature set by the noise strength $\Gamma/2$.

Emergent equilibrium behaviour is in fact possible more widely, if $\mathrm{Re}[f[\psi]] = r \mathrm{Im}[f[\psi]]$ for some constant $r$; i.e. the real and imaginary parts are proportional to each other [31, 445, 451, 452]. The purely imaginary functional discussed above corresponds to the special case $r = 0$. This generalised equilibrium condition can be related [445,452] to the KMS condition for equilibrium as discussed above.

In general, the functional $f[\psi]$ will not obey the conditions noted above, and so behaviour will not be equivalent to an equilibrium ensemble. However, in many cases (including the WIDBG model discussed above), a renormalisation group analysis of the action shows [242, 445, 446] that at long wavelengths, one reaches an effective action for which $\mathrm{Re}[f[\psi]] = r \mathrm{Im}[f[\psi]]$; i.e. this symmetry is an emergent property of the long wavelength theory.

One may observe that even without a full renormalisation analysis, the semiclassical limit of the WIDBG action, Eq. (136), takes the expected form of Eq. (137) with the noise strength $\Gamma = \kappa_- + \kappa_+$. As such, this yields a stochastic complex Gross-Pitaevskii equation, with the form:

$$i \frac{d}{dt} \psi_{cl} = \left( -\frac{\nabla^2}{2m} + \frac{U}{2} |\psi_{cl}|^2 + \frac{i}{2} \left[ \kappa_+ - \kappa_- - \frac{\kappa_-^{(2)}}{2} |\psi_{cl}|^2 \right] \right) \psi_{cl} + \xi, \qquad (142)$$

with the appropriate white-noise spectrum for $\xi$, with strength $\kappa_- + \kappa_+$.

## 5 Dissipative state engineering

After having discussed the general theoretical and experimental frameworks for the study of open quantum many-body systems, we now focus more specifically on one of its applications: the idea of environment or reservoir engineering. As already seen in several examples in these lecture notes, the steady state of a many-body system is generically mixed. Nevertheless, as we will extensively discuss in this section, under certain conditions, it is possible to use dissipation in order to engineer non-trivial, pure, many-body states in the stationary regime. This however is possible only through a careful tailoring of the coupling of the system to the environment. Originally proposed and successfully applied in the field of quantum optics and quantum information (see [453] for a recent review), where the focus was on single- or few-particle properties, be they atoms, molecules or photons, this idea is by now very well explored. The conceptual novelty that we want to discuss in this section is the recent extension to the realm of many-body quantum systems, which poses a set of novel problems, such as that of locality or of quantifying the typical timescale necessary to relax to the stationary state in the thermodynamic limit. This extension has been pioneered by two research articles appeared
in the late 2000's [8, 9], which outlined the potential for studies of dissipation beyond single-particle effects, proposing the development of universal quantum computers in a truly multi-qubit setting solely based on dissipation.

After a brief review on optical pumping in Sec. 5.1, we present in Section 5.2 the theoretical aspects of dark-state engineering [8, 9, 454] and focus then on the dissipative preparation of Bose-Einstein condensates [454] (Section 5.3.1) and of topological states hosting Majorana zero modes [455–458] (Section 5.3.2). We then discuss two examples that are of particular experimental relevance: the creation of entangled states using two-body losses [459–461], in Section 5.3.3, and the preparation of correlated many-body states of photons, such as Mott insulators [462], in Section 5.3.4. We continue with a discussion in Sec. 5.4 on the use of dark-state engineering for quantum-information purposes [9] and in Sec. 5.5 on the use of measurements as a manipulation technique, the so-called quantum many-body steering [463]. Some space will finally be devoted in Sec. 5.6 to the discussion of their stability against perturbations [464, 465].

In this section we focus on only some aspects of dissipative state engineering. There are however many other topics that could well fit into it and that will not be discussed for the sake of brevity. For instance, an enormous literature has focused on the problem of entanglement generation assisted by dissipation, even with many-body quantum systems; the topic has been recently reviewed in Ref. [453] and will not be discussed here. In the context of many-body fermionic fluids, the possibility of creating η-pairing superconductivity has recently attracted a lot of attention [454, 466, 467], whereas other works have focused on the dissipative realisation of matrix-product-states [391, 454, 468, 469]. We refer the interested reader to the primary literature for those topics.

## 5.1 From optical pumping to dark states

In the context of quantum physics, the study of states that cannot absorb nor emit photons—so-called *dark states*—dates back to the early studies on three-level systems, where laser light was used to couple short- and long-lived atomic energy levels (Lambda systems) [184, 470, 471]. The existence of such states plays a crucial role in the appearance of dark absorption lines, and has eventually led to the study of phenomena such as *electromagnetic induced transparency*, *coherent population trapping*, and *slow light* [184, 472]. Before discussing dark states in the many-body context, it is thus useful to start with a discussion of a simple single-atom manipulation via techniques that leverage both closed- and open-system dynamics. Two paradigmatic examples can be mentioned: *sideband cooling* [473] and *optical pumping* [474]. Since this latter phenomenon is particularly simple, we briefly detail it here below.

The goal of optical pumping, in its simplest form, is to prepare a cloud of atoms in some chosen internal hyperfine state. It relies on the angular-momentum conservation of the light-matter interaction, but also crucially on the possibility of manipulating the atom using spontaneous emission, which is an open-quantum-system effect. Optical pumping is best explained considering the simple situation of a transition connecting a ground state with angular momentum $J_g = 1/2$ to an excited state with angular momentum $J_e = 1/2$ (see the sketch in Fig. 20). If such a transition is driven with coherent light that is right-hand circularly polarised, the only coupling that is possible is that between the states with $m_g = \downarrow$ and $m_e = \uparrow$ (here $m_g$ and $m_e$ are the magnetic quantum numbers of the ground and excited states, respectively; we will use $\downarrow, \uparrow$ to label the corresponding quantum numbers $\pm 1/2$). The excited state is short lived and decays via spontaneous emission to both ground-state levels, $m_g = \downarrow, \uparrow$. Once the laser light illuminates the atomic vapour, a part of the population decays to the $m_g = \uparrow$ and is trapped: only undesired effects (not shown in the figure) can modify the state of those atoms. The atoms decaying to $m_g = \downarrow$, instead, are re-pumped to the excited state via the coherent field and the entire process is restarted. This eventually brings all the atomic population to the
<figure><img src="image_23.png" /><figcaption>Figure 20: Sketch of the simple optical-pumping process described in the text. A ground state ($J_g = 1/2$) is connected to a short-lived excited state ($J_e = 1/2$) with a coherent field with right-handed circular polarisation (solid arrow). The excited state $m_e = \uparrow$ can decay to both ground states $m_g = \uparrow, \downarrow$ (dashed arrows). The state $m_g = \uparrow$ is a dark state of this dynamics.</figcaption></figure>

state $m_g = \uparrow$, which is a dark state of this dynamics. Optical pumping of this kind is routinely performed in many cold-atom experiments.

Mathematically, we can describe the dynamics with the following single-atom Lindblad master equation:

$$ \frac{d\rho}{dt} = -i[\hat{H}, \rho] + \sum_{\sigma=\downarrow,\uparrow} \sum_{\sigma'=\downarrow,\uparrow} \left( \hat{L}_{\sigma,\sigma'} \rho \hat{L}_{\sigma,\sigma'}^{\dagger} - \frac{1}{2} \{ \hat{L}_{\sigma,\sigma'}^{\dagger} \hat{L}_{\sigma,\sigma'}, \rho \} \right), \quad (143) $$

with Hamiltonian

$$ \hat{H} = \sum_{v=e,g} \sum_{\sigma=\downarrow,\uparrow} E_v |m_v = \sigma\rangle\langle m_v = \sigma| + (\Omega e^{-i(E_e-E_g)t} |m_e = \uparrow\rangle\langle m_g = \downarrow| + \text{H.c.}), \quad (144) $$

and jump operators

$$ L_{\sigma, \sigma'} = \sqrt{\gamma_{\sigma, \sigma'}} |m_g = \sigma\rangle\langle m_e = \sigma'|. \quad (145) $$

The Hamiltonian describes the four energy levels and the coherent field that has a right-handed circular polarisation and couples selectively between only two of them. Four decay channels are however possible, from both excited states $\sigma'$ to both ground states $\sigma$, and are described by four different jump operators. It is possible to show that the pure state $|m_g = \uparrow\rangle$ is a stationary state of the dynamics, as its time-derivative equals zero. This follows from the fact that the density matrix $|m_g = \uparrow\rangle\langle m_g = \uparrow|$ commutes with the Hamiltonian and is annihilated by all four jump operators. In the next section we will see how this kind of reasoning can help us in the study of many-body problems.

## 5.2 Many-body dark states

The idea of generalising dark-state physics to the dissipative many-body scenario is more recent [8, 9, 475, 476] and while extensive theoretical work has been done, it has not yet been fully investigated experimentally. The starting point is a quantum many-body system coupled to an environment, whose dynamics is described by a Lindblad master equation, as introduced above. In the following we first discuss the nature of the steady (i.e. dark) states, and then discuss the asymptotic approach to such states.
### 5.2.1 Nature of the steady states

In the discussion presented in these lecture notes, the physics of dark states appears when we have one pure state $|\Psi\rangle$ such that (i) it is an eigenstate of the Hamiltonian, so that $[\hat{H}, |\Psi\rangle\langle\Psi|] = 0$, and (ii) it is in the kernel of all jump operators, $\hat{L}_\mu |\Psi\rangle = 0$. If conditions (i) and (ii) are satisfied, it is not difficult to show that $|\Psi\rangle\langle\Psi|$ is a stationary state of the dynamics:

$$ \frac{d}{dt} |\Psi\rangle\langle\Psi| = -i[\hat{H}, |\Psi\rangle\langle\Psi|] + \sum_{\mu} \hat{L}_{\mu} |\Psi\rangle\langle\Psi| \hat{L}_{\mu}^{\dagger} - \frac{1}{2} \{\hat{L}_{\mu}^{\dagger} \hat{L}_{\mu}, |\Psi\rangle\langle\Psi|\} = 0. \quad (146) $$

In such cases, $|\Psi\rangle$ is commonly dubbed “dark”. This is exactly what happens in the master equation that we studied above in Sec. 5.1 on optical pumping.

This construction can be easily generalised to dark subspaces in the case of many linearly independent dark states but will not be further considered here.

The most general conditions for the appearance of dark states are broader than those defined above. It has been formulated in Ref. [454]. A pure dark state $|\Psi\rangle$ exists if and only if: (i) it is an eigenstate of the effective non-Hermitian Hamiltonian associated to the problem, $\hat{H}_{\text{nH}} = \hat{H} - \frac{i}{2} \sum_{\mu} \hat{L}_{\mu}^{\dagger} \hat{L}_{\mu}$, with eigenvalue $\Lambda \in \mathbb{C}$, (ii) it is an eigenstate of all quantum-jump operators, $\hat{L}_{\mu}$ with eigenvalue $\lambda_{\mu}$, and (iii) the following relation holds: $\text{Im} \Lambda = -\frac{1}{2} \sum_{\mu} |\lambda_{\mu}|^2$. We refer to the original article, Ref. [454], for the proof of this. The more-restricted statement presented above corresponds to the case $\lambda_{\mu} = 0$ and $\Lambda \in \mathbb{R}$.

The class of restricted dark states, i.e. $\lambda_{\mu} = 0, \Lambda \in \mathbb{R}$ is clearly less general, yet, it is of interest as it leads to various non-trivial many-body dark states. Let us first observe that any such dark state satisfies the identities $\langle\Psi|\hat{L}_{\mu}^{\dagger} \hat{L}_{\mu}|\Psi\rangle = 0$ for all $\mu$; thus, since every operator $\hat{L}_{\mu}^{\dagger} \hat{L}_{\mu}$ is positive semi-definite, $|\Psi\rangle$ must minimise $\langle\Psi|\hat{L}_{\mu}^{\dagger} \hat{L}_{\mu}|\Psi\rangle$. Consequently, $|\Psi\rangle$ is a zero-energy ground state of the parent Hamiltonian:

$$ \hat{H}_p = \sum_{\mu} \hat{L}_{\mu}^{\dagger} \hat{L}_{\mu}. \quad (147) $$

This result is rather important because it gives a positive characterisation of the pure dark states corresponding to this situation. Indeed, Hamiltonians of the form (147) featuring a zero-energy ground state are well-known in the condensed-matter literature under the name of *frustration-free Hamiltonians*; they are the sum of non-negative operators that in general do not commute but that admit a common ground space. Notable examples include the toric-code Hamiltonian that plays a major role in the context of topological order [477] and the Rokhsar-Kivelson Hamiltonian for valence-bond states [478].

The concept of dark state has recently experienced a vigorous revival within the framework of dissipative quantum state preparation. If $|\Psi\rangle$ is the unique dark state of the dissipative dynamics and for some reasons one is interested in preparing it (for instance, it could be an entangled state with some relevance for quantum metrology), then the experimentalist may prepare the state via an accurate engineering of an appropriate environment, that will then create the state through dissipation [8, 9, 476]. This is what experimentalists working in quantum optics have long done at the level of single-particle systems with techniques like that of optical pumping, and that is currently being extended to the many-body realm. Clearly, interesting dark states do not appear in all quantum setups, but the set of systems where this takes place is wider than one could imagine at first sight. Before discussing explicit examples, however, it is important to stress a fundamental requirement that should be taken into account when dealing with many-body quantum systems: locality. A well-defined and physical many-body Lindbladian should have a local Hamiltonian and local jump operators, namely they should be the sum of operators with a support whose extension does not scale with the system size in
the thermodynamic limit. Relaxing this condition would give rise to an uncountable number of forms of unphysical dissipative dynamics that could be used to engineer literally any quantum state. Even with this constraint of locality, the list of theoretical studies where dissipative quantum state preparation with a physical Lindbladian has been discussed is now rather long.

### 5.2.2 Asymptotic decay toward the dark states

Whereas statements on the stationary states can be obtained with relative ease, the characterisation of the dynamics, and in particular of the long-time dynamics, is more difficult. The scaling with the system size of the asymptotic decay rate, which is the gap of the Lindbladian operator, has attracted a lot of attention because it determines the time scale that is necessary for the system to relax to the stationary state that is being engineered. In this section we will use the length of the system, $L$ to denote system size—for lattice based systems, one may note that $N_s \propto L^d$. When the asymptotic decay rate does not scale to zero with increasing system size, as in the example discussed in Sec. 5.3.2, the dissipative engineering is particularly efficient. In other situations, the asymptotic decay rate decreases polynomially with $L$, the size of the system; this is a less advantageous situations which can, nonetheless, be of interest.

There is no general recipe for the analytical calculation of the asymptotic decay rate in these setups. However, in several of the examples that will be presented in the remainder of the section we will detail specific techniques that have been designed for particular problems. Ideally, it would be desirable to link the asymptotic decay rate of the dissipative dynamics to the spectral properties of the parent Hamiltonian $\hat{H}_p$, which is a Hermitian operator, and for which, typically, more results are known. (Note we are only proposing here to use the parent Hamiltonian only as a simplified mathematical tool for establishing theorems, not as a fundamental description of the open system dynamics.) As an example of the possibilities of this approach, we remind here a result proven in Ref. [458] in the context of Lindbladians with no Hamiltonian. According to this work, when the parent Hamiltonian, Eq. (147), is gapless then the Lindbladian is also gapless. In particular, if $\hat{H}_p$ has a gap closing as $L^{-\beta}$, where $L$ is the size of the system, then the asymptotic decay rate of the Lindbladian closes at least as $L^{-\beta}$. The result is easily obtained by considering the first excited state of the parent Hamiltonian, dubbed $|\Psi_{ex}\rangle$, with energy $E_{ex}$. The linearly independent operators $|\Psi_{ex}\rangle\langle\Psi|$ and $|\Psi\rangle\langle\Psi_{ex}|$ are eigenoperators of the Lindbladian with eigenvalue $-E_{ex}/2$.

To conclude, we briefly comment on the accuracy of the many-body state by dissipative state engineering; as expected, waiting for longer times is beneficial to the quantum state preparation. If we define $p_0(t) = \text{Tr}[\hat{P}_0 \hat{\rho}(t)]$, where $\hat{P}_0 = |\Psi\rangle\langle\Psi|$ is the projector onto the ground space of the parent Hamiltonian $\hat{H}_p$, then the following monotonicity property holds: $\frac{d}{dt}p_0(t) \ge 0$. Indeed, with an explicit calculation, $\frac{d}{dt}p_0(t) = \sum_\mu \text{Tr}[L_\mu^\dagger \hat{P}_0 L_\mu \hat{\rho}(t)]$, which is the sum of the expectation values of non-negative operators, and thus of non-negative numbers.

## 5.3 Examples

Having exposed some general considerations on dissipative state preparations, we now move to the presentation of several concrete examples. We consider dissipative quantum systems that leverage the notion of dark states, based on different experimental setups such as trapped ions [66, 67, 479], Rydberg atoms [134], ultracold atoms [8, 480]. We note that this is not an exhaustive summary of the important literature that has appeared in recent years.

### 5.3.1 BEC condensate

As a first example, we present here the dissipative preparation of Bose-Einstein condensates with long-range order [8, 454]; the fermionic case is discussed in Ref. [480]. Let us consider a
<figure><img src="image_24.png" /><figcaption>Figure 21: An example of topological physics induced by dissipation. (Left) When a Hilbert space $\mathcal{H}$ is dissipatively constrained to a subspace $S$, the dynamics is governed by the projected Hamiltonian $\hat{H}_S = \hat{P}_S \hat{H} \hat{P}_S$, where $\hat{H}$ and $\hat{P}_S$ are the bare Hamiltonian and the projection operator onto $S$, respectively. (Middle) In the situation discussed here, the authors propose a dissipative mechanism that constraints the dynamics in the lower band of a lattice, that has a non-trivial Berry curvature $\Omega_{xy}$, plotted over the first Brillouin zone. (Right) As a final experimentally observable phenomenology associated to topological physics, it is shown that a wave packet undergoes transverse motion in response to a potential gradient $\mathbf{F}$. Adapted from Ref. [481] [Copyright (2017) by the American Physical Society].</figcaption></figure>

lattice populated by bosons, described by the operators $\hat{b}_i^{(\dagger)}$ that satisfy canonical commutation relations, where $i$ labels the lattice site. For simplicity, we set the Hamiltonian to zero and consider only the dissipative dynamics induced by the following jump operators:

$$ \hat{L}_{i,j} = \sqrt{\gamma} (\hat{b}_i^{\dagger} + \hat{b}_j^{\dagger}) (\hat{b}_i - \hat{b}_j), \qquad (148) $$

acting on a pair of adjacent sites, where $\gamma$ is a decoherence rate. Intuitively, this dissipative operator eliminates from the lattice a boson that is antisymmetric component of the wave-function on two neighbouring lattice sites and re-introduces it in a symmetric mode. Let us now consider the $N$-particle condensate state $|BEC_N\rangle = \frac{1}{\sqrt{N!}} (\hat{b}_{k=0}^{\dagger})^N |v\rangle$ where $\hat{b}_{k=0} = \frac{1}{\sqrt{N_s}} \sum_j \hat{b}_j$ is the $k=0$ bosonic mode and $|v\rangle$ the vacuum. The state $|BEC_N\rangle$ is the $N$-particle ground state of a problem of free bosons, and represents an exact idealised Bose-Einstein condensate where all particles have condensed in one bosonic mode. In order to show that $|BEC_N\rangle$ is a dark state of the dissipative dynamics induced by the jump operators in Eq. (148), it is sufficient to show that $(\hat{b}_i - \hat{b}_j)|BEC_N\rangle = 0$ for any pair of neighbouring sites. The thesis follows from the fact that the commutator $[\hat{b}_i, \hat{b}_{k=0}^{\dagger}] = \frac{1}{\sqrt{N_s}}$ does not depend on $i$ and can be proven by induction over $N$.

### 5.3.2 Topological states and Majorana dark states

Here we want to present the possibility of creating topological states of matter using dissipation, for which there are many examples. In Ref. [481], from which Fig. 21 is extracted, strong dissipation is exploited to confine the dynamics to a lattice band with non-trivial Berry curvature [481]. Related works have extensively discussed the possibility of using non-Hermitian dynamics, losses or other decoherence sources to engineer p-wave superfluids [482], integer or fractional quantum Hall states [483, 484], as well as integer [485] or fractional Chern insulators [486–489]. The use of Rydberg atoms has been instead proposed for the dissipative engineering of Kitaev’s toric code [134]; related to this subject is also the study of open quantum systems where strong dissipation is employed to create lattice gauge theories with cold atoms [490]. Recently, there has been very intriguing progress in realising Laughlin states composed of few photons, which is a setup that is intrinsically dissipative [146, 491].
In the following we will discuss in some details the case of topological superconductors that host Majorana zero-energy modes [492] (although, as we will see, in this out-of-equilibrium scenario, it will be more appropriate to speak of Majorana dark modes or, equivalently, of Majorana zero-damping modes). We discuss in particular the results of Ref. [455] and use them to illustrate how the notion of dark state can be used to guide the dissipative state preparation of Majorana zero modes.

The simplest model displaying zero-energy unpaired Majorana modes is the Kitaev model for a one-dimensional p-wave superconductor. We focus on one specific point of the phase diagram where its solution is particularly straightforward [492]. We introduce the fermionic operators $\hat{c}_j^{(\dagger)}$ that satisfy canonical anticommutation relations and describe the annihilation and creation of a spinless fermion at site j. The Hamiltonian reads

$$ \hat{H}_K = -t_H \sum_j [\hat{c}_j^\dagger \hat{c}_{j+1} + \hat{c}_j \hat{c}_{j+1} + \text{H.c.}], \quad t_H > 0. \quad (149) $$

The model can be solved with the Bogoliubov-de Gennes transformation, and, when considered on a chain of length L with open boundaries, it takes the form:

$$ \hat{H}_K = E_0 + \frac{t_H}{2} \sum_{j=1}^{L-1} \hat{l}_j^\dagger \hat{l}_j, \quad \text{with} \quad \hat{l}_j = \hat{C}_j^\dagger + \hat{A}_j, \quad \hat{C}_j^\dagger = \hat{c}_j^\dagger + \hat{c}_{j+1}^\dagger, \quad \hat{A}_j = \hat{c}_j - \hat{c}_{j+1}. \quad (150) $$

The modes $\hat{l}_j$ are fermionic operators that obey the canonical anticommutation relations, $\{\hat{l}_i, \hat{l}_j\} = 0$ and $\{\hat{l}_i, \hat{l}_j^\dagger\} = \delta_{ij}$, and that can be used to obtain the entire spectrum of the model. The ground state has energy $E_0$ and is two-fold degenerate: there are two linearly independent states $|\psi_e\rangle$ and $|\psi_o\rangle$ which satisfy $\hat{l}_j|\psi_\sigma\rangle = 0$ for all j. The quantum number distinguishing the two states is the parity of the number of fermions, $\hat{P} = (-1)^{\sum \hat{c}_j^\dagger \hat{c}_j}$, which is a symmetry of the model (indeed, the subscripts e and o stand for even and odd). Both states $|\psi_\sigma\rangle$ are p-wave superconductors, as can be explicitly proven by computing the expectation value of the corresponding order parameter: $\lim_{L\to\infty} \langle\psi_\sigma|\hat{c}_j\hat{c}_{j+1}|\psi_\sigma\rangle = 1/4$. The rest of the spectrum is obtained by subsequent application of the $\hat{l}_j^\dagger$ operator, knowing that in this precise and simple problem each excitation has an energy cost of $t_H/2$.

It is not difficult to design a master equation which features $|\psi_e\rangle$ and $|\psi_o\rangle$ as steady states and has local jump operators. This is possible because, as we have just seen, the two states lie in the kernel of the entire set of $\hat{l}_j$ operators. Introducing a Markovian dynamics with jump operators $\hat{L}_j = \sqrt{\gamma}\hat{l}_j$ ensures that the states $|\psi_\sigma\rangle$ are stationary and that in the long-time limit the system is brought into a subspace described in terms of p-wave superconducting states. Additionally, the parent Hamiltonian of this Markov process coincides with $\hat{H}_K$ in Eq. (150) up to an additive constant. The obtained dynamics satisfies an important property, namely locality: the jump operators $\hat{l}_j$ only act on two neighbouring fermionic modes.

In such a model, there are two unpaired Majorana modes that are not affected by the dissipative dynamics:

$$ \hat{\gamma}_L = \hat{c}_1 + \hat{c}_1^\dagger, \quad \hat{\gamma}_R = i(\hat{c}_L - \hat{c}_L^\dagger). \quad (151) $$

In the Hamiltonian context, these boundary operators are typically called Majorana zero modes, because they commute with the Hamiltonian and force the entire spectrum of the Hamiltonian to be two-fold degenerate. This can be shown by creating a complex fermion $\hat{d} = \hat{\gamma}_L + i\hat{\gamma}_R$ that squares to zero ($\hat{d}^2 = 0$), satisfies canonical anticommutation relations ($\{\hat{d}, \hat{d}^\dagger\} = 1$), and that anticomutes with all the Bogoliubov modes $\hat{l}_j$ and $\hat{l}_j^\dagger$. Since it commutes with the Hamiltonian, it is said to be a zero-energy operator. This is the operator that differentiates the two ground states: $\hat{d}|\psi_e\rangle = 0$ and $\hat{d}^\dagger|\psi_e\rangle = |\psi_o\rangle$. More generally, the application of $\hat{d}^\dagger$ or $\hat{d}$ on one of the eigenstates listed above allows for the construction for
<figure><img src="image_25.png" /><figcaption>Figure 22: Proposed implementation scheme for the Lindblad master equation featuring Majorana boundary modes. The quantum wire is represented by the lower sites of an optical superlattice for spin polarised atomic fermions. They are coherently coupled to auxiliary upper sites by lasers with Rabi frequencies ±Ω, alternating from site to site. Dissipation results from spontaneous Bogoliubov phonon emission via coupling of the system to a BEC reservoir (light grey). As shown in the original article, this setting reduces to the jump operators $\hat{L}_j = \hat{C}_j^\dagger \hat{A}_j$ at late times, see Eq. (152). Inspired from Ref. [455].</figcaption></figure>

another orthogonal state that has the same energy. It is thus responsible for the exact two-fold degeneracy of the entire spectrum of the Hamiltonian that we claimed, and of course also of the topologically protected degeneracy of the ground state.

In the dissipative context, the interpretation of the $\hat{\gamma}_{L,R}$ operators is slightly different, because now they represent operators that are not affected by the dissipation. As has been anticipated, it is more correct to speak of dark modes or of zero damping modes. Interestingly, as was also shown in the Hamiltonian context, Majorana zero damping modes also satisfy a non-trivial braiding statistics [455].

We may also mention that the jump operators $\hat{l}_j$ are linear in the fermionic operators, which yields a quadratic Lindbladian that can be exactly diagonalised, as we discussed in Sec. 4.1. A simple analysis shows that the asymptotic decay rate of the problem is $\gamma$ and that it does not scale with the size of the system.

The approach discussed so far suffers from the difficulty of an experimental implementation of a jump operator $\sqrt{\gamma} \hat{l}_j$ that does not conserve the number of particles. However, as was already noted in the original reference [455], another approach exists based on the experimental implementation of a number-conserving jump operator:

$$ \hat{L}_j = \sqrt{\gamma} \hat{C}_j^\dagger \hat{A}_j. \qquad (152) $$

Similarly to the jump operators in Eq. (148), these operators can be interpreted as the annihilation of an antisymmetric state and creation of a symmetric one, while the number of fermions is conserved by the ensuing dissipative dynamics.

It is possible to characterise the stationary states of the dynamics of this model [493]. We first observe that $\hat{l}_j |\psi_\sigma\rangle = 0$ implies $\hat{C}_j^\dagger |\psi_\sigma\rangle = -\hat{A}_j |\psi_\sigma\rangle$, so that:

$$ \hat{L}_j |\psi_\sigma\rangle = \hat{C}_j^\dagger \hat{A}_j |\psi_\sigma\rangle = -\hat{C}_j^\dagger \hat{C}_j^\dagger |\psi_\sigma\rangle = 0. \qquad (153) $$

Thus, the $|\psi_\sigma\rangle$ are also steady states of this new dynamics. In fact, we can say more: any projection of $|\psi_\sigma\rangle$ that has a well-defined number of particles is stationary:

$$ |\psi_N\rangle = \frac{\hat{\Pi}_N |\psi_\sigma\rangle}{\sqrt{\langle\psi_\sigma|\hat{\Pi}_N|\psi_\sigma\rangle}}, \qquad (154) $$
where $\hat{\Pi}_N$ is the projector onto the subspace of the global Hilbert (Fock) space with $N$ fermions. In Ref. [458] the authors show that there is only one dark state $|\psi_N\rangle$ once the value of $N$ is fixed. They further conclude that the dynamics induced by the jump operators in (152) conserves the number of particles and drives the system into a quantum state with the properties of a p-wave superconductor, that is, in the thermodynamic limit, $|\psi_e\rangle$ and $\hat{\Pi}_N |\psi_e\rangle$ have the same bulk properties. However, since the steady states of the system for open boundary conditions are unique, they do not display any topological edge states. Figure 22 illustrates some ideas about the experimental implementation of this scheme.

Finally, it is interesting to inspect the parent Hamiltonian associated to this alternate dissipative process:

$$ \hat{H}_p = \sum_j \hat{L}_j^\dagger \hat{L}_j = \gamma \sum_j \hat{A}_j^\dagger \hat{C}_j \hat{C}_j^\dagger \hat{A}_j = -2\gamma \sum_j [\hat{c}_{j+1}^\dagger \hat{c}_j + \hat{c}_j^\dagger \hat{c}_{j+1} + 2\hat{n}_j \hat{n}_{j+1} - \hat{n}_j - \hat{n}_{j+1}] . \quad (155) $$

This parent Hamiltonian can be mapped—via a Jordan-Wigner transformation—to the ferromagnetic spin-1/2 Heisenberg Hamiltonian, which is gapless for any $N$ scaling linearly with $L$. We can thus employ results discussed above to conclude that the dissipative process is characterised by an asymptotic decay rate that scales at least as $L^{-2}$, the scaling of the gap of the parent Hamiltonian.

Following these foundational results reviewed here, the literature on this subject has been particularly rich. The two-dimensional case was also studied, and the dissipative creation of a vortex trapping a Majorana mode is discussed in Ref. [494]. A more general study of topological phases induced by dissipation in fermionic systems is presented in Ref. [456]. The number-conserving case is discussed in Ref. [458]. These results have opened the way to a more general investigation of the notion of topology in open quantum systems and for density matrices [495–512]. Finally, we mention that the topic of topological dark states is also interwined with the long-standing quest of engineering fractional quantum Hall states in photonic setups [146, 457, 513–524].

### 5.3.3 Spin-Dicke states of fermionic gases

In this section we will focus on fermionic gases, where the interplay of statistics and spin produces intriguing steady states, as first pointed out in Ref. [459]. Our presentation follows in several points the more recent Ref. [525]. In contrast to the setup in the previous section that required careful engineering, it is interesting to discuss here the dark states associated to one of the most natural open-system processes in cold-atom experiments: atom losses.

We consider a $d$-dimensional hyper-cubic optical lattice in the single-band approximation populated by a fermionic spin-1/2 gas. We introduce the annihilation and creation operators $\hat{c}_{j,\sigma}^{(\dagger)}$ associated to the site $j$ and to the spin $\sigma$ that satisfy the canonical anticommutation relations. The Hamiltonian of the system is the standard Hubbard model:

$$ \hat{H} = -t_H \sum_{\langle i,j \rangle} \sum_{\sigma=\uparrow,\downarrow} [\hat{c}_{i,\sigma}^{\dagger} \hat{c}_{j,\sigma} + \text{H.c.}] + U \sum_j \hat{n}_{j,\uparrow} \hat{n}_{j,\downarrow}. \quad (34) $$

We will focus on the experimentally relevant case of local two-body losses with rate $\gamma$, described by the following set of jump operators:

$$ \hat{L}_j = \sqrt{\gamma} \hat{c}_{j,\uparrow} \hat{c}_{j,\downarrow}. \quad (35) $$

Fermionic statistics here poses a constraint on the single-site loss process, that it can only remove two particles with opposite spins because it is not possible to have two particles with the same spin sitting at the same site. It follows immediately that the $z$ component of magnetisation of the gas is conserved during the dynamics. In fact, the model is SU(2) invariant and we
<figure><img src="image_26.png" /><figcaption>Figure 23: (Top) Sketch of the experiment on two-body losses. A one-dimensional array of double-well potentials is initially filled with one fermionic particle per site; six spin components are possible but in a double-well only two of them appear, so that the description of a single double-well can be carried out in terms of an effective spin-1/2 setup, that can be in a spin singlet or triplet state. Singlet states are characterised by a spatially symmetric wavefunctions, this allows the two particles to sit on the same site and to be lost. Triplet states, characterised by an antisymmetric wavefunction, cannot have a double site occupation, and are never lost. In the experiment, on-site two-body losses are induced with a photoassociation process. (Bottom) Dynamics of (b, left) the singlet fraction and of (c, right) the triplet fraction in all double-wells as a function of the time to which they are exposed to photoassociating light inducing losses. Whereas the singlet fraction decreases exponentially with time, the triplet fraction remains constant. The three set of data (i), (ii) and (iii) correspond to three different initial conditions Adapted from Ref. [461] [Copyright (2023) by the American Physical Society.].</figcaption></figure>

could have expressed the jump operator in Eq. (35) in a basis with a different spin-quantisation axis, always with the result that the loss process removes two fermions with opposite spin. We can thus conclude that all components of magnetisation of the gas are a strong symmetry of the dynamics.

More remarkably, we can now observe that the two particles lost by the dissipative process in Eq. (35) are in a spin singlet. This follows easily from the fact that they share the same orbital wavefunction. A loss process that only annihilates spin singlets conserves the total spin of the system: the expectation value of $\hat{S}^2$ is conserved during the dynamics. Using the classification introduced in Section 3.1.1, it follows from what we just wrote that also $\hat{S}^2$ is a strong symmetry of the dissipative dynamics. We can visualise this by considering the dissipative dynamics of an initial state with only two fermions, that will generically be a linear superposition of a spin-singlet $S=0$ state and of a spin-triplet $S=1$. For the $S=1$ state the orbital wavefunction must be antisymmetric with respect to the exchange of particles, and thus the two fermions can never occupy the same site: the $S=1$ spin channel is thus not lossy. Conversely, the $S=0$ spin channel can have particles on the same site and can dissipate them. Thus, as the average number of particles decreases, the spin $\langle \hat{S}^2 \rangle$ of the system remains constant. This simple observation has recently been the object of an experimental measurement [461]. As
shown more clearly in Fig. 23, the authors have prepared a fermionic gas into several double wells, each filled with exactly two atoms, and have observed a decrease with time of the number of atoms in the singlet channel, while the number of atoms in the triplet channel was held constant.

Such spin-selective decay properties can be generalised to an arbitrary number N of initial particles. In this case the relation between spin and dynamics is more complicated, but we can easily draw some general conclusions from the simple observation that during the loss dynamics the number of particles decreases while the spin of the gas remains constant. In fact, due to the spin conservation the number of particles of the gas cannot decrease arbitrarily, as a gas with N particles can have a spin quantum number S that is upper bounded by N/2 and thus:

$$ \langle \hat{S}^2 \rangle \le \left\langle \frac{\hat{N}}{2} \left( \frac{\hat{N}}{2} + 1 \right) \right\rangle. \qquad (156) $$

Since $\langle \hat{S}^2 \rangle$ is conserved, this equation places a lower bound on the number of particles. Stationary states thus are expected to be incoherent mixtures of states with maximal spin, whose wavefunction has a spin part that is a Dicke state. It is known that Dicke states have a fully symmetric spin wavefunction: as a consequence, the orbital wavefunction is completely anti-symmetric and cannot have double occupancies.

We can give a more rigorous description of these dark states of the loss dynamics. We begin by constructing the set of the N-particle eigenstates of the Hubbard Hamiltonian in Eq. (34) that are fully polarised. In an hypercubic plane with periodic boundary conditions, the eigenmodes of the kinetic energy part of the Hamiltonian (34) are plane waves with wavevector k. We thus introduce the operator:

$$ \hat{c}_{k,\sigma}^{\dagger} = \frac{1}{\sqrt{N_s}} \sum_{j} e^{i\mathbf{k}\cdot\mathbf{r}_j} \hat{c}_{j,\sigma}^{\dagger}, \qquad (157) $$

and the states:

$$ \hat{c}_{k_1, \uparrow}^{\dagger} \hat{c}_{k_2, \uparrow}^{\dagger} \cdots \hat{c}_{k_N, \uparrow}^{\dagger} |0\rangle, \qquad (158) $$

where $|0\rangle$ is the empty state, $\hat{c}_{k,\sigma} |0\rangle = 0$. These states are the dark states of the loss dynamics, as can be proven by showing that they are eigenstates (with real eigenvalue) of the non-Hermitian Hamiltonian associated with the Hamiltonian in Eq. (34) and the jump operators in Eq. (35):

$$ \hat{H}_{\text{nH}} = \hat{H} - \frac{i}{2} \sum_j \hat{L}_j^\dagger \hat{L}_j = -t_H \sum_{\langle i,j \rangle} \sum_{\sigma=\uparrow,\downarrow} \left[ \hat{c}_{i,\sigma}^\dagger \hat{c}_{j,\sigma} + \text{H.c.} \right] + \left( U - i \frac{\gamma}{2} \right) \sum_j \hat{n}_{j,\uparrow} \hat{n}_{j,\downarrow}. \qquad (159) $$

This non-Hermitian Hamiltonian is a Hubbard model with complex interaction, where the interacting term has a conservative and non-conservative parts, and plays an important role in the physics of fermionic gases with two-body losses. The states in Eq. (158) have wavefunctions that factorise the orbital and spin parts. The orbital part is a Slater determinant constructed with the single-particle eigenmodes of the kinetic energy. The dark states that are not fully polarised in the z direction can be obtained by applying the spin-lowering operator

$$ \hat{S}^{-} = \sum_{\mathbf{k}} \hat{c}_{\mathbf{k}, \downarrow}^{\dagger} \hat{c}_{\mathbf{k}, \uparrow}, \qquad (160) $$

to the states in Eq. (158). The SU(2) invariance of the non-Hermitian Hamiltonian can then be invoked to argue that they are also dark states of the dynamics.

These results have been recently generalised to SU(3) gases subject to two-body losses, where the stationary states are discussed using the representation theory of the SU(3) symmetry group and are dubbed generalised Dicke states [460, 525].
The study of the asymptotic decay rates of the loss dynamics that we described so far are actually related to the spectral properties of the non-Hermitian Hubbard model in Eq. (159). This is a general property of Lindblad master equations with only losses that follows from the fact that the Lindbladian operator can be put in a block triangular form [526]. In a few words, this result can be understood in the following way. We first observe that loss operators are characterised by a weak symmetry, generalising that of the toy model introduced in Eq. (15), namely symmetry under the transformations $e^{i\theta \sum_{j,\sigma} \hat{n}_{j,\sigma}}$. As noted in Sec. 3.1.1, weak symmetries mean the Lindbladian has a block diagonal structure. If one labels blocks of the density matrix as $\hat{\rho}_{n_L, n_R}$, where $n_L, n_R$ correspond to eigenstates of the total-particle-number operator $\hat{N} = \sum_{j,\sigma} \hat{n}_{j,\sigma}$, the dynamics of the part of the density matrix with $n_L = n_R$ decouples from those parts where $n_R \neq n_L$. To reach the triangular form, a further requirement is needed: the dissipation must consist of only loss terms. This means that the master equation has the following block form:

$$ \frac{d}{dt} \hat{\rho}_{n,n} = \hat{\mathcal{L}}_{n}^{(0)}[\hat{\rho}_{n,n}] + \hat{\mathcal{L}}_{n}^{(+2)}[\hat{\rho}_{n+2,n+2}]. \qquad (161) $$

Physically this states that the time evolution of the block corresponding to $n$ particles depends only on itself, or on a source term from the block with $n + 2$ particles. This yields a block triangular form as noted above. The number conserving part, $\hat{\mathcal{L}}_{n}^{(0)}$ corresponds to the non-Hermitian Hamiltonian, while the jump operators lead to the coupling from the $n + 2$ block. Because $\hat{\mathcal{L}}_{n}^{(+2)}$ are just loss processes, they appear on the upper half of the triangle. This means that the spectrum of the Lindbladian is determined only by the diagonal blocks, i.e. the non-Hermitian Hamiltonian. Note that this behaviour would not arise if there were incoherent pumping to balance losses, even if the pumping respected the weak symmetry $\hat{S}$.

Using the above, one can show that for one-dimensional gases, the Lindbladian is gapless and the asymptotic decay rate scales as $1/L^3$ for states with total momentum close to $K \sim 0$ and with a finite density of excitations. We refer the reader to the specialised publications for further details [525, 527]. This result is derived analytically exploiting the fact that the Hamiltonian in (159) can be solved with the techniques of Bethe Ansatz (see also Refs. [319, 323, 528, 529] for other examples).

### 5.3.4 Dissipative preparation of Mott insulators of photons

A recent example of successful experimental dissipative state preparation has been the realisation of a Mott insulator of microwave photons using superconducting circuits [57, 462]. Although this case does not deal with an exact dark state of the dissipative dynamics, the result has many analogies with the spirit of the discussion presented out so far. Crucially, it relies on dissipation to achieve the state preparation and is of significant experimental relevance.

A major challenge in realising strongly correlated fluids of light is the fact that photons are not conserved particles whose density can be tuned by an external chemical potential. Starting from early results on lasing and polariton condensation (See Sec. 2.6) and the weakly interacting Bose-Einstein condensates of photons, first reported in Ref. [188], several theoretical proposals have been put forward to realise an effective chemical potential for strongly interacting photons using tailored dissipative reservoirs [530–533]. In Ref. [462] the authors realised a 1D chain of eight coupled superconducting qubits (transmons [534] in their implementation), realising a Bose-Hubbard Hamiltonian. Each site has intrinsic losses which would lead to an empty state. In order to stabilise a Mott insulating phase with exactly one photon on each site, the authors implemented two versions of a dissipative stabiliser, respectively involving one and two auxiliary transmons. We will describe the working principle of the stabiliser, in the scheme involving one auxiliary transmon for simplicity, first discussing the case in which the system consists of a single site coupled to the auxiliary transmon.
<figure><img src="image_27.png" /><figcaption>Figure 24: Scheme for the one-transmon dissipative stabiliser implemented in the experiment. (a) Level scheme showing a two-photon drive resonant with the $0 \leftrightarrow 2$ transition of the stabiliser site, s. A reservoir site r is tuned so that the $n_s = 2 \rightarrow n_s = 1$ transition of the stabiliser site is resonant with the $n_r = 0 \rightarrow n_r = 1$ transition. The reservoir site then quickly decays at rate $\kappa_r$. The overall effect is to stabilise the $n_s = 1, n_r = 0$ state, as highlighted in pink. (b) Integration of the stabiliser into a Bose-Hubbard lattice. Inspired by Ref. [462].</figcaption></figure>

We start by considering the situation where we want to prepare exactly one photon on a single site, assuming for simplicity that all the *n*-photon states have a long but not infinite lifetime. The solution that is routinely applied in isolated quantum systems cannot be applied: a direct transition via a $\pi$ pulse to the $n=1$ state cannot work because the state will decay at a certain unknown time, so that it is not obvious a priori when one should reapply the $\pi$ pulse. The envisaged solution is sketched in Fig. 24 and is based on coherent transitions from the $n=0$ state the $n=2$ state with a classical field of frequency $\omega_{02}/2$: the $n=2$ state is thus populated via a two-photon process. Thanks to the on-site Hubbard interaction $U$, the drive frequency is detuned from the transition from the $n=0$ to the $n=1$ state by $U/2$: $\omega_{02}/2 \neq \omega_{01}$ and hence the $n=1$ state remains unpopulated. The single site is close to an auxiliary and strongly dissipative transmon that has energy $\omega_{12}$: the two-photon state is thus resonant with the state where one photon is in the site to be stabilised and one photon in the auxiliary site. Once the coherent transition to this state has taken place, the auxiliary transmon dissipates its photon in the environment and the target site is left with exactly one photon. The interesting point of this process is that it is initiated by a classical driving field that is left switched on during the entire procedure. When the $n=1$ state is populated, the classical drive is off-resonance and cannot induce any transition in the system. However, as soon as the $n=1$ state eventually decays, the classical drive becomes resonant again and rapidly refills the one-photon state.

When considering the whole chain instead of a single site, the scheme is generalised stabilising only the first site of the chain, so that it acts as a narrow-band photon source for the rest of the chain. Photons from this site tunnel into the other sites until filling them all. When all the sites are filled with one photon, then the stabiliser will be unable to inject additional photons, as this would require an energy higher than the source energy, because of the incompressibility of the Mott phase. If instead a photon is lost from the chain, then the stabilised site pumps one photon in the chain. The experimental results are presented in Ref. [462] showing that within a few microseconds, the photon fluid has reached a steady state with near-unity occupancy.
A driving scheme in the same spirit of having an energy selective source of particles was proposed in [531], which is part of a general set of studies on reservoir engineering for many-body photonic states [530, 533, 535, 536] and on the theoretical characterisation of the dissipatively prepared Mott state [368, 537]. Other experiments involving strongly correlated photonic states have been performed [538, 539] and more experiments are to be expected; for instance, similar approaches are currently under discussion for the realisation of topological strongly correlated photonic fluids [146, 457, 513–524] (see also the previous section).

## 5.4 Quantum computation using dissipation

The notion of dark states is more versatile than what we have discussed so far and has applications that reach beyond the idea of dissipative state preparation. In a seminal article, F. Verstraete, M. M. Wolf and J. I. Cirac have discussed the computational power of the dissipative creation of many-body dark states, and have shown that the techniques explained above can be used to perform universal quantum computation [9]. The core of the reasoning proposed by the authors is rather simple and goes as follows. A set of N qubits is considered over which a set of unitary gates {$\hat{U}_t$}$_{t=1}^T$ needs to be applied; each of the unitary gates $U_t$ is assumed to act only on two neighbouring spins. The system is initialised in $|\psi_0\rangle = |0\rangle_1 \otimes |0\rangle_2 \cdots \otimes |0\rangle_N$, the state where all N qubits are in the state $|0\rangle$. The final outcome of the quantum computation is obtained by the ordered application of all gates, $|\psi_T\rangle = \hat{U}_T \cdots \hat{U}_2 \hat{U}_1 |\psi_0\rangle$. The main point of the authors is that $|\psi_T\rangle$ can be obtained as a dark state of a properly engineered Lindblad master equation.

The Lindbladian proposed to realise this desired state does not have any Hamiltonian part but only a set of jump operators that need to be carefully tailored to implement the gates $\hat{U}_t$. In order to define the jump operators, it is necessary to introduce an additional register, i.e. an additional quantum system, that in this case is defined on a T-dimensional Hilbert space with basis {$|t\rangle$}$_{t=0}^T$. The quantum jump operators are split into two sets and read:

$$ \hat{L}_i = |0\rangle\langle 1|_i \otimes |0\rangle\langle 0|_r, \qquad i = 1, 2, \dots, N, \qquad (162) $$

$$ \hat{L}_t = \hat{U}_t \otimes |t\rangle\langle t-1|_r + \hat{U}_t^\dagger \otimes |t-1\rangle\langle t|_r, \qquad t = 1, 2, \dots, T. \qquad (163) $$

The N operators $L_i$ act on the register and on the i-th qubit. The T operators $L_t$ act on the register and those qubits on which the gate $\hat{U}_t$ acts.

The stationary state of such a Lindbladian reads:

$$ \hat{\rho}_0 = \frac{1}{T+1} \sum_{t=0}^{T} |\psi_t\rangle\langle\psi_t| \otimes |t\rangle\langle t|_r, \quad \text{with} \quad |\psi_t\rangle = \hat{U}_t \cdots \hat{U}_2 \hat{U}_1 |\psi_0\rangle. \qquad (164) $$

In order to prove this, let us first consider the jump operators $\hat{L}_i$. With a simple calculation $\hat{L}_i \hat{\rho}_0 \hat{L}_i^\dagger = 0$; moreover $\hat{L}_i^\dagger \hat{L}_i = |1\rangle\langle 1|_i \otimes |0\rangle\langle 0|_r$, and thus {$\hat{L}_i^\dagger \hat{L}_i, \hat{\rho}_0$} = 0. This shows that the state $\hat{\rho}_0$ is stationary with respect to the action of the first set of jump operators. Let us now consider the second set of $\hat{L}_t$ operators. With a bit of algebra, one shows that:

$$ \begin{aligned} \hat{L}_t \hat{\rho}_0 \hat{L}_t^\dagger &= \hat{U}_t |\psi_{t-1}\rangle\langle\psi_{t-1}| \hat{U}_t^\dagger \otimes |t\rangle\langle t|_r + \hat{U}_t^\dagger |\psi_t\rangle\langle\psi_t| \hat{U}_t \otimes |t-1\rangle\langle t-1|_r \\ &= |\psi_t\rangle\langle\psi_t| \otimes |t\rangle\langle t|_r + |\psi_{t-1}\rangle\langle\psi_{t-1}| \otimes |t-1\rangle\langle t-1|_r. \end{aligned} \qquad (165) $$

After having computed that $\hat{L}_t^\dagger \hat{L}_t = I \otimes (|t\rangle\langle t| + |t-1\rangle\langle t-1|)$ it is not difficult to show that the anticommutator in the Lindblad master equation from this jump operator takes the same form as Eq. (165):

$$ \frac{1}{2} \{ \hat{L}_t^\dagger \hat{L}_t, \hat{\rho}_0 \} = |\psi_t\rangle\langle\psi_t| \otimes |t\rangle\langle t|_r + |\psi_{t-1}\rangle\langle\psi_{t-1}| \otimes |t-1\rangle\langle t-1|_r. \qquad (166) $$
It follows that the state $\hat{\rho}_0$ is stationary also with respect to the second set of jump operators.

A remarkable result is that the asymptotic decay rate for reaching the stationary state decays as $T^{-2}$ and thus the computation time diverges polynomially with the number of gates to be applied during the computation. The complexity of the computation is therefore hidden in the scaling of the number of gates $T$ with the number of qubits $N$. The uniqueness of the steady state can be proved by diagonalising the Lindbladian.

It is interesting to conclude this Section by making contact with earlier works on quantum information processing assisted by open-quantum-system techniques, as for instance the work by J.P. Paz and W.H. Zurek dating from 1998 where the idea of continuous quantum error correction is put forward [540]. They consider a quantum error-correction scheme where the system is repeatedly corrected and show that as the correction rate goes to infinity the dynamics of the system is actually described by a Lindblad master equation. Recent advances in the development of synthetic systems have revitalised these ideas, both from the viewpoint of error correction [541, 542] or of quantum memories [543]. Moving beyond the problem of safely storing quantum information with the help of dissipation, it was proposed in the early 2000's to perform quantum computation assisted by the presence of an environment [475, 544]. The idea differs from that of Ref. [9] in that the dark subspace is used to confine the entire quantum computation process, rather than engineering a dark state that contains the answer to the quantum computation.

## 5.5 Quantum many-body steering

Within this Section on dissipative state preparation we may also consider a different, although related, protocol based on quantum steering. The idea of steering is as old as quantum mechanics. Originally introduced by Schrödinger [545], its meaning evolved progressively to indicate questions related to the understanding of non-locality and the existence (or rather impossibility) of hidden-variable models. Quantum steering refers to the non-classical correlations that can be observed between the outcomes of measurements on entangled state and the properties of the state after the measurement, see Ref. [546] for a comprehensive review on this topic.

In a broader sense, but in line with the initial formulation, quantum steering protocols can be also related to quantum state preparation. The underlying idea is to control the dynamics of a quantum system by means of a sequence of measurements; as a result, the system is forced to evolve towards a desired quantum state. There is a rather large variety of cases depending on how the preparation is achieved; for instance, the measurement does not need to be projective but it can be a generic positive operator-valued measure (POVM) [190]. The protocol can either involve post-selection or be non-selective, as in the rest of this section. Here we briefly discuss the (non-selective) approach developed in Ref. [463] because of its connection with the Lindblad dynamics and its applications to many-body state preparation. In the language we used throughout these lecture notes, the steering protocol designed in [463] can be considered as a (time-periodic) Floquet-Lindblad engineered dynamics.

Following [463], the protocol consists in a sequence of discrete steps in which the detector is initialised and then unitarily coupled to the system for a given time and then decoupled. The detector is then re-initialised. In spirit, this protocol is equivalent to the collisional models of system-bath dynamics (for further discussion of this topic see the review [547]). Indicating the system-detector interaction Hamiltonian with $\hat{H}_I$, and the initial state of the detector by $\hat{\rho}_d$, the state of the system (after tracing over the detector degrees of freedom) evolves as

$$ \hat{\rho}(t + \delta t) = \mathrm{Tr}_d \left[ e^{-i\hat{H}_I t} \hat{\rho}(t) \otimes \hat{\rho}_d e^{i\hat{H}_I t} \right]. \qquad (167) $$

This elementary step plus the re-initialisation reduces, in the continuous time limit, to a Lindblad dynamics (the proof is done with techniques that are reminiscent of the standard methods
that are used to derive the Lindblad master equation, as we do in appendix A).

It is useful to consider first the case of a spin-1/2 system. By choosing a $\hat{H}_I = J\hat{o}_s^+\hat{o}_d^- + \text{H.c.}$ (the indices d, s refer to the system and detector qubits) it is possible to show that the state of the system will converge to an eigenstate along the z-direction in the Bloch sphere.

The same scheme, with a more complex choice of the number of detector qubits and their interaction with the system, can be followed to generate the ground state of the AKLT Hamiltonian. As detailed in Ref. [463], for each pair of neighbouring spins, one has to couple five detector qubits initialised in a fully polarised state. The number of detector qubits is dictated by the complex form of the interaction Hamiltonian that will lead to transition in the different spin sectors. Without entering into details, it is worth recalling that the choice of the interaction Hamiltonian $\hat{H}_I$ should satisfy the following condition

$$ \langle \Psi | \hat{\rho}(t + \delta t) | \Psi \rangle > \langle \Psi | \hat{\rho}(t) | \Psi \rangle, $$

meaning that at each step of the protocol, the state of the system get closer to the target state. This guarantees that the state approaches progressively the target state but does not determine the speed at which this occurs. To this end, active control strategies that choose the appropriate system-detector interaction, based on the measurement record, may speed-up the convergence [548]. We note that related type of protocols have been recently implemented with superconducting circuits with the goal of cooling a many-body system [549].

To conclude, it may be useful to draw some comparison of this steering protocols with the dissipative state preparation discussed in the previous sections. In all the approaches, the system is driven towards the target states by a sequence of measurements, these are either performed by a measurement apparatus (as discussed here) or by the environment (previous sections). While dissipative engineering is entirely based on the identification of the form of the jump operators that have the desired (many-body) state as a target state, in steering protocols the measurement is generically adaptive (in some sense the Kraus operators may depend on time). Furthermore, another difference is that bath engineering is a non-selective protocol (the state is progressively purified), whereas steering may require post-selection. When this last condition applies, steering becomes a special case and a first example of monitored quantum system, whose dynamics will be addressed in more detail in Section 8.

## 5.6 Stability to perturbations

We conclude this section by discussing an important problem associated to the physics of many-body dissipative state engineering, namely that of the stability of the protocol to perturbations. In Ref. [464] the authors discuss the effect of a perturbation on the preparation of a Bose-Einstein condensate discussed in Sec. 5.2, finding a phase diagram that displays a dissipative phase transitions. Dissipative phase transitions are the subject of Sec. 6, so this subsection, that summarises some of the results in Ref. [464], constitutes a bridge toward the material that will be presented there.

We consider the situation in which the purely dissipative dynamics generated by the jump operators in Eq. (148) is complemented by the unitary evolution associated with the standard Bose-Hubbard Hamiltonian for a d-dimensional lattice, first introduced in Eq. (36):

$$ \hat{H} = -t_H \sum_{\langle i,j \rangle} (\hat{a}_i^\dagger \hat{a}_j + \text{H.c.}) + \frac{U}{2} \sum_i \hat{n}_i (\hat{n}_i - 1). \quad (168) $$

In this case, the interacting term proportional to U breaks the conditions discussed above for the appearance of a pure dark state: $|BEC_N\rangle$ is not an eigenstate of $\hat{H}$. In this situation one cannot expect a pure dark state to appear, yet, as shown by the authors of Ref. [464] with a mean-field time-dependent Gutzwiller analysis, for weak enough values of U and for weak
loss rates the system remains in a condensed phase with a macroscopic occupation of the $\mathbf{k}=0$ bosonic mode. When the interaction parameter $U$ is increased, however, the system is driven in a thermal phase that does not display off-diagonal long-range order. The result of this study also provides a qualitative understanding of the stability of the system: although the perturbation breaks the condition that is necessary for having a pure dark state, the mixed state that is obtained displays the same features. When the perturbation exceeds a critical value, the stationary state changes structurally and becomes a trivial thermal phase.

Note however that the thermodynamic limit and the $t \to \infty$ limit do not commute, as we will discuss in more detail in Section 6. Indeed, if the system is discussed at any finite time, the properties of the state $\rho(t)$ evolve smoothly in the phase diagram as a function of the parameters. Yet, if one looks at the properties of the stationary states, they feature sharp discontinuities in the expectation values of the $\mathbf{k}=0$ mode.

Other articles have addressed this stability issue in other models; we mention here Ref. [465], which has addressed the competition between two dissipative mechanisms featuring incompatible many-body dark states. Calling $\eta$ the ratio of the decay rates of the two dissipative mechanisms, it is found that the steady state is pure only in the two limiting cases $\eta = 0$ and $\eta \to \infty$. In between, the stationary state is mixed and thus one cannot speak of dark state. Focusing on dissipative mechanisms that engineer a paramagnetic spin state and ferromagnetic spin states, a mean-field analysis reveals the existence of a nonequilibrium phase diagram composed of two phases separated by a second order phase transition that parallels the properties of the thermal phase diagram.

# 6 Dissipative phase transitions

An essential feature of the Lindbladian many-body dynamics is the competition between coherent evolution, generated by the Hamiltonian, and incoherent processes generated by the dissipation. This competition—as well as competition amongst different dissipative or coherent processes—can result in sharp, non-analytic changes in the properties of the stationary state of the system. This leads to the concept of Dissipative Phase Transitions (DPT), corresponding to non-analytic changes to the steady state as one varies some parameter, $\eta$, through a critical value $\eta_c$. We specifically assume the Lindbladian is a smooth function of $\eta$ (as would occur if $\eta$ controls the strength of some specific term in the Hamiltonian or dissipation). A non-analytic change of the steady state despite smooth evolution of the Lindbladian is the open-quantum-system analogue of the behaviour well known in equilibrium phase transitions. For most dissipative phase transitions one can define an order parameter $O$ in terms of an observable $\hat{O}$ computed on the steady state density matrix, $\hat{\rho}_{ss}$ for given $\eta$, i.e. $O(\eta) = \text{Tr}[\hat{\rho}_{ss}(\eta)\hat{O}]$. A DPT occurs when $(\frac{d}{d\eta})^k O(\eta)$ shows singular behaviour (for some value of $k$) at $\eta \to \eta_c$. Since the observable is independent of the control parameter $\eta$, a DPT must imply a change in the structure of the stationary state.

As with other phase transitions, DPT can be first order (discontinuous, see e.g. [237, 283, 426, 450, 550–554]) or second order (continuous, see e.g. [231, 234, 464, 555, 556]). Examples of these are illustrated in Fig. 25, through bifurcation diagrams describing the order parameter for different stationary states (fixed points) of the system. Figures such as these occur in the theory of dynamical systems [557]; they can also apply to describing the stable, metastable, or unstable points of a many-body system, often found from mean-field or semiclassical analyses. Here we use “unstable” to refer to points where any small fluctuation will grow exponentially. The distinction between “stable” and “metastable” is more subtle, and regards the relative rates of switching between such states—we will discuss this in detail in Sec. 6.1.

The form of curves shown in Fig. 25 are very typical for mean-field theories of dissipative
<figure><img src="image_28.png" /></figure>

Figure 25: Examples of types of discontinuous behaviour, in terms of order parameter $O$ as a function of parameter $\eta$. Thick solid lines indicate stable states, thin solid lines metastable states, while dashed lines indicate unstable fixed points. (See text below for further discussion). Case (a) presents the prototypical first-order transition without spontaneous symmetry break (SSB). Case (b) is the prototypical second-order transition with SSB. However, other scenarios can arise, such a first-order transition with SSB, shown in panel (c), or a second-order transition without SSB, shown in panel (d).

phase transitions, showing multiple solutions at each value of $\eta$. Mean-field theory can distinguish stable and unstable fixed points (i.e. the dashed and solid lines in the figure). However, mean-field approaches do not distinguish between the stable (thick) and metastable (thin) lines, thus predicting bistability when multiple stable and metastable states exist. We will discuss below how the spectral theory of phase transitions refines this understanding.

As with other transitions, DPT can be associated with spontaneous symmetry breaking (SSB) in the thermodynamic limit. In the open quantum system context this means there are symmetries under which the equation of motion is invariant but the steady state is not. We will discuss further below how this relates to the spectral properties of phase transitions, as well as how it emerges as one approaches the thermodynamic limit. The prototypical first-order phase transition, analogous to optical bistability, shown in Fig. 25(a) does not involve SSB, while the prototypical second-order transition, Fig. 25(b) does. However, as illustrated in Fig. 25(c,d), one can also have first-order transitions with SSB and second-order transitions without, for examples see [558]. Moreover, phase boundaries can change between first- and second-order through tricritical points, i.e. Fig. 25(b) can evolve to Fig. 25(c); see Refs. [559–561] for some examples of this.

A feature of DPT that is not generally available [562] for other types of transition is the possibility of cases where there is no stationary state, but instead persistent oscillations. This corresponds to a spontaneous symmetry breaking of continuous time-translation symmetry, as in the case of Boundary Time Crystal phases [563] or dissipative limit cycle phases [284,564–567].

Dissipative phase transitions have been the subject of experimental investigations in various platforms. Notable examples include solid-state quantum optics [568–572], ultracold atoms in optical cavities in the context of Dicke superradiance [59, 110, 573], and Rydberg gases demonstrating an absorbing state transition and self-organised criticality [145]. Transitions to time crystal or limit cycle phases have been observed in experiments with ultracold atoms in cavities [574] or free space [575, 576] as well as nanolasers [577].
In this Section we review DPT from different theoretical perspectives, we provide several key examples and discuss experimental evidence in different platforms. A general classification of DPT has been proposed in terms of the properties of the Lindblad superoperator generating the dynamics, see Refs. [203, 204]. A complementary perspective on DPT can be obtained using Schwinger-Keldysh Field Theory [31, 32] which also allows to assess their relation with field-theoretic classification of other—i.e. classical (thermal) and quantum—phase transitions.

## 6.1 General classification: Spectral theory and Schwinger-Keldysh field theory

The formulation of the Lindblad master equation as an eigenvalue problem for the Lindbladian superoperator is suggestive of similarities with the spectral theory of quantum phase transitions: in equilibrium at zero temperature the existence of a quantum phase transition (QPT) is associated to the closing of an energy gap upon changing a control parameter. If this is the case the ground-state becomes degenerate with the first excited state and the nature of the system can change abruptly, either through a first order or higher order non-analyticity. As we discuss below, the same can be said for DPT in terms of eigenvalues of the Lindbladian; in particular, we will see this relates to the properties of eigenvalues with zero real part. This spectral theory of DPT has been the subject of significant interest recently, see for example Refs. [203, 204] and references therein.

As we will also discuss below, DPT (like other phase transitions) exist only in a thermodynamic limit of large system size. We will discuss below how finite-sized systems still show features of the thermodynamic-limit phase transition, and how the spectral properties depend on system size.

### 6.1.1 Transitions without symmetry breaking

We start by considering the case of a first-order DPT at $\eta = \eta_c$ without symmetry breaking, associated to the existence of a discontinuity in some physical observable, e.g. Fig. 25(a), and discuss the consequences of this on the eigenspectrum of the Lindbladian. We first discuss behaviour in the thermodynamic limit and focus on the ensemble-averaged steady-state density matrix. We then consider effects of finite size and finite duration of experiments. As we discuss further below, there are subtle points about how these two limits—which do not commute—behave.

**Spectral theory in the thermodynamic limit.** It is tempting to try to naively apply the construction that works for a ground-state QPT, where the critical point is associated with the exchange of two eigenstates—i.e. the ground-state and the first excited state swap their role at the critical point—with the closing of the energy gap. Adapting this to the Lindbladian however requires care, since eigenstates of the Lindbladian do not generally correspond to states of the system. Specifically, while $\hat{\rho}_0 = \hat{\rho}_{ss}$ is the physical steady state, all other Lindbladian eigenvectors $\hat{\rho}_{\mu \ge 1}$ are traceless and so, as discussed in Sec. 3.1.2, do not describe physical density matrices. The appropriate picture, as developed in Ref. [203], addresses this issue as we discuss below.

One may start by considering a transition in which there are two different stationary states, one for $\eta < \eta_c$ and one $\eta > \eta_c$. We define these states to be called $\hat{\rho}^<$ and $\hat{\rho}^>$ as follows:

$$ \hat{\rho}^< = \lim_{\eta \to \eta_c^-} \hat{\rho}_{ss}(\eta), \quad \hat{\rho}^> = \lim_{\eta \to \eta_c^+} \hat{\rho}_{ss}(\eta), \qquad (169) $$

and crucially $\hat{\rho}^< \neq \hat{\rho}^>$. This suggests an expression such as

$$ \hat{\rho}_{ss}(\eta) = \theta(\eta_c - \eta)\hat{\rho}^< + \theta(\eta - \eta_c)\hat{\rho}^>. \qquad (170) $$
We discuss the behaviour exactly at $\eta = \eta_c$ further below, as the behaviour at this point will only become clear when we consider the effects of finite size and finite time.

In order that the steady state can change discontinuously at $\eta = \eta_c$ it must be the case that at this point there are two steady states of the Lindbladian, i.e. $\hat{\mathcal{L}}_{\eta=\eta_c}[\hat{\rho}^{<}] = \hat{\mathcal{L}}_{\eta=\eta_c}[\hat{\rho}^{>}] = 0$. We will discuss below how the eigenvalues must behave as one moves away from this point. However, based on the statements already made, we may make statements about the behaviour of the eigenstates of the Lindbladian, $\hat{\rho}_0, \hat{\rho}_1$ in the vicinity of $\eta = \eta_c$.

As we have already seen, the steady state density matrix, $\hat{\rho}_0$ is discontinuous across the transition, jumping between $\hat{\rho}^<$ and $\hat{\rho}^>$ either side of the transition. As we will argue below, exactly at the transition we can expect $\hat{\rho}_0$ is a mixture of these states, such as the equal weight mixture ($\hat{\rho}^< + \hat{\rho}^>)/2$. One may note that the equal-weight mixture is not the only possible solution, and whether this steady state takes exactly this form may vary from problem to problem. In contrast to the discontinuous behaviour of $\rho_0$, the first excited state $\hat{\rho}_1$ is continuous, as $\hat{\rho}_1 = \hat{\rho}^> - \hat{\rho}^<$ is the only traceless form that can be made out of the two density matrices $\hat{\rho}^{<,>}$.

**Behaviour of Lindbladian eigenvalues.** By definition, the steady state $\hat{\rho}_0$ is an eigenvector of the Lindbladian with eigenvalue $\lambda_0 = 0$, while $\lambda_1$ determines the asymptotic decay rate as defined above, $\lambda_{ADR} = -\text{Re}\,\lambda_1$. For simplicity, we will focus here on the case where $\lambda_1$ is real; see Sec. 6.2.3 below for the case where this is not true. Near to the critical point $\eta = \eta_c$ one can in general expect that these two eigenvectors $\hat{\rho}_{0,1}$ fully control the long time dynamics [578]. That is, one can expect that no other eigenvectors have such small eigenvalues. Further away from $\eta_c$, one may have that the eigenvalue $\lambda_{ADR}$ associated with $\hat{\rho}_1$ grows to be comparable with other modes, and the two lowest modes can no longer be considered in isolation.

A key question regarding the behaviour near $\eta = \eta_c$ is whether $\lambda_{ADR}$ vanishes only at a point, or vanishes over an extended region. In Fig. 26 we have suggested the latter case. This would correspond to allowing bistability, where e.g. a system prepared in $\hat{\rho}^<$ can survive even when $\eta > \eta_c$ or vice versa. Bistability is a specific case of metastability, where there exists one stable point and one metastable point. We will discuss below how such metastability can be consistent with a unique steady state density matrix $\hat{\rho}_{ss}$ by considering large but finite system sizes. We further illustrate this with a specific example taken from Ref. [203] shown in Fig. 30. We note that this argument differs from that in Ref. [203], which suggests that the spectral gap vanishes only at the point $\eta = \eta_c$.

Exactly at the critical point, one may note that the arguments in the previous section imply that both $\hat{\rho}_0$ and $\hat{\rho}_1$ have eigenvalue zero. As such, this means one cannot uniquely determine what is the steady state. As we discuss below, this is an issue with considering the Lindbladian in the thermodynamic limit. We thus next turn to consider finite system size and finite time, which help resolve some points raised in the discussion so far.

**Effects of finite system size and finite time.** In the discussion so far, we have focused on the behaviour in the thermodynamic limit, since it is in that limit that phase transitions are strictly defined. We comment here on how this relates to behaviour in finite systems. This discussion will also clarify the significance of long-lived metastable states as introduced above.

In a finite system we may note that in general, even at $\eta = \eta_c$, one can expect a unique steady state to exist. Exceptions to this exist for models with strong symmetries where conserved quantities exist, in this case even finite-sized models have multiple steady states. To understand how this unique state arises from the two states $\hat{\rho}^{<,>}$, and the evolution of the eigenvalues of the Lindbladian in the finite-sized case, one can consider switching rates between the two states. The term “switching rates” imply a considering in terms of a trajectory-based approach (as discussed in Sec. 3.3). As noted there, different approaches to unravelling
<figure><img src="image_29.png" /><figcaption>Figure 26: Dissipative phase transition without symmetry breaking, showing evolution of the Lindbladian gap, $\lambda_{ADR}$ (top) and the steady state density matrix (bottom) in the thermodynamic limit. At the critical point $\eta = \eta_c$ the steady state density matrix becomes a mixture of two states, corresponding to the steady state right before and after the transition. Over a range near this transition, the Liouvillian gap vanishes in the thermodynamic limit (see Fig. 30 for an actual example). This leads to a bistable regime, where a system prepared in the wrong state can persist for a long time (a time diverging as system size approaches the thermodynamic limit). Inspired by Ref. [203].</figcaption></figure>

can lead to different results; here we assume an approach such as unravelling to describe homodyne detection. In such an approach one would see trajectories that stay around one of the two states, $\hat{p}^{<,>}$, but then occasionally switch between these two states—for an example, see Ref. [579]. By assuming this switching is a Poissonian process, one can then extract the switching rates from each state $<,>$ to the alternate state. As we discuss below, in the large system size limit, switching rates become small, but also the ratios of switching rates between the states also become small or large. Metastable states arise when switching from a state has a low rate, so that the state can live for a long time, but when the ratio of switching to and from a state means that the metastable state is disfavoured. In contrast, stable states are those favoured by the switching process.

The trajectory picture is not however required to understand the rates. To make the discussion more concrete, one can consider a classical rate equation the time evolution of probabilities of two states $p_{<,>}$ determined by:

$$ \frac{d}{dt} \begin{pmatrix} p_{>} \\ p_{<} \end{pmatrix} = \begin{pmatrix} -r_1 & r_2 \\ r_1 & -r_2 \end{pmatrix} \begin{pmatrix} p_{>} \\ p_{<} \end{pmatrix}, \qquad (171) $$

where $r_1$ describes the rate of switching from $>$ to $<$, and $r_2$ vice-versa. This has the solution:

$$ \begin{pmatrix} p_{>} \\ p_{<} \end{pmatrix} = \frac{1}{r_1 + r_2} \begin{pmatrix} r_2 \\ r_1 \end{pmatrix} + C \begin{pmatrix} 1 \\ -1 \end{pmatrix} e^{-(r_1+r_2)t}, \qquad (172) $$

where $C$ depends on initial conditions. This can be compared to the ensemble-averaged density matrix, which will decay to the stable state at a rate $\lambda_{ADR}$. This simple model thus has $\lambda_{ADR} = r_1 + r_2$, and a steady state determined by the ratio $r_1/r_2$.

If one considers how these switching rates vary with system size, one expects the switching rates to decay with increasing system size, so that all rates vanish in the thermodynamic limit. Physically this occurs as states become macroscopically distinct in the large system size limit, thus suppressing transition rates between them. To give a specific example of how this might occur, consider the case where $r_{1,2} = \exp(-c_{1,2}N)$ where $N$ denotes system size. In this case, both $r_1, r_2$ (and thus $\lambda_{ADR}$) vanish as $N \to \infty$, while $r_1/r_2 = \exp(-(c_1 - c_2)N)$ will either
vanish or diverge as $N \to \infty$ depending on the sign of $c_1 - c_2$. In this scenario, for $c_1 < c_2$, $r_1/r_2 \to \infty$, so $\hat{\rho}_<$ is the stable state, while $\hat{\rho}_>$ is metastable, while for $c_1 > c_2$ the situation is reversed. This example explains how a vanishing Lindbladian gap $\lambda_{ADR}$ over a finite region—which thus implies metastability—can be compatible with a unique steady state. That is, in the thermodynamic limit, metastability does not imply any continuous transfer of weight between the two states or any bimodality. Indeed, only if $r_{1,2}$ have the same functional dependence on $N$ would one expect $r_1/r_2$ to be finite as $N \to \infty$. The distinct behaviour exactly at $\eta = \eta_c$ is that at this point the switching rates become equivalent, making the steady state a mixture of the two competing states. For the example here, where $c_1 = c_2$ one has $r_1 = r_2$ and thus has exactly the equal weight mixture at this critical point.

Considering this example, we may also note that this shows that the limits of long times and large system size do not commute. We may consider the two cases:

**Large system limit first.** In the large system limit, $N \to \infty$ the switching rates between metastable and stable states vanish. Thus there can be multistability where multiple stable or metastable states exist. That is, for any finite $t$ at $N \to \infty$, the solution of Eq. (172) depends on the initial conditions via the constant $C$

**Long time limit first.** In this case, as long as $\lambda_{ADR} = r_1+r_2$ is non-zero one will reach a unique steady state, independent of the initial conditions. The discussion of the thermodynamic limit given above (based on Ref. [203]) assumes this limit.

There are several consequences of the arguments above that are worth noting:

* The fact that the switching rates vanish in the thermodynamic limit means that multi-stability can occur in experiments, i.e. bistability in this specific example. That is, if one prepares a state $\hat{\rho}^<$ but for parameters where $\eta > \eta_c$ (or vice versa), the state $\hat{\rho}^<$ will persist for a long time, and this timescale will diverge as the system size increases. This means that the bistability predicted by mean-field approaches can have a physical meaning. Metastable states are distinguished from unstable states in that they can appear as long-lived states in a finite-sized system, with a lifetime growing with system size.

* This observation also highlights an important point to note about the density matrix. The density matrix represents the ensemble average over many realisations. This can be related to time averages for ergodic systems. The notion of ergodicity means that the dynamics of a system allow it to explore all possible states. As such, for ergodic systems, long time averages of the dynamics will be equivalent to ensemble averages. However, even when ergodic, the timescale required for the system to explore all states may become large—i.e. the time to switch between distinct solutions may grow (even exponentially) with the size of the system. As such single experiments do not necessarily correspond to the ensemble-averaged density matrix. We discuss other consequences of the distinction between ensemble averages and individual trajectories in Sec. 8.

* Even when considering ensemble-averaged measurements, there are still signatures of such metastability. These arise if one considers the properties of two-time correlations. These allow one to ask the question: "If I measure the system in one state at time $t_0$, how likely am I to measure it in the same state at time $t_1$?". This correlation function can give direct access to the Lindbladian gap $\lambda_{ADR}$.

* For a finite system, the finite switching rates (in both directions) lead to a smooth crossover in the finite-size steady state. The width (i.e. range of parameter $\eta$) over which this crossover occurs will shrink as the system size increases. One can thus extrapolate to find the expected behaviour in the thermodynamic limit from finite-size experiments.
**Extracting $\hat{\rho}^{<,>}$ from finite-size simulations.** While the finite-size system has a unique steady state, it can be useful to have an estimate for the metastable states $\hat{\rho}^{<,>}$. Reference [203] provides a way to extract this from the first non-zero Lindbladian eigenvector, $\hat{\rho}_1$. This is based on writing this eigenvector as follows:

$$ \hat{\rho}_1 = \sum_{\mu} p_{\mu} |\psi_{\mu} \rangle \langle \psi_{\mu}| \approx \hat{\rho}^{<} - \hat{\rho}^{>} \quad (173) $$

The first expression is just the eigendecomposition of any Hermitian operator, while the second is the expected form of $\hat{\rho}_1$ in terms of the metastable states. Note that in writing Eq. (173) we have assumed an overall sign in $\hat{\rho}_1$. Since $\hat{\rho}_1$ is a traceless operator, some of its eigenvalues $p_{\mu}$ must be negative, and some positive. However, since $\hat{\rho}^{<,>}$ should be valid density matrices, they must both have a positive spectrum. This then suggests separating the sum over $\mu$ into those terms for which $p_{\mu}$ is positive (or negative) and identify this with $\hat{\rho}^{<}$ (or $\hat{\rho}^{>}$. That is:

$$ \hat{\rho}^{<} = \sum_{\mu | p_{\mu} > 0} p_{\mu} |\psi_{\mu} \rangle \langle \psi_{\mu}|, \quad \hat{\rho}^{>} = \sum_{\mu | p_{\mu} < 0} (-p_{\mu}) |\psi_{\mu} \rangle \langle \psi_{\mu}| \quad (174) $$

After normalisation, these will allow one to extract two valid density matrices $\hat{\rho}^{<,>}$ from a calculation at finite $N$. As noted above, in this argument we assumed the signs of $\hat{\rho}^{<}$, $\hat{\rho}^{>}$ such as in Eq. 173. If we had chosen the opposite signs the role of the two density matrices would be swapped, but the overall conclusions would be unchanged. One can determine which of $\hat{\rho}^{<}$, $\hat{\rho}^{>}$ is which by comparison with the behaviour away from the critical point.

A concrete example of finite-size scaling of the Lindbladian gap, the evolution from a crossover to a sharp transition, and the extraction of $\hat{\rho}^{<,>}$ is given in Sec. 6.2.2, in particular Fig. 30.

### 6.1.2 Symmetry-breaking phase transitions

We now turn to phase transitions that show spontaneous symmetry breaking, such as the examples shown in Fig. 25(b,c). When discussing symmetries in open quantum systems, one must distinguish between strong and weak symmetries, as discussed in Sec. 3.1.1. In this section we focus on cases with weak symmetries. Once again we first discuss the behaviour in the thermodynamic limit and then discuss how this emerges from the behaviour at large but finite system sizes.

**Spectral theory in the thermodynamic limit.** We first discuss discrete symmetries, such as those shown in Fig. 25(b,c) which corresponds to breaking a $\mathbb{Z}_2$ symmetry. Here, for $\eta < \eta_c$ there is a unique steady state, but for $\eta > \eta_c$ the symmetry is spontaneously broken, leading to multiple steady states. While Fig. 25(b,c) suggests two steady states, which we denote $\hat{\rho}^{\pm}$, the full density matrix solution allows any linear combination of these states. Hence, a degenerate family of steady states is possible. This in turn implies that for $\eta > \eta_c$ the Lindbladian must have at least one additional zero eigenvalue (additional to the zero eigenvalue that always exists associated with there being a steady state). This extra zero eigenvalue corresponds to an eigenvector that describes moving within the space of possible steady states.

This argument is similar to that predicting the existence of a Goldstone mode when a continuous symmetry is broken [254,448], however in the open quantum system a zero eigenvalue exists even if a discrete symmetry is broken. When considering systems with continuous symmetries, a similar situation arises, although now one must consider the possibility of infinitesimal changes to the order parameter symmetry, rather than just switching between discrete states. This leads to a larger degenerate subspace, and thus the possibility of multiple zero eigenvalues.
<figure><img src="image_30.png" /><figcaption>Figure 27: Dissipative Phase Transition with symmetry breaking. In this case, multiple steady states (corresponding to different ways to break the symmetry) exist beyond the critical point. This means the Lindbladian gap vanishes throughout the ordered phase [Inspired by Ref. [203]].</figcaption></figure>

To summarise, whenever symmetry breaking occurs, the Lindbladian gap must vanish throughout this phase. This is shown in Fig. 27. The vanishing gap implies that all symmetry broken states are equally likely; this contrasts with the bistable system where (apart from exactly at $\eta = \eta_c$) one or other solution is favoured. One may also note that this is what distinguishes the similar structures in Fig. 25(b,d); in panel (b) the two solutions are equivalent so must have equal tunnelling rates, while in panel (d) they are not equivalent.<fnref n="4" /> Figure 25(c) shows both symmetry breaking and a possible metastable regime. Here zero eigenvalues can be expected to arise from both these features.

**Effects of finite system size.** Similarly to transitions without symmetry breaking discussed above, we may note that in a finite size system, for weak symmetries, there is a unique steady state solution. Moreover, this unique state always respects the symmetries of the equation of motion. This means that in a finite system, as viewed from the properties of the steady state density matrix, there is no spontaneous symmetry breaking. We can however again gain further insight into how the thermodynamic limit emerges by considering the behaviour of the Lindbladian spectrum.

For the case of discrete symmetry breaking we may consider the behaviour of the two lowest eigenvectors of the Lindbladian. As noted above, the steady state will respect the symmetry, and so must take the form $\hat{\rho}_0 = (\hat{\rho}^+ + \hat{\rho}^-)/2$ in terms of the two symmetry-broken solutions $\hat{\rho}^\pm$ introduced above. The other mode (that would have zero eigenvalue in the thermodynamic limit) is the traceless anti-symmetric combination $\hat{\rho}_1 = (\hat{\rho}^+ - \hat{\rho}^-)/2$. This mode corresponds to an eigenvalue that vanishes as system size diverges. The consequence of this mode is that if one prepares one symmetry-broken state at a point where symmetry breaking would occur in the thermodynamic limit, then the finite-size system will only recover the symmetric state very slowly. One can think of this in the same way as for the first-order transition, with switching rates between the two symmetry broken states. A key difference in the symmetry breaking case is that due to symmetry, the switching rates between different metastable states $\hat{\rho}^{<,>}$ are exactly equal. This means that $\hat{\rho}_0$ is always the symmetric sum of symmetry-broken states for all $\eta > \eta_c$, not only at the critical point.

As with the case without symmetry breaking, multi-time correlations can again be used to directly probe the Lindbladian gap. Here, in the ordered state, one can consider preparing one symmetry broken state at $t_0$ and asking whether one finds the same state at a later time $t_1$. A paradigmatic example of this (for continuous symmetry breaking) is the Schawlow-Townes linewidth of a laser [184, 580] (See also Sec. 2.6). Here, the thermodynamic limit of a laser corresponds to breaking a continuous symmetry corresponding to the phase of the emitted

<footnote n="4">⁴We are not aware of particular examples that show the behaviour of panel (d), however there is no obvious reason they should not exist.</footnote>
light; the linewidth depends on how phase correlations decay with time, thus are a signature of the Lindbladian gap. The fact the linewidth is non-zero indicates a finite correlation time, as expected away from the thermodynamic limit.

### 6.1.3 Critical behaviour and relation to Schwinger-Keldysh field theory

Schwinger-Keldysh field theory – as introduced in Sec. 3.2 and Sec. 4.7 – provides a complementary approach to understanding dissipative phase transitions, in particular continuous (i.e. second-order) transitions, and a powerful formalism to study the resulting critical behaviour. Such an approach is reviewed extensively elsewhere [31, 32]. Here we briefly mention some key points, as well as discussing how such an approach links to the spectral theory discussed above.

The key approach to describing second-order phase transitions is the renormalisation group (RG) as introduced in Sec. 4.7.3. The idea underlying this is that at a critical point (i.e. at a second-order phase transition), correlation lengths and timescales diverge [254, 448]. Correlation and response functions (as discussed below) typically depend on a dimensionless length or time, set by the ratio of separation to a typical length and timescale. As such, when these scales diverge, correlation functions take universal forms. The RG approach describes these universal forms on long length and timescales, by integrating out the effects of short wavelength modes. Furthermore, one finds that under the RG, the asymptotic long wavelength form the Keldysh action is defined by the symmetry of the system and the dimensionality. This leads to a notion of Universality: the long wavelength effective actions can be divided into “universality classes” (determined by symmetry and dimensionality), and these in turn determine critical behaviour such as critical exponents [32].

In discussing critical exponents and correlation lengths, these refer specifically to properties of correlation and response functions. That is, the field theory is considered as a route to calculate the (long wavelength) form of a set of response functions. Because of universality, there are only a finite set of response functions that need to be found to fully determine the behaviour. One such correlation function is susceptibility $\chi$, which determines how the order parameter $O$ behaves as a function of a symmetry breaking term in the Hamiltonian, $\hat{H} \rightarrow \hat{H} + f\hat{O}$, i.e. $\chi = (\partial O / \partial f)|_{f \rightarrow 0}$. For a second-order transition, this susceptibility diverges algebraically as $\eta \rightarrow \eta_c$, allowing one to define a critical exponent $\gamma$ via $\chi \propto (\eta_c - \eta)^{-\gamma}$. Further exponents can be defined by considering time-dependent fields and the frequency dependence, or finite-size systems, and the dependence on lengthscales. Other susceptibilities may also be important.

To connect such critical behaviour to spectral theory, we must consider how correlation functions or response functions such as susceptibility are related to the Lindbladian spectrum. For correlation functions, we may note that (as discussed above) multi-time correlation functions can in principle reveal information about the Lindbladian gap. Using the quantum regression theorem [184, 191, 216], a correlation function $\langle \hat{X}(t) \hat{X}(0) \rangle$ in steady state can be found by taking an initial state $\hat{X} \hat{\rho}_{ss}$, time evolving this state to time $t$ and then evaluating the expectation value. If one has a full eigensystem of the Lindbladian, i.e. eigenvalues and left- and right-eigenvectors $\lambda_\mu$, $\langle\langle l_\mu | , | r_\mu \rangle\rangle$, as described in Sec. 3.1.2, the time evolution will take the form:

$$ \langle \hat{X}(t) \hat{X}(0) \rangle = \sum_{\mu} e^{\lambda_{\mu} t} \langle \langle \mathbb{1} | \hat{X} r_{\mu} \rangle \rangle \langle \langle l_{\mu} | \hat{X} \rho_{ss} \rangle \rangle. \quad (175) $$

One thus sees that only those eigenvalues $\lambda_\mu$ for which $\langle\langle l_\mu | \hat{X} \rho_{ss} \rangle\rangle \neq 0$ will contribute. For response functions, such as susceptibility, these can be calculated using (non-Hermitian) perturbation theory, again making use of the eigensystem of the Lindbladian. If one considers the vectorised Lindbladian in the form $\hat{L} \rightarrow \hat{L} + \hat{L}_{pert}$, and calculates the expectation of an
observable $\hat{O}$, this gives a steady state expectation:

$$ \langle \hat{o} \rangle = \langle \langle 1 | \hat{o} \rho_{ss} \rangle \rangle - \sum_{\mu \ge 1} \frac{1}{\lambda_{\mu}} \langle \langle 1 | \hat{o} r_{\mu} \rangle \rangle \langle \langle l_{\mu} | \hat{\mathcal{L}}_{pert} \rho_{ss} \rangle \rangle . \quad (176) $$

This shows directly how diverging susceptibility is related to a vanishing Lindbladian gap. Naively this might suggest the critical exponent is determined by how $\lambda_{ADR}$ depends on $\eta_c - \eta$. However, in general this expression can contain a contribution from many Lindbladian eigenvalues that vanish (or become small) at the critical point. Summing over their contributions can lead to non-trivial (non-integer) critical exponents.

It is worth noting that the spectral theory of transitions discussed in previous sections focuses on the behaviour of a few eigenvectors of the Lindbladian; in particular the properties of the eigenvector corresponding to the Lindbladian gap. This single eigenvalue will determine the critical properties of transitions in zero-dimensional systems, or those where mean-field theory applies. In general, the critical properties will depend on the spectrum of eigenvalues corresponding to low energy effective modes—i.e. on a set of many eigenvalues, with the number diverging with system size.

One well studied example of this occurs when considering the long-wavelength modes of the driven-dissipative weakly-interacting dilute Bose gas (WIDBG), Eq. (38). These can be found in various ways, including by considering fluctuations around the steady state of the complex Gross-Pitaevskii equation, Eq. (109), which—as discussed earlier—can also be seen as the saddle point of the Schwinger–Keldysh path integral. To see this, one starts by writing the Bogoliubov–de Gennes parameterisation

$$ \psi(r, t) = \left( \psi_0 + \sum_k u_k e^{-i\nu_k t + i\mathbf{k}\cdot\mathbf{r}} + \nu_k^* e^{i\nu_k^* t - i\mathbf{k}\cdot\mathbf{r}} \right) e^{-i\mu t}, $$

where $\psi_0 e^{-i\mu t}$ is a stationary solution, so satisfies $\mu = N U |\psi_0|^2$ and $\kappa_+ - \kappa_- - \kappa_-^{(2)} N |\psi_0|^2 = 0$. One then substitutes this into Eq. (109) and expands up to linear order in the fluctuations $u_k, \nu_k$ giving coupled equations for the fluctuations:

$$ \nu_k u_k = \left( \frac{k^2}{2m} + \mu - i\kappa_{\text{eff}} \right) u_k + (\mu - i\kappa_{\text{eff}}) \nu_k, \quad (177) $$

$$ -\nu_k \nu_k = \left( \frac{k^2}{2m} + \mu + i\kappa_{\text{eff}} \right) \nu_k + (\mu + i\kappa_{\text{eff}}) u_k, \quad (178) $$

where $\kappa_{\text{eff}} = \frac{1}{2} \kappa_-^{(2)} N |\psi_0|^2 = \frac{1}{2}(\kappa_+ - \kappa_-)$. This then leads to the complex dispersion:

$$ \nu_k = -i\kappa_{\text{eff}} \pm \sqrt{\frac{k^2}{2m} \left( \frac{k^2}{2m} + 2\mu \right) - \kappa_{\text{eff}}^2}. \quad (179) $$

In the limit $\kappa_{\text{eff}} = 0$ this recovers the standard Bogoliubov dispersion of the equilibrium WIDBG. For the dissipative case, one finds that the dispersion at small $k$ is purely imaginary and diffusive, i.e. $\nu_k \approx -iD_{\text{eff}}k^2 + \mathcal{O}(k^4)$ for the + root, and $\nu_k \approx -2i\kappa_{\text{eff}}$ for the – root. As $k$ increases, there is an exceptional point, beyond which the imaginary parts of the two roots coalesce, and a real part of the dispersion arises [158, 159]. Such behaviour is characteristic of the Goldstone mode in many driven-dissipative systems.

Recent experiments on exciton polaritons have observed such a diffusive Goldstone mode [581]. Similar results also apply for cold atoms in multimode optical cavities, where the Goldstone mode associated with crystallisation—i.e. a phonon mode—has been observed [115]. The distinct properties of the diffusive Goldstone mode can change the critical behaviour associated with symmetry breaking transitions in such driven dissipative systems.

Renormalisation group approaches are a key tool to account for the effects of this continuum of low energy modes. This is discussed further in the reviews [31, 32, 163].
<figure><img src="image_31.png" /><figcaption>Figure 28: Nonequilibrium mean-field phase diagram of the XYZ model as a function of the anisotropy $J_x \neq J_y$, obtained respectively with $J_z = \gamma$ (left) and $J_z = 0$ (right). We see that along the line $J_x = J_y$ the model displays a regular paramagnetic phase which gives rise to a rich pattern of broken symmetry phases in presence of anisotropy. From Ref. [234] [Copyright (2013) by the American Physical Society].</figcaption></figure>

## 6.2 Examples

Having discussed the general properties of dissipative phase transitions, we here give a few illustrative examples of models in which dissipative phase transitions have been studied.

### 6.2.1 Dissipative phase transitions in spin chains

The first example we consider is spin-chain models. A notable feature these models illustrate is that DPT can display patterns of symmetry breaking which are completely different from the equivalent Hamiltonian in thermal equilibrium. A particularly clear application of this is presented in Ref. [234] for a dissipative XYZ model. This model (as already introduced in Sec. 3.1.5) is given by the Hamiltonian

$$ \hat{H} = \sum_{\langle i,j \rangle} \sum_{\alpha} J_{\alpha} \hat{\sigma}_{i}^{\alpha} \hat{\sigma}_{j}^{\alpha}, \qquad (28) $$

along with jump operators

$$ \hat{L}_{i-} = \sqrt{\gamma} \hat{\sigma}_{i-}. \qquad (180) $$

The ground state of this Hamiltonian is a standard example of quantum magnetism in condensed matter physics. Except in one dimension, the ground state of this model shows long-range magnetic order. The nature of that order—ferromagnetic vs antiferromagnetic—is determined by the sign of the dominant coupling $J_x, J_y, J_z$. While some patterns of coupling can lead to more complex phenomena such as ferrimagnetic order, small changes to ratio of $J_x/J_y$ do not generally change the nature of the order unless the signs change. As we discuss below, and as illustrated by the phase diagram in Fig. 28, this is quite different for the dissipative model.

In Ref. [234], this model was studied using a real-space mean-field factorisation. This leads
to the following equations of motion for the local magnetisation $\langle \hat{\sigma}_i^\alpha \rangle$:

$$ \frac{d\langle\hat{\sigma}_i^x\rangle}{dt} = -\frac{\gamma}{2}\langle\hat{\sigma}_i^x\rangle + \sum_j (J_y\langle\hat{\sigma}_i^z\rangle\langle\hat{\sigma}_j^y\rangle - J_z\langle\hat{\sigma}_i^y\rangle\langle\hat{\sigma}_j^z\rangle), \quad (181a) $$

$$ \frac{d\langle\hat{\sigma}_i^y\rangle}{dt} = -\frac{\gamma}{2}\langle\hat{\sigma}_i^y\rangle + \sum_j (J_z\langle\hat{\sigma}_i^x\rangle\langle\hat{\sigma}_j^z\rangle - J_x\langle\hat{\sigma}_i^z\rangle\langle\hat{\sigma}_j^x\rangle), \quad (181b) $$

$$ \frac{d\langle\hat{\sigma}_i^z\rangle}{dt} = -\gamma \sum_i (1 + \langle\hat{\sigma}_i^z\rangle) + \sum_j (J_x\langle\hat{\sigma}_i^y\rangle\langle\hat{\sigma}_j^x\rangle - J_y\langle\hat{\sigma}_i^x\rangle\langle\hat{\sigma}_j^y\rangle). \quad (181c) $$

This dynamics always admits a steady state, corresponding to $\langle\hat{\sigma}_i^x\rangle = \langle\hat{\sigma}_i^y\rangle = 0$ and $\langle\hat{\sigma}_i^z\rangle = -1$. That is, a state with all the spins down along the z axis, with this axis chosen by the dissipative processes. However linear stability analysis around this steady state reveals a rich pattern of symmetry broken phases, distinct from those in the ground state of the Hamiltonian.

As discussed in Sec. 3.1.5, to reach a non-trivial state requires we consider $J_x \neq J_y$. This is because for $J_x = J_y = J$ the Hamiltonian can be rewritten as

$$ \hat{H} = \sum_{ij} J (\hat{\sigma}_i^+ \hat{\sigma}_j^- + \text{H.c.}) + J_z \hat{\sigma}_i^z \hat{\sigma}_j^z, \quad (182) $$

which conserves the total number of excitations $\sum_i \hat{\sigma}_i^z$. Since the dissipation leads to loss of excitations, this means that for $J_x = J_y$ the only steady state is the dark state $\hat{\rho} = |\downarrow \cdots \downarrow\rangle\langle\downarrow \cdots \downarrow|$.

For $J_x \neq J_y$ excitation number is not conserved because of the reduced symmetry. This allows coherent excitation creation to balance the loss. This process can be seen from the mean-field dynamics which can be described as a competition between dissipative processes, taking the local spin toward pointing down, and precession around a time-dependent field. The effective field seen by each site comes from the neighbouring sites and takes the form

$$ h_i^{\text{eff}} = \sum_{j \in \partial_i} (J_x \langle \hat{\sigma}_j^x \rangle, J_y \langle \hat{\sigma}_j^y \rangle, J_z \langle \hat{\sigma}_j^z \rangle), $$

using the notation of Sec. 4.4.1. Now, the key point is that for isotropic exchange and homogeneous order this field is always parallel to the spin direction itself, so there is no precession and dissipation always wins. However for $J_x \neq J_y$ the spin can have an arbitrary angle with respect to the field and therefore counterbalance the losses, so that new stationary states can emerge.

The discussion so far has been based on the mean-field theory, which, as discussed in Sec. 4.5 is expected to be valid in sufficiently high dimensions. Following Ref. [234] there have been many attempts to explore the phase diagram of this model beyond mean-field theory. One such example, Ref. [283], used the Cluster Mean-Field techniques described in Sec. 4.6.1. This revealed a striking effect arising from the nonequilibrium nature of the problem: short-ranged correlations, disregarded in the single-site mean-field approximation leading to Eq. (181a)–(181c), induce crucial changes to the topology of the phase diagram. In particular, the ferromagnetic phase that previously extends to arbitrarily large $J_y, J_x$ is instead found only for a finite range of couplings, and instead a re-entrant paramagnetic phase is seen at larger couplings. This is shown in Fig. 29.

Such significant changes due to short-range interactions would be highly unusual in ground-state or thermal-equilibrium phase diagrams. The origin of this effect can be seen to be rooted in the driven dissipative nature of the problem. In fact, in the limit of $J_y \to \infty$, the steady-state of the system becomes fully mixed and the purity drops to zero. This in turn leads to suppression of the magnetisation.