# Many-body open quantum systems

© Rosario Fazio¹,², © Jonathan Keeling³, © Leonardo Mazza⁴ and © Marco Schirò⁵

1 The Abdus Salam International Center for Theoretical Physics (ICTP),
Str Costiera 11, I-34151 Trieste, Italy

2 Dipartimento di Fisica “E. Pancini”, Università di Napoli Federico II,
Monte S. Angelo, I-80126 Napoli, Italy

3 SUPA, School of Physics and Astronomy, University of St Andrews,
KY16 9SS St Andrews, United Kingdom

4 Université Paris-Saclay, CNRS, LPTMS, 91405, Orsay, France

5 JEIP, UAR 3573 CNRS, Collège de France, PSL Research University,
11 Place Marcelin Berthelot, 75321 Paris Cedex 05, France

<figure><img src="image_1.png" /></figure>

Part of the ICTP Lecture Notes Series

## Abstract

These lecture notes discuss the recent theoretical advances in the understanding of open quantum many-body physics in platforms where both dissipative and coherent processes can be tuned and controlled to a high degree. We start by reviewing the theoretical frameworks and methods used to describe and tackle open quantum many-body systems. We then discuss the use of dissipative processes to engineer many-body stationary states with desired properties and the emergence of dissipative phase transitions arising out of the competition between coherent evolution and dissipation. We review the dynamics of open quantum many body systems in the presence of correlated many-body dissipative processes, such as heating and many-body losses. Finally we provide a different perspective on open quantum many-body systems by looking at stochastic quantum trajectories, relevant for the case in which the environment represents a monitoring device, and the associated measurement-induced phase transitions.

<figure><img src="image_2.png" /></figure>

Copyright R. Fazio et al.

This work is licensed under the Creative Commons Attribution 4.0 International License.

Published by the SciPost Foundation.

Received 2024-10-15

Accepted 2025-05-12

Published 2025-08-14

doi:10.21468/SciPostPhysLectNotes.99

<figure><img src="image_3.png" /></figure>

Check for
updates

## Contents

1 **Introduction** 4

1.1 Scope and structure of the lecture notes 4

1.2 Historical perspective 6

1.3 Relation with thermal equilibrium 8

1.4 Notation 12

2 **Experimental platforms** 13

2.1 Trapped ions 13

2.2 Ultracold atoms 15
<table><tr><td>2.3</td><td>Cold atoms in optical cavities</td><td>17</td></tr><tr><td>2.4</td><td>Rydberg gases</td><td>20</td></tr><tr><td>2.5</td><td>Superconducting circuits</td><td>21</td></tr><tr><td>2.6</td><td>Polaritons and photons</td><td>23</td></tr></table>

## 3 Theoretical frameworks 27

<table><tr><td>3.1</td><td>Lindblad master equation</td><td>27</td></tr><tr><td>3.1.1</td><td>Symmetries</td><td>28</td></tr><tr><td>3.1.2</td><td>Vectorisation, eigenstates and steady states</td><td>29</td></tr><tr><td>3.1.3</td><td>Multi-time correlations</td><td>31</td></tr><tr><td>3.1.4</td><td>Thermalisation</td><td>32</td></tr><tr><td>3.1.5</td><td>Examples of models</td><td>34</td></tr><tr><td>3.2</td><td>Schwinger-Keldysh path-integral formalism</td><td>37</td></tr><tr><td>3.3</td><td>Quantum trajectories</td><td>40</td></tr><tr><td>3.4</td><td>Additional theoretical frameworks</td><td>43</td></tr><tr><td>3.4.1</td><td>Phase space approaches</td><td>43</td></tr><tr><td>3.4.2</td><td>Heisenberg-Langevin and adjoint equations</td><td>44</td></tr><tr><td>3.4.3</td><td>Third quantisation</td><td>45</td></tr><tr><td>3.4.4</td><td>Effective non-Hermitian Hamiltonian</td><td>46</td></tr></table>

## 4 Overview of theoretical methods 47

<table><tr><td>4.1</td><td>Exact solutions</td><td>47</td></tr><tr><td>4.2</td><td>Cluster expansions</td><td>51</td></tr><tr><td>4.3</td><td>Perturbative expansion</td><td>54</td></tr><tr><td>4.4</td><td>Mean-field approximations</td><td>55</td></tr><tr><td>4.4.1</td><td>Real space factorisation and Gutzwiller mean-field theories</td><td>56</td></tr><tr><td>4.4.2</td><td>Momentum space factorisation</td><td>58</td></tr><tr><td>4.4.3</td><td>Particle-based factorisation: Hartree-Fock theory and Gross-Pitaevskii equation</td><td>60</td></tr><tr><td>4.5</td><td>Large connectivity limit and dynamical mean field theory</td><td>60</td></tr><tr><td>4.5.1</td><td>Large connectivity and validity of mean-field theories</td><td>61</td></tr><tr><td>4.5.2</td><td>DMFT for driven-dissipative bosonic lattices</td><td>62</td></tr><tr><td>4.5.3</td><td>DMFT for driven-dissipative fermions</td><td>65</td></tr><tr><td>4.6</td><td>Systematic expansions beyond mean-field theory.</td><td>65</td></tr><tr><td>4.6.1</td><td>Cluster mean field theory.</td><td>66</td></tr><tr><td>4.6.2</td><td>Tensor-network methods.</td><td>66</td></tr><tr><td>4.6.3</td><td>Neural network Ansatz</td><td>71</td></tr><tr><td>4.6.4</td><td>Corner-space renormalisation</td><td>73</td></tr><tr><td>4.6.5</td><td>Cumulant expansion.</td><td>73</td></tr><tr><td>4.6.6</td><td>Multi-configuration time-dependent Hartree (MCTDH).</td><td>74</td></tr><tr><td>4.7</td><td>Schwinger-Keldysh path integral methods.</td><td>75</td></tr><tr><td>4.7.1</td><td>Saddle point and mean field theory.</td><td>76</td></tr><tr><td>4.7.2</td><td>Quadratic fluctuations and ladder series expansions.</td><td>77</td></tr><tr><td>4.7.3</td><td>Renormalisation group.</td><td>77</td></tr><tr><td>4.7.4</td><td>Stochastic interpretations.</td><td>78</td></tr></table>

## 5 Dissipative state engineering 79

<table><tr><td>5.1</td><td>From optical pumping to dark states</td><td>80</td></tr><tr><td>5.2</td><td>Many-body dark states</td><td>81</td></tr><tr><td>5.2.1</td><td>Nature of the steady states</td><td>82</td></tr><tr><td>5.2.2</td><td>Asymptotic decay toward the dark states</td><td>83</td></tr></table>
<table><tbody><tr><td>5.3 Examples</td><td>83</td></tr><tr><td>5.3.1 BEC condensate</td><td>83</td></tr><tr><td>5.3.2 Topological states and Majorana dark states</td><td>84</td></tr><tr><td>5.3.3 Spin-Dicke states of fermionic gases</td><td>87</td></tr><tr><td>5.3.4 Dissipative preparation of Mott insulators of photons</td><td>90</td></tr><tr><td>5.4 Quantum computation using dissipation</td><td>92</td></tr><tr><td>5.5 Quantum many-body steering</td><td>93</td></tr><tr><td>5.6 Stability to perturbations</td><td>94</td></tr><tr><td><b>6 Dissipative phase transitions</b></td><td><b>95</b></td></tr><tr><td>6.1 General classification: Spectral theory and Schwinger-Keldysh field theory</td><td>97</td></tr><tr><td>6.1.1 Transitions without symmetry breaking</td><td>97</td></tr><tr><td>6.1.2 Symmetry-breaking phase transitions</td><td>101</td></tr><tr><td>6.1.3 Critical behaviour and relation to Schwinger-Keldysh field theory</td><td>103</td></tr><tr><td>6.2 Examples</td><td>105</td></tr><tr><td>6.2.1 Dissipative phase transitions in spin chains</td><td>105</td></tr><tr><td>6.2.2 Coherent drive and bistability</td><td>107</td></tr><tr><td>6.2.3 Boundary time-crystal and dissipative limit-cycles</td><td>110</td></tr><tr><td>6.2.4 Absorbing state phase transitions</td><td>113</td></tr><tr><td><b>7 Dynamics of approach to the steady state</b></td><td><b>115</b></td></tr><tr><td>7.1 Heating dynamics under dephasing</td><td>115</td></tr><tr><td>7.1.1 Quantum spin chains under dephasing</td><td>116</td></tr><tr><td>7.1.2 Bose-Hubbard model under dephasing</td><td>118</td></tr><tr><td>7.2 Depletion dynamics of bosonic gases under many-body losses</td><td>121</td></tr><tr><td>7.2.1 Ideal Bose gas and weak two-body losses</td><td>122</td></tr><tr><td>7.2.2 Fermionised regime and quantum-Zeno effect</td><td>123</td></tr><tr><td>7.2.3 Related results</td><td>127</td></tr><tr><td><b>8 Monitored quantum systems</b></td><td><b>127</b></td></tr><tr><td>8.1 Trajectory-resolved averaging</td><td>127</td></tr><tr><td>8.2 Measurement-induced phase transitions</td><td>129</td></tr><tr><td>8.3 Post-selection problem</td><td>130</td></tr><tr><td>8.4 Experimental evidence for MIPT</td><td>131</td></tr><tr><td><b>9 Summary and outlook</b></td><td><b>131</b></td></tr><tr><td><b>A Derivation of Lindblad master equation</b></td><td><b>135</b></td></tr><tr><td>A.1 From quantum channels to the Lindblad equation</td><td>135</td></tr><tr><td>A.2 Born-Markov master equation</td><td>136</td></tr><tr><td><b>References</b></td><td><b>138</b></td></tr></tbody></table>
# 1 Introduction

## 1.1 Scope and structure of the lecture notes

Understanding the effect of an external environment on the dynamics of quantum system has been a problem of paramount importance since the early days of quantum mechanics [1] because it touches questions that are at the heart of its fundamental principles [2,3]. Much more recently, with the birth of quantum information processing, the properties of open quantum systems have become very relevant, as decoherence and damping are the ultimate limiting factors for the development of quantum technologies [4,5]. Quantum systems in contact with an external environment are referred to as *open quantum systems* and are described by a quantum state that is, generically, mixed. For this purpose, several techniques, both in the design and in the operating mode of qubits and gates, have been implemented with the scope of reducing unavoidable environmental effects [6,7].

If on one side an external environment typically tends to hinder quantum coherence effects, on the other side several new phenomena emerge from the competition between unitary evolution and dissipative effects. This phenomenology gets considerably richer in the case of many-body systems. Exotic phases of matter that are not allowed at equilibrium, novel dynamical regimes in the evolution of complex quantum systems or critical phenomena controlled by dissipation are only few of the most exciting examples. Furthermore, and contrary to what is generally believed, dissipation is not always detrimental for manipulating quantum systems: in the same way as optical pumping is a powerful tool for tailoring the properties of quantum systems at the single-atom level, it is also possible to control the dynamics of many-body quantum systems by suitably engineered dissipation. In fact, with the help of only dissipative dynamics, one can perform a universal quantum computation in a way that is completely equivalent to the standard formulation in terms of unitary quantum circuits [8,9].

This set of lecture notes aims to give a pedagogical introduction to the physics of *many-body open quantum systems*. Obviously this is a broad research field, that essentially includes every aspect of the equilibrium and nonequilibrium statistical mechanics of quantum systems coupled to baths. One could even go further by noting that, for ergodic systems, the system itself can play the role of the bath for one of its subparts and conclude that the recent studies on the unitary dynamics of closed quantum systems and the problem of thermalisation belong to this subject. Given the potential breadth of such a field, these lecture notes will necessarily have a more narrow focus, as we outline next. First of all, our interest will be on *Synthetic Quantum Matter*—artificial quantum complex systems which can be used to realise tunable and controllable phases of matter. In such platforms, one naturally faces the presence of a bath, but one could also engineer the bath itself and its coupling with the setup, and mix its action with an external drive, i.e. *driven-dissipative systems*. Monitoring the system by performing measurements is also a way to perturb the unitary dynamics that we will consider and that is gaining increasing experimental interest. Secondly, in these lecture notes, we will be studying those situations in which the external bath *modifies* the collective behaviour of synthetic quantum systems, for instance by inducing new phases, provoking novel kinds of critical properties, or transforming the conditions under which various states of matter can be found. Thirdly, we will consider specifically Markovian systems, whose dynamics could be described by a time-local Lindblad master equation, with a term that accounts for the unitary evolution due to the Hamiltonian and a contribution that describes the effect of the external bath, responsible for the relaxation to a steady state and for the loss of quantum coherence.

Synthetic open quantum systems also allow one to explore phenomena that are otherwise rare (or non-existent) in nature, even though they are possible according to the laws of quantum physics. They can also be seen, following Feynman's vision [10], as open-system quantum
simulators. In the last decades there has been a considerable experimental progress in realis-
ing platforms where many-body effects in the presence of dissipation can be explored, ranging
from optical lattices and trapped ions to Bose-Einstein condensates in cavities, from cavity
arrays to quantum circuits subject to measurements. In essentially all the cases mentioned
above, the dynamics of the density matrix of the system can be accurately described by a Lind-
blad equation. The unitary part typically represents a many-body Hamiltonian (for example
the exchange Hamiltonian for interacting spin-systems, or the Bose-Hubbard Hamiltonian for
atoms in optical lattices) or a sequence of quantum gates; on its own, it will induce correla-
tions and entanglement between local degrees of freedom. The presence of the environment
or of measurements will result in a local incoherent dynamics that may compete against such
ordering: this interplay will strongly affect both the steady state and the transient dynamics.
Our goal is to give a broad and pedagogical introduction to these phenomena.

Understanding the influence of decoherence and dissipative effects in many-body quantum
systems touches many questions that are beyond the scope of these lecture notes. Here we
note some of these topics, and provide references to reviews that discuss these points. We will
concentrate on systems whose dynamics is time-local (Markovian) and specifically of Lindblad
form, thus we will not consider non-Markovian open systems [11]. We will further focus
our attention on those phenomena where dissipation is an essential element, rather than a
problem to be avoided. This means that we will not cover the problems related to the reliability
of quantum simulators [12, 13] in the presence of imperfections and noise [5-7]. Focusing
on many-body systems we will also put to one side the rich topic of spin-boson dynamics
and associated models [14, 15] (other than brief remarks in the next section), and similarly
we will ignore questions about the role of dissipation and decoherence in the emergence of
classical behaviour [2, 3]. While we will discuss some aspects of exciton-polariton physics—in
particular, methods first developed in the context of modelling polariton condensates—we will
not aim to provide a comprehensive overview of this large field. Reviews of this field can be
found in a number of places, e.g. Refs. [16-19].

Likewise, we will not discuss those situations where few qubits are coupled to two or more
reservoirs (as in quantum thermal machines [20,21] or transport setups [22,23]), nor will we
explore the whole field of quantum thermodynamics [24,25]. We also refer to more specialised
reviews for the exciting problem of developing efficient numerical tools for the simulation of
open quantum systems [26]. Open quantum systems are often described in terms of effective
non-Hermitian models for pure states: we will not enter into this subject and stick to the
formalism based on density matrices and Lindblad master equations, referring the interested
reader to existing reviews [27]. Although very interesting and very important for quantum
technological applications, these topics will not be considered further below, except when of
relevance to the many-body realm.

These lecture notes are organised as follows. The remainder of this introduction contains a
brief historical perspective of the field of open-system dynamics in synthetic quantum matter,
and a discussion of the conceptual differences between considering equilibrium vs. driven-
dissipative systems. In Sec. 2 we present a brief overview of the experimental platforms used
to realise synthetic quantum systems. In Sec. 3 we introduce the general frameworks that are
currently employed to study the open-system dynamics of these setups: the Lindblad dynamics
for the density matrix, the Schwinger-Keldysh path integral, and quantum trajectories. Sec-
tion 4 summarises the key theoretical approaches to determine the dynamics in these various
frameworks.

We then proceed by discussing the key properties of open many-body systems. In Sec. 5
we review the applications of dissipative state preparation in the many-body context, where
a tailored design of dissipative processes and many-body jump operators can stabilise states
with non-trivial properties. As one system parameter is tuned, these states can undergo sharp,
non-analytic changes in their observable properties, leading to dissipative phase transitions, that we discuss in Sec. 6. We review their general classification, their relation with classical and quantum phase transitions and symmetry breaking and provide several examples of recent theoretical and experimental interest. Dissipative phase transitions can be also associated to the absence of a steady-state and to time-translation symmetry breaking, such as in boundary or dissipative time crystals. The focus of Sec. 7 shifts from the steady state to the dynamical transient approach to it. This is particularly relevant for many-body dephasing problems leading to heating dynamics or many-body losses associated to the emergence of the quantum Zeno effect. Further insight on open quantum many-body systems can be obtained by following the dynamics along quantum trajectories and the unravelling of the Lindblad master equation. The case of monitored dynamics unveils new phenomenology that is invisible to the average state, i.e. the density matrix. In Sec. 8 we will consider the dynamics under different unravelling corresponding to continuous weak monitoring and review the emergence of measurement-induced entanglement transitions. Finally, in Sec. 9 we will draw our conclusions. For completeness, Appendix A presents a self-contained derivation of Lindblad master equations.

Finally, these lecture notes should not be intended as a comprehensive review, but rather as a pedagogical introduction to the field. Many of the topics discussed in these lectures notes were also covered in the past in several, interesting reviews. The reader is encourage to read Refs. [28–33] to acquire a more comprehensive picture of the physics of many-body open systems.

## 1.2 Historical perspective

It is useful to frame the field of open quantum systems in the broader perspective of dissipative quantum systems that have been investigated for a long time. We focus our discussion here on Josephson junction arrays (see e.g. the review [34]), which were one of the first successful attempts to observe dissipative phase transitions in synthetic (i.e. artificially realised) systems. There are also a number of other topics where there was early work that could be understood in this context. This includes work on the analogy between lasers and phase transitions [35–37]—a topic we discuss further in Sec. 2.6 in relation to polaritons and photons. In addition a great deal of work has been done on phase transitions of two-level systems coupled to external baths, as reviewed in Refs. [14, 38].

Turning to Josephson junction arrays, early works initially prompted by proposals of P. W. Anderson [39] and A. J. Leggett [40] identified small Josephson junctions as very promising candidates to observe macroscopic quantum phenomena. The original motivation was mainly driven by the possibility to explore macroscopic quantum dynamics and the classical-to-quantum boundary in meso- or nano-structures. Parallel to these studies, Josephson junction arrays (regular lattices of small superconducting islands connected by tunnel junctions) were identified as promising platforms to study classical and quantum statistical mechanics models. The first experimental evidence of Berezinskii–Kosterlitz–Thouless transitions was observed [41–43], and in the late 80s a superconductor-insulator quantum phase transition was seen in quantum arrays using devices with small enough geometrical capacitance in order to have charging energies comparable to Josephson couplings [44]. Furthermore, several theoretical works pointed out that by including an external environment à la Caldeira–Leggett the zero- and finite-temperature phase diagram would be considerably enriched. In particular, on increasing the strength of dissipation the whole array would undergo a transition to a phase of global coherence solely controlled by dissipation, as in [45]. For further details, see the review [34] and references cited there.

The experimental verification of a transition induced by dissipation came much later from the independent work of two research groups. The group of J. Clarke in Berkeley [46] realised
<figure><img src="image_4.png" /><figcaption>Figure 1: (a-b) Experimental measurement of a dissipative phase transition in a Josephson junction array, from [46] [Copyright (1997) by the American Physical Society]. (a) The Josephson array (on the top of the device) is in electrostatic contact with a two-dimensional electron gas (2DEG) that is responsible for a voltage noise on the superconducting islands. The level of the noise is modulated by changing the resistance $R_g$ of the 2DEG. (b) The current-voltage characteristics for different levels of the noise, i.e. different values of the resistance $R_g$. On increasing $R_g$, i.e. decreasing the dissipation, the current-voltage characteristics change from superconducting-like (finite current at zero voltage) to insulator-like (zero current below a threshold voltage). (c) Phase diagram—sketch of the zero-temperature phase diagram, above a certain value of the dissipation strength the array acquires global coherence for any value of the Josephson coupling [The phase diagram is sketched from [45]].</figcaption></figure>

the Josephson array on top of a 2D electron gas, controlling the strength of damping through a gate voltage that allowed to vary the conductance of the 2DEG. Y. Takahide and coworkers [47], in the group of S. Kobayashi in Tokyo, instead engineered shunt resistances in parallel to each junction in a one-dimensional array (they had to fabricate several different samples with different shunts in order to observe the increasing effects of dissipation!). In Fig. 1 we show the results of the Berkeley experiment. The experimental setup is shown in panel (a): dissipation is controlled by coupling the arrays to a two-dimensional electron gas and modulating its conductivity (hence the voltage noise acting on the array) by means of the back gate voltage. Panel (b) shows the current-voltage characteristics for different level of dissipation (i.e. for different values of the resistance of the 2DEG $R_g$). At low $R_g$ the current-voltage characteristics are that of a superconductor (with global phase coherence). On increasing $R_g$ there is a qualitative change to the current-voltage characteristics of an insulator (Coulomb blockade regime) with zero current below a threshold voltage. This is, to our knowledge, the first experimental observation of a dissipative phase transition in a quantum many-body synthetic system.
In the language of quantum information processing, Josephson junction arrays can be considered as the ancestors of synthetic systems, although with limited control on their quantum evolution due to the presence of imperfections and noise. The impressive technical advances over the last decades have enabled the design of new platforms to study open many-body systems. As we will see in the forthcoming Sections, these new systems offer an unprecedented control both in the realisation of synthetic quantum matter and in the study of its dynamical evolution (for instance, the coherent dynamics has always proven impossible to study in Josephson arrays).

An important feature characterising the early works and differentiating them from the most recent experiments, is that Josephson arrays operated in a regime where the detailed balance relation is always satisfied; as we discuss in Sec. 1.3, this is a hallmark of the global thermal equilibrium of the system and of its environment. As already mentioned, an excellent description of these systems could be based on Caldeira-Leggett models, provided they are properly generalised to the many-body context. Furthermore, in many cases the out-of-equilibrium properties of the superconducting order parameter of Josephson junctions would be described by dynamical equations typical of classical dynamical critical phenomena [48]. Instead, in the modern platforms that are currently being developed (funnily enough, some of them still use superconductors and Josephson junctions), the dynamics is governed by the external driving and by dissipation, so that the state cannot be found by considering an equilibrium partition function—not even including the environment within the partition sum. We discuss these points further in Sec. 1.3 below. In the Sec. 2 we will review many experimental platforms of more recent interest which are currently playing a major role in the study of synthetic quantum matter and which enable the study of open quantum many-body physics in controlled settings.

## 1.3 Relation with thermal equilibrium

In this section we discuss some key concepts about the distinction between quantum systems in equilibrium and driven-dissipative systems. Note that we use the term driven-dissipative here to identify systems where there is either some driving term in the Hamiltonian along with coupling to an environment or coupling to multiple environments. The term “open quantum systems” is broader, and can include the case of a system coupled to a single environment. One point we discuss in particular here is to expand on the discussion in the previous section, distinguishing the behaviour of a system coupled to a single environment from the behaviour of a driven-dissipative system.

To set this discussion in context we first provide a brief summary of the key ideas that arise in equilibrium statistical physics, in order to then comment on how these are changed in driven-dissipative systems. For a full discussion, the reader may see textbooks on thermodynamics and statistical physics, e.g. [49, 50]. Equilibrium statistical physics starts from the postulate of equal equilibrium probabilities (PEEP), which states that all microstates occur with equal probabilities. This then leads to the idea that the probability of a macrostate—where some macroscopic variables are specified—is determined by the number of ways in which such a macrostate can arise, i.e. the number of microstates corresponding to a given macrostate. Defining the entropy as the logarithm of this number of microstates then leads to the notion that the most probable macrostates are those that maximise entropy.

To get from the PEEP to the familiar Boltzmann distribution, one may consider a system that is weakly coupled to an environment in thermal equilibrium. In particular, we consider the case where the system plus environment is in a microcanonical state, i.e. the total energy of the system plus environment is fixed. Weak coupling allows one to consider the energy of the system and environment additively, i.e. to ignore any energy associated with their coupling. Thus, the energy in the environment is determined by the (fixed) total energy and the energy
of the system. If the system exchanges only energy with the environment (as opposed to e.g. particles, or any other conserved quantity), then one can consider the entropy of the environment as being a function of the energy of the system. Since the environment can be assumed to be large, and thus weakly affected by this transfer of energy, one can consider the entropy of the environment to depend linearly on the energy of the system. The coefficient in this expression is the inverse temperature. This then leads to the standard Boltzmann distribution: $P \propto \exp(-\beta(E - TS_S))$ with $\beta = 1/(k_B T)$ and $S_S$ being the entropy of the system.

To understand what changes in driven-dissipative systems we highlight a key assumption to obtain the Boltzmann distribution, namely the existence of a well defined entropy cost to increasing the energy of the system. This concept applies for coupling to a single environment, but cannot apply when a system is coupled to multiple environments that are not in thermal equilibrium with each other. For multiple environments with multiple temperatures, the entropy cost of increasing the system energy depends on which environment is involved. As such, for systems coupled to multiple environments, the relative coupling rates and dynamics of the system will matter, and the probability that the system is in a given state cannot only depend on its energy.

While we have introduced the Boltzmann distribution from the assumption of equilibrium with an external environment, the same expression can also be derived for a large system: here, when considering one part of the system, one can regard the rest as acting as an effective environment. In a finite system there will be a distinction between the result derived above, which are correct for the canonical ensemble (ensemble at fixed temperature) vs the microcanonical ensemble (ensemble at fixed system energy). However, in the thermodynamic limit, these ensembles become equivalent for extensive quantities.

In the discussion so far we assumed weak coupling between the system and its environment. As long as one is restricted to a single environment (or multiple environments in thermal equilibrium), one can extend the above equilibrium analysis further. One may thus discuss how behaviour changes when a system is strongly coupled to the environment, as is of relevance to the discussion in the previous section. We discuss first an approach that works for static equilibrium properties, and then discuss dynamical responses.

For static properties, one may consider the Boltzmann distribution for the combined system and environment. This means the system density matrix is given by the combined Gibbs state:

$$ \hat{\rho}_S = \frac{1}{Z} \mathrm{Tr}_E \left[ \exp(-\beta(\hat{H}_S + \hat{H}_E + \hat{H}_{S-E})) \right]. \qquad (1) $$

Here we have written the equilibrium system density matrix, $\rho_S$, in terms of the Hamiltonian of the system, $\hat{H}_S$, the environment, $\hat{H}_E$, and the system-environment interaction, $\hat{H}_{S-E}$. This expression generalises the Boltzmann distribution discussed above, by assuming internal equilibrium of the combined system and environment. The partition function $Z$ appears here as normalisation, required to ensure $\mathrm{Tr}[\rho_S] = 1$.

If the coupling to the environment is weak, then $\hat{H}_{S-E}$ can be neglected. In this case, after tracing over the environment, one recovers $\hat{\rho}_S \propto \exp(-\beta\hat{H}_S)$. Assuming $\hat{H}_{S-E}$ vanishes is reasonable when considering the equilibrium state of the system. However, one should note that if the system environment coupling vanishes, the timescale for the system and environment to come into equilibrium will diverge. The equilibrium state can generally be assumed to be equivalent to the long time average—this is the notion of ergodicity. Formally this equilibrium ensemble with the system Hamiltonian corresponds to first taking the limit of long times, and only then considering the limit of vanishing system-environment coupling.

If the coupling to the environment is not weak, then the equilibrium state defined by Eq. (1) will not equal $\exp(-\beta\hat{H}_S)/Z_S$ (known as the system Gibbs state), but will be modified by the coupling to the environment. Since an equilibrium density matrix, $\hat{\rho}_S$, does however exist
(as defined by Eq. (1)), one can however define an effective “Hamiltonian of Mean-Force” (HMF) [51, 52] by writing:

$$ \hat{\rho}_S = \frac{e^{-\beta \hat{H}_{\text{HMF}}}}{\mathcal{Z}_{\text{HMF}}}. \qquad (2) $$

The state $\rho_S$ as written here is sometimes referred to as the Mean-Force Gibbs state; as discussed here, this is equivalent to the combined Gibbs state. Mathematically the definition of $\hat{H}_{\text{HMF}}$ just corresponds to extracting the logarithm of the equilibrium density matrix. It is easy to verify, by using the identity

$$ e^{-\beta(\hat{H}_S + \hat{H}_E + \hat{H}_{S-E})} = e^{-\beta(\hat{H}_S + \hat{H}_E)} T_{\tau} e^{-\int_0^{\beta} d\tau \hat{H}_{S-E}(\tau)}, $$

with the coupling Hamiltonian expressed in the time-ordered interaction representation, that the HMF

$$ \hat{H}_{\text{HMF}} = \hat{H}_S - \frac{1}{\beta} \ln \langle T_{\tau} e^{-\int_0^{\beta} d\tau \hat{H}_{S-E}(\tau)} \rangle_E. $$

The form of $\hat{H}_{\text{HMF}}$ may be complicated with environment-mediated interactions, and thus potentially nonlocal, and its form may depend on temperature. Nonetheless, this formulation identifies a natural generalised of the concept of potential of mean force [51, 52] sometimes studied in classical statistical physics. There can also exist simple forms that arise in the limit of strong system-environment coupling [53]. In many situations, it is also convenient to use a path integral formulation and consider the effective action deriving after tracing out the environmental degrees of freedom [38]. Most of the early works studying dissipative phase transitions in Josephson arrays, see Section 2.5, approached the problem in this way.

While the above suggests an effective mean-force Hamiltonian can be found to determine static properties, the situation is more complicated when considering dynamic properties. No effective system Hamiltonian can correctly capture the dynamical response of a system strongly coupled to an environment. This is because the timescales of the dynamics of the environment become important. One can model the effects of the environment through an effective Keldysh action [54] (a concept that we discuss further below). The Keldysh action however describes a frequency-dependent response of the environment, and this cannot be described by any effective Hamiltonian. It is this physics that led to the interest in dissipative phase transitions discussed in the historical perspective above.

One question that arises when considering systems strongly coupled to their environment is how one can tell whether the system is in thermal equilibrium or not. As noted above, the Hamiltonian of mean force can always be defined from the logarithm of the density matrix, and since $\hat{H}_{\text{HMF}} \neq \hat{H}_S$, one cannot simply use the equilibrium density matrix to practically determine whether a given system is in equilibrium or not. What can however be used is the fluctuation dissipation theorem [55, 56] (FDT). The FDT is defined by considering the two-time correlation functions of some observable $\hat{A}(t)$. Specifically, one considers the power spectral density $S_A(\omega)$ and the response function $\chi_A(\omega)$, that are the Fourier transforms of the following two-time correlation functions:

$$ S_A(\tau) = \frac{1}{2} \langle \{\hat{A}(t + \tau), \hat{A}(t)\} \rangle, \quad \chi_A(\tau) = i\Theta(\tau) \langle [\hat{A}(t + \tau), \hat{A}(t)] \rangle, \qquad (3) $$

where {} indicates the anticommutator, [, ] indicates the commutator, and $\Theta(\tau)$ the step function. Note that in a stationary state, these expressions depend only on the time difference, $\tau$ and not the initial time $t$, hence the possibility of Fourier transforming only with respect to $\tau$. What the FDT proves [55, 56] is that in thermal equilibrium, the Fourier transforms of these two functions obey:

$$ \frac{S_A(\omega)}{\text{Im}[\chi_A(\omega)]} = \coth(\beta \omega). \qquad (4) $$
We have assumed here the operator $\hat{A}$ corresponds to an observable, so $\hat{A}$ is Hermitian and cannot be a fermionic operator. Beyond this restriction, the expression is general: it holds for any observable $\hat{A}$ and for any form of system-environment coupling. One can understand that this expression works because $\chi_A(\omega)$ measures the density of states of the system, while $S_A(\omega)$ measures the noise due to thermally occupied modes. As a result, the combination given here allows the density of states (and thus effects of strong system-environment coupling) to "cancel out" from the expression.

Summarising, in systems coupled to a single environment, no matter how strongly, one can expect the FDT to hold in full at late enough times, i.e. once the system has reached equilibrium. In general open quantum systems, such as driven-dissipative setups, the FDT may approximately hold over some range of frequencies, but it will not do so in general, since such systems are not in thermal equilibrium. Other theoretical frameworks then need to be developed, and the Lindblad master equation is the simplest and most widespread that describes a system reaching a stationary state that is, in general, not at thermal equilibrium. The study of the varied phenomenology that can appear in these setups is the object of these lecture notes.
1.4 Notation

It is convenient to summarise in a table the acronyms and the notation principally used throughout these lecture notes; some exceptions exist and we highlight these when they occur.
We set $\hbar = k_B = 1$ throughout.

<table><tr><td>Asymptotic Decay Rate</td><td>ADR</td></tr><tr><td>Dynamical Mean-Field Theory</td><td>DMFT</td></tr><tr><td>Dissipative Phase Transition</td><td>DPT</td></tr><tr><td>Mean-Field</td><td>MF</td></tr><tr><td>Steady State</td><td>ss</td></tr></table>

<table><tr><td>Commutator</td><td>[A, B] = A B - B A</td></tr><tr><td>Anticommutator</td><td>{A, B} = A B + B A</td></tr><tr><td>Lattice sites</td><td>i, j</td></tr><tr><td>Coordination number</td><td>z</td></tr><tr><td>Space coordinate in continuous systems</td><td>r</td></tr><tr><td>Time</td><td>t</td></tr><tr><td>Number of particles</td><td>N</td></tr><tr><td>Number of lattice sites</td><td>N<sub>s</sub></td></tr><tr><td>Length of system</td><td>L</td></tr><tr><td>System dimensionality</td><td>d</td></tr><tr><td>Volume</td><td>L<sup>d</sup></td></tr><tr><td>System dimensionality</td><td>d<sub>H</sub></td></tr><tr><td>Cartesian components</td><td>α, β = x, y, z</td></tr><tr><td>Wave-vector in continuous systems</td><td>k, q</td></tr><tr><td>Labels (not referring to space/lattice)</td><td>Greek letters δ, ..., μ, ν, ...</td></tr><tr><td>Operators and Superoperators</td><td>Ô, ˆÔ</td></tr><tr><td>Density matrix of the system</td><td>ρ̂(t)</td></tr><tr><td>Hamiltonian of the system</td><td>Ĥ</td></tr><tr><td>Lindblad jump operators</td><td>L̂</td></tr><tr><td>Lindbladian superoperator</td><td>L̃</td></tr><tr><td>μ-th eigenvalue of the Lindbladian superoperator</td><td>λ<sub>μ</sub></td></tr><tr><td>non-Hermitian Hamiltonian</td><td>Ĥ<sub>nH</sub></td></tr><tr><td>Bose operators</td><td>â, â†</td></tr><tr><td>Fermi operators</td><td>ĉ, ĉ†</td></tr><tr><td>Spin operators on a lattice</td><td>σ̂<sup>α</sup><sub>i</sub></td></tr><tr><td>Field operators</td><td>ψ̂</td></tr><tr><td>Magnetic interaction</td><td>J<sup>α</sup><sub>i,j</sub></td></tr><tr><td>On-site interaction</td><td>U</td></tr><tr><td>Hopping coupling constant</td><td>t<sub>H</sub></td></tr><tr><td>Light-matter interaction</td><td>g</td></tr><tr><td>Detuning</td><td>Δ</td></tr><tr><td>Dissipative rates</td><td>κ, γ, Γ, ...</td></tr><tr><td>Asymptotic Decay Rate</td><td>λ<sub>ADR</sub></td></tr><tr><td>Control parameter for phase transition</td><td>η</td></tr></table>
## 2 Experimental platforms

In this Section we briefly discuss some of the experimental setups which are relevant for the field of many-body open quantum systems. In line with the general goal of these lecture notes, we focus on synthetic quantum matter where many-body effects and open-system dynamics are particularly relevant: superconducting circuits and circuit QED arrays [28,57], ultracold atoms with controlled dissipation [58], ultracold atoms in high-finesse cavities [59, 60], and trapped ions [61, 62]. These experimental platforms are often presented as paradigmatic examples of isolated quantum systems. However, in most experimental realisations there will be effects due to the presence of an environment, typically represented by the experimental apparatus itself. Although the effect of such an environment can often be made rather small, it is never completely negligible. On the other hand, what is very important here is that it can be made highly controllable.

We emphasise that there are several reviews already existing on the related topic of quantum simulation, such as Ref. [12, 13, 63], therefore we will limit ourselves to discuss the key aspects that are particular relevant in the present context.

### 2.1 Trapped ions

Single trapped ions represent elementary quantum systems with degrees of freedom that are well isolated from the environment [61,64]. Experimentalists working with trapped ions have for long time manipulated the quantum state of a single ion using coherent fields, as well as using the opportunities offered by spontaneous emission (see e.g. Refs. [65-67]), which is one of the simplest open-quantum-system effects that can be considered. Using this, it is possible to bring the ions nearly to rest by laser sideband cooling, and to prepare the ion in a desired electronic quantum state by the use of optical pumping. [64]. This makes trapped ions well-suited for the study of open-quantum-system dynamics under well-controlled conditions.

As a simple example, we briefly summarise the key results reported in Ref. [65], which achieved the first experimental observation of the quantum Zeno effect [68], i.e. strong coupling to an environment inhibiting coherent transitions between the energy levels of a quantum system. Working with $Be^+$ ions, this work considered two hyperfine levels coupled by a radio-frequency field while a strong laser pulse coupled one of those levels to an excited unstable state which would undergo spontaneous emission (see the sketch in Fig. 2, left panel). By observing the spontaneously emitted photons, it was possible to quantitatively measure the reduction of coherent transitions between the hyperfine states.

The early results on the manipulation of single ions can be seen as foundational for contemporary research. Indeed, in quantum optics and atomic physics, the techniques of optical pumping and laser cooling are now routinely employed for the dissipative manipulation of quantum states, although on a single-particle level. The current studies on engineering dissipative dynamics that we discuss in Sec. 5 can be seen as a generalisation of these techniques to a many-particle context: it is not by accident that the quantum Zeno effect is currently experiencing a revival in the many-body setting, as we discuss more extensively in Sec. 7.2.

In more recent years, trapped ions have started to play a significant role in the context where the role of the environment is played by a measurement apparatus and the quantum system is *monitored*. The possibility of repeatedly and frequently measuring the ions—which is at the heart of the quantum Zeno effect—is also what places these setups on the forefront of this research field. As an example, in Ref. [69] the authors consider one-dimensional chains of ions of length spanning from 4 to 8 ions and engineer a unitary dynamics (technically, a random circuit composed of unitary gates) with interspersed projective measurements. When the measurement rate is low, the evolution reaches a mixed phase: the multiple measurement outcomes create a mixed density matrix. On the contrary, for a measurement rate exceeding
<figure><img src="image_5.png" /><figcaption>Figure 2: (Left) Energy levels of $^9$Be$^+$ ions in a magnetic field B employed in the experiment on the quantum Zeno effect. The two hyperfine levels, 1 and 2, are coupled by a radio-frequency field; a strong laser pulse couples the level 1 with the excited level 3, which undergoes spontaneous emission. In this way, the quantum mechanical properties of a single ion are manipulated using coherent and dissipative processes. From Ref. [65] [Copyright (1989) of the American Physical Society]. (Right) Experimental tools for the simulation of open quantum systems with ions. n + 1 ions are used to simulate the open-system dynamics of n ions (the system, S) in a controlled fashion. Spontaneous emission is induced in a controlled fashion on the first ion (the environment, E), inducing the desired open-system-dynamics on the other n ions. Figure based on experiment in Ref. [67].</figcaption></figure>

a given threshold, the evolution reaches a pure phase: measurements project the system on to a deterministic quantum state. We note that other experimental platforms have also been used to explore measurement-induced physics, including superconducting circuits and ultra-cold atoms, even though trapped ions are particularly well suited due to quality and speed of measurements. We will discuss the theoretical and experimental considerations associated to monitored quantum systems in Sec. 8.

Regarding the development of experimental platforms, it is interesting to mention one early approach to creating a many-body synthetic quantum system based on trapped ions employing the coupling to an environment [67]. Building on a well-tested quantum computing architecture, which is designed to undergo purely unitary (i.e. closed-system) dynamics, this work combined multi-qubit unitary gates with optical pumping to implement dissipative processes. In particular, the experimental platform is composed of an array of $n+1$ ions, with the first ion serving as an ancilla, and the remaining $n$ ions being the system (see the sketch in Fig. 2, right panel). The ancilla is coupled to the system through application of specific gates, but also undergoes dissipative dynamics ultimately due to spontaneous emission events. This in turn induces a controlled and tunable dissipative dynamics on the rest of the system. As an example, this platform was employed to prepare entangled Bell states composed of $n=2$ ions or GHZ states composed of $n=4$ ions. In the first case, the environment acts on two ions at the same time, and its effect depends on the reduced density matrix of the two particles. In the second case, the coupling to the environment is mediated by four ions. Given the relatively small number of qubits involved, it was possible to perform full quantum state tomography of the system and demonstrate that, thanks to the dissipative processes, the desired state had been realised.

In more recent years, several other experiments have been reported which have actively exploited the judicious control of dissipation for the study of quantum matter or dynamics;
we mention here a few of them. One example is the phonon laser realised in [70] where the phase diagram and the dissipative transition have been detected. More recently the dissipative generation of an entangled state of two trapped-ion qubits, based on a scheme developed in Ref. [71], has been experimentally realised in Ref. [72]. The possibility of measuring and reset- ting the ions with high fidelity and addressability is crucial in the results reported in Ref. [73], where entanglement entropy of a spin chain is directly measured. Finally, the possibility of using dissipation to simulate the quantum dynamics of electronic transfer, a process at the heart of several biochemical processes, has been reported in Ref. [74].

## 2.2 Ultracold atoms

In recent years, experiments with ultracold atoms have evolved towards the development of platforms where the coupling to an environment can be precisely characterised and tuned, so that it can be reduced or enhanced over several orders of magnitude; consequently, open- quantum-system dynamics can be studied with these platforms in the most genuine spirit of quantum simulation. One key example of this is cold atoms coupled to optical cavities, where the light in the cavity represents an open quantum system. In the absence of coupling to an optical cavity, there are two main open-quantum-system effects that can take place: losses of particles from the trapping potential, and heating due to the interaction with the electromag- netic fields that trap or control the gas. In the following we first present these experiments and then discuss in a separate Section the case of cavity QED with ultracold atoms. Finally, we will introduce Rydberg states of ultracold atoms.

**Atom loss.** Losses have been a major concern from the early days of Bose-Einstein conden- sation. Several articles dating from the 1990's and early 2000's have characterised the dy- namics of the number of bosons in the context of evaporative cooling and for Bose-Einstein condensates, proposing simple theories that are mean-field in spirit as they neglect correla- tions among the atoms [75–77]. Generally, there are three kinds of processes that can give rise to atom loss, which are classified according to the number of condensate particles that are involved. One-body losses are due to the interaction of condensate particles with thermal background particles that float in the vacuum chamber [78]. Two-body losses typically char- acterise experiments with atoms trapped in a metastable excited energy level or, in general, with molecules [79–82]. Other situations where the two-body loss process is relevant have also been reported in the literature, for instance with magnetic atoms like Chromium, or with alkali atoms trapped in the state of largest hyperfine spin—here this process is typically irrele- vant for Rubidium but important for Caesium. In such cases, collisions between two particles can induce transitions towards internal levels with lower energy: the internal energy that is gained is converted into kinetic energy of the particles, that can then escape the trap. This contrasts to atoms in their lowest internal state; in such cases, even if a lower-energy bound state exists, the combined effects of energy and momentum conservation typically prevent re- laxation to a bound state. Three-body losses can however always occur [76, 80, 82, 83], and are the result of recombination processes where two atoms bind into a molecule and transfer part of their binding energy to a third atom in the form of kinetic energy, so that both the molecule and the atom can escape the trap. In recent years, there has been work focused on the engi- neering of loss processes, in order to have a higher degree of control, focusing in particular on one-body [84, 85] and two-body [86, 87] loss processes. A series of experimental works has also focused on the possibility of inducing losses in a localised region of space, typically a single site of an optical lattice [88]. Some aspects of the intriguing physics that can arise from loss processes will be discussed in details later in the review.

As a paradigmatic example, let us discuss the experiment reported in Ref. [89], where a bosonic Yb isotope is cooled in a three-dimensional optical lattice. Two-body inelastic atom
<figure><img src="image_6.png" /><figcaption>Figure 3: Time-of-flight images taken with different lattice depths, from the very deep case of 20 recoil energy $E_R$ to the case of no lattice $0 E_R$; the first line corresponds to an atomic level that is not subject to dissipation, the second instead corresponds to a metastable energy level that is subject to strong two-body losses; the transition to a superfluid becomes less clear. From Ref. [91] [Copyright (2019) by the American Physical Society].</figcaption></figure>

losses at a controllable rate are induced via a single-photon association process: two atoms in the same lattice site are photoassociated into a molecular state and immediately dissociated. This process releases a high kinetic energy to the dissociated atoms, which eventually escape the optical lattice. By changing the intensity of the photoassociating laser, the rate of two-body losses is controlled with high fidelity. This then allows one to study the effect of two-body losses on the zero-temperature quantum phase transition between a Mott insulator and a superfluid, first reported in Ref. [90].

Figure 3 is taken from Ref. [91], where the authors present a series of time-of-flight absorption images obtained by changing the final lattice depth from the Mott insulator to the superfluid regime for two atomic energy levels characterised by a significantly different intrinsic two-body loss rate. For the dissipative energy level, the interference pattern that characterises the superfluid phase in the shallow lattice regime becomes blurred.

References [89, 91] opened a new avenue for quantitative studies of the depletion of coherence in an atomic gas, as characterised by the formation of off-diagonal long-range order due to two-body loss events. We mention here Ref. [92] where a theoretical semiclassical description of the dissipative dynamics is proposed.

In Secs. 5.3.3 and 7.2 we will present a more thorough discussion of the theoretical and experimental developments associated to the open-quantum-system dynamics induced by losses.

**Atom heating.** Heating processes are the second open-quantum-system process that needs to be characterised, primarily because they drive the system towards featureless high-temperature states and wash out all interesting quantum effects. There are many microscopic processes that can give rise to heating effects. These include fluctuations of the laser field responsible for the trapping of the particles in an optical lattice, which can induce a noise dynamics of the trapping potential; collisions among the atoms, particularly those that release energy from excited internal states (the difference with loss processes is that here, after the collision, the reaction products are still trapped); incoherent scattering of the lattice light. The study of the heating dynamics due to incoherent light scattering for single atom systems is well understood [93,94]. A key question of current interest is the study of how this heating dynamics occurs for many-body quantum states, where the interplay of few body heating processes and many-body physics can give rise to interesting effects [95–97]. Related to that question is the study of the effect of making frequent measurements of the position of the atoms in optical lattices [98, 99]. Similar to the discussion of losses given above, there have recently been experiments where the heating dynamics could be studied in a controlled way. This can
<figure><img src="image_7.png" /><figcaption>Figure 4: (a) An ultracold gas of bosonic $^{174}$Yb atoms is loaded in a two-dimensional array of one-dimensional systems. The green arrow represents a laser propagating through the systems and inducing spontaneous emission in a controlled fashion. (b) Details on the energy levels and laser couplings. The orange arrow represents a laser close to an atomic resonance that induces spontaneous emission at a rate $\Gamma \sim 520s^{-1}$. The atom is randomly recoiled and this destroys the quantum coherence of the gas. Figure inspired by Ref. [100].</figcaption></figure>

be done by illuminating the atoms with laser light that is resonant with an atomic transition, so that spontaneous emission can take place at a controlled rate.

As an example, we consider the experiments reported in Refs. [100, 101]. A bosonic Yb isotope is cooled in a stack of two-dimensional optical lattices, and exposed to a laser tuned close to an atomic resonance which induces spontaneous emission at a controllable rate. The spontaneous emission of photons induces a random recoil of the atom in the trap, that is responsible for the loss of coherence (see the sketch in the left panel of Fig. 4). The experiment shows the standard time-of-flight picture with coherent Bragg peaks that are washed out in time by the induced spontaneous emission processes. In this controlled experiment, the appearance of an anomalous decay time is also reported, and will be discussed in Sec. 7.1.

## 2.3 Cold atoms in optical cavities

As noted above, one particular case where many-body open quantum behaviour arises is with ultracold gases trapped in optical cavities. While there can be experiments where one uses the bare coupling between cold atoms and a cavity light mode [102], these coupling strengths are not tunable and it is hard to realise strong coupling. Nevertheless, this setting has been used for example to generate squeezing of the collective spin of an atomic ensemble by coupling to a cavity resonator [103, 104]. One can however use an alternate scheme, as discussed next, to realise a tunable form of matter-light interaction, which has led to a growing research field (see e.g. the reviews [59, 60]).

The key idea, proposed in Ref. [106], is to realise coupling between atoms and cavity modes via a two photon process, involving one cavity photon and one pump photon (see Fig. 5(a)). This idea has several advantages. Firstly, it allows control over the effective light-matter coupling strength, by tuning the laser pump strength. An additional benefit is that this scheme overcomes any “No-go theorem” associated with bounds on the strength of matter-light coupling compared to energy splitting [107]; this is because the coupling to light is controlled
<figure><img src="image_8.png" /><figcaption>Figure 5: (a) Cartoon showing self organisation of atoms in an optical cavity, with transverse pumping. (b) Energy levels for Raman pumping scheme. Transitions between two low-lying atomic states are driven by a two photon process, involving scattering light between the pump and cavity. (c) Theoretical phase diagram of self organisation of atoms in an optical cavity, predicted by solving the Gross-Pitaevskii equation for atoms, coupled to a self-consistent cavity field. The coupling (horizontal axis) is controlled by the pump strength. The colour scale indicates the intensity of light in the cavity. Figure inspired by Ref. [105].</figcaption></figure>

by the dipole matrix element of the virtual transition, rather than that of the transition between the two low-lying states. As a result it is possible to realise the Hepp–Lieb–Dicke superradiance transition [108, 109], where the cavity acquires a macroscopic photon field above a critical coupling strength. Such a transition was indeed observed in a realisation where the two low-lying atomic states correspond to different motional states of the atoms [110] (see Fig. 5(b)).

The Raman pumping scheme in fact permits even greater control over the effective Hamiltonian, beyond controlling the strength of matter-light coupling. In the rotating frame of the pump, the effective cavity frequency becomes the energy cost of scattering a photon from the pump to the cavity, so is tunable. Further, if one uses two separate pump beams for the two legs of the Raman transition, the effective energy splitting of the low-lying levels becomes the two-photon detuning [106], allowing control over this via the beat frequency between the two pump beams.

When combined with cavity geometries that permit control over multiple separate cavity modes [111], this Raman pumping scheme allows for considerable control. For example, experiments with two competing cavity modes have been able to realise a model with a continuous symmetry (associated with relative amplitude of the two cavity modes), and thus to explore spontaneous symmetry breaking of a continuous symmetry [112, 113]. When going from two cavity modes to many, further new features arise. By using a confocal cavity—a cavity where sets of transverse modes come in degenerate families—one can build transversely localised wavepackets from the different degenerate modes. By slightly adjusting the mirror spacing away from the confocal degeneracy point, one can then realise tunable-range interactions between atoms [114]. Recent work has further shown how the two ideas of short-range interactions and continuous symmetry breaking can be combined to realise a compliant optical lattice, that can distort, and thus has sound modes [115] (see Fig. 6). Several other interesting directions have been explored with similar experiments. By using positioning of atomic ensembles along the cavity axis and using a (nearly) concentric cavity, one can engineer programmable interactions between clumps [116], or one can use this to explore chaotic dy-
<figure><img src="image_9.png" /><figcaption>Figure 6: Cartoon of atomic gas in a multimode cavity. Self organisation into a density-wave pattern occurs, but because of the freedom allowed by the multimode cavity, the phase of the density wave can vary across the cavity. To allow full phase freedom for the density wave, a pump scheme with two colours of light is required (see Ref. [121] for details. [Image credit B. Marsh, Y. Guo. Based on experiment in Ref. [115]]).</figcaption></figure>

namics [117]. There are also proposals for how an associative memory and spin-glass physics
can be realised with a confocal cavity [118–120].

In the context of a single mode cavity, when the atomic states involved are different mo-
tional states, the transition to a superradiant state is associated with the formation of an optical
lattice. As such, there can also be a Mott insulator transition that occurs as the depth of this
emergent lattice grows [122]. Very recently, there have also been work realising such spatial
ordering using a degenerate Fermi gas [123].

As noted above, such experiments are intrinsically open quantum systems, due to the
combination of the cavity loss, and the driving via the transverse pump. While many of the
phenomena mentioned above (self organisation transitions in various geometries) can be un-
derstood equally in a closed system, there have also been clear experimental signatures that
can only be understood by including dissipation. These include the observation of limit cy-
cles [124, 125]. Much of the theoretical work on such systems—particularly those consider-
ing Bose-Einstein condensates in single-mode cavities—has been able to make use of mean-
field treatments. However, there are specific questions—particularly when considering cor-
related states of atoms, such as Mott insulating states, where beyond-mean field approaches
become necessary. In this context various approaches have been used, including matrix prod-
uct states [126] and multi-configuration time-dependent Hartree (MCTDH) [127]. Beyond
mean-field methods also become necessary when considering spatially extended systems in
nearly confocal cavities, where fluctuations can significantly modify phase boundaries [128].

The design of new platforms may lead to the realisation of different forms of couplings to
the bath that, while retaining the Lindblad form, involve more complex (long-range and/or
correlated) jumps and feedback. A dissipative term in the Lindbladian with power-law-
decaying correlations could in principle be realised in cold atoms inside a cavity. The scheme,
proposed theoretically in [129, 130], requires three-level atoms in the presence of a magnetic
field gradient, and a Raman beam. The flexibility in the design of suitable jump operators
means, in the context of quantum state preparation, the possibility to have access to a wider
range of out-of-equilibrium many-body states or phases.
<figure><img src="image_10.png" /><figcaption>Figure 7: (a) Illustration of Rydberg blockade and facilitation for two atoms. The transition frequency between the state with a single Rydberg excitation, $|r, g\rangle$ and the state with two excitations, $|r, r\rangle$ depends on the spatial separation between the atoms. As a result, for red-detuned driving there is a blockade, as the transition is shifted out of resonance. In contrast for blue-detuned driving there can be a specific separation at which the transition is resonant, leading to facilitated dynamics. (b,c,d) Bimodal counting statistics for an ensemble of Rydberg atoms. Panel (b) shows the experimental setup, with insets illustrating excitation processes for increasing pump frequency. Panels (c,d) show a measured phase diagram, found by measuring (c) the mean number of Rydberg excitations $\bar{n} = \langle n \rangle$, and (d) the Mandel Q parameter of the distribution of excitations $Q = \langle (n - \bar{n})^2 \rangle / \bar{n} - 1$. [Panels (b,c,d) adapted from Ref. [133] - Copyright (2014) by the American Physical Society].</figcaption></figure>

## 2.4 Rydberg gases

Another special class of experiments with cold atoms is that involving atoms in highly excited “Rydberg states”. When atoms are in highly excited states, this can significantly modify their optical properties as compared to ground state atoms. Most notably it leads to enhanced dipole-dipole interactions, due to the increased radius of the electronic orbit, and to enhanced lifetimes, as reviewed in Ref. [131].

Direct transitions from the atomic ground state to a Rydberg state are forbidden by selection rules, so that transitions are typically driven by a two-photon process, via driving with two laser frequencies. As discussed further below, one can also consider replacing one of these with an optical cavity mode. If the intermediate state is far off resonance, the fact that the process is a two-photon transition becomes unimportant, and one can consider an effective coherent driving of the ground state to Rydberg state. Because of the large dipole-dipole interaction between atoms in Rydberg states, the presence of one excited atom shifts the energy required to excite other atoms to the Rydberg state (see Fig. 7). When the drive frequencies are resonant with the optical transition, or detuned slightly below resonance, then the energy shift due to an excited atom moves the transition for a second atom further away from resonance. This leads to the “Rydberg blockade” [132], where the presence of one atom prevents others being excited over some range of distances. For a small atom cloud, this blockade means that only a single Rydberg excitation may be possible in the cloud: this allows one to gain the collective enhancement in matter-light coupling that comes from considering large numbers of atoms while still restricting to a two-level system of zero or one excitations. For larger clouds of atoms, multiple excitations will be allowed but separated by a blockade radius.

The ability to drive a single Rydberg transition in an ensemble of atoms provides a potential route to engineer many-body gates between atoms, and there have been proposals to use this for digital quantum simulation [134] of both closed-system and dissipative models. In addition to proposals for experiments that involve directly exciting a Rydberg state, there have been a number of proposals involving atoms in a controllable superposition of ground and
Rydberg states. This induces “dressing” of the atom properties by the Rydberg state, providing a controllable route to enhance atom-atom interactions [135–137].

Including dissipative processes in Rydberg experiments is straightforward, by considering the effect of spontaneous radiative decay from the excited atomic states, including the Rydberg state. This gives a natural source of dissipation that is counterbalanced by the external laser drive. Depending on the parameter range considered (e.g. for Rydberg dressed atoms as noted above), one may sometimes adiabatically eliminate excited states to yield an effective dissipative model within a ‘slow’ subspace. This has allowed for a wide range of theoretical proposals using dissipative dynamics with Rydberg atoms [138, 139]. In several of these examples the role of the Rydberg blockade is inverted. This can be done by detuning the driving laser to frequencies above the transition frequency. In such cases, the presence of one excited atom moves the transition frequency for a subsequent atom closer to resonance. Both regimes (Rydberg blockade and facilitated dynamics) are a natural context to study “Kinetically constrained models” [138, 140–142], where the rate at which a process happens depends on the state of adjacent sites.

Experiments on arrays of Rydberg atoms [133, 143–145] have begun to explore these concepts. In the experiment in Ref. [145], a key role is played by dissipative processes which take atoms to inactive states. This leads to a dynamics that drives the system toward the critical point between a low excitation density “absorbing” state and an “active” state. These phenomena are discussed further in Sec. 6.2.4.

There has also been work that combines the use of Rydberg states with coupling to an optical cavity as described in the previous section. Here one of the two laser beams is replaced by coupling to an optical cavity, to create Rydberg polaritons [146]. By performing these experiments in a twisted optical cavity [147], this produced polaritons which behave as if in the presence of an effective magnetic field. Combined with the strong Rydberg interaction this led to the creation of a Laughlin state of polaritons. Since such experiments involve both (dissipative) cavities and excited states of atoms, there is again a clear potential to explore many-body physics of dissipative systems.

## 2.5 Superconducting circuits

Superconducting nano-circuits have played a fundamental role in the field of quantum technologies since the late 1990s, when the first designs based on charge/flux qubits where theoretically studied and experimentally realised. For an extensive review of the early developments in this platform, see Ref. [150]. Superconducting qubits consist of circuits which contain one or a few Josephson junctions in the quantum regime, with each Josephson junction (JJ) formed by two superconducting electrodes separated by a thin oxide layer. In the quantum regime a JJ is described by the Hamiltonian

$$ \hat{H}_{JJ} = \frac{\hat{Q}^2}{2C} - E_J \cos \hat{\varphi}, \qquad (5) $$

where $E_J$ is the Josephson energy associated to the superconducting flow through the junction, $C$ is the capacitance of the junction. The first term in the Hamiltonian is the electrostatic energy while the second is the phase-dependent contribution due to supercurrents. The non-trivial dynamics of this element arises from the commutation relation between the (canonically conjugate) charge $\hat{Q}$ and phase difference $\hat{\varphi}$

$$ \left[ \frac{\hat{Q}}{2e}, e^{\pm i\varphi} \right] = \pm e^{\pm i\varphi}. \qquad (6) $$

A sketch of a JJ and the corresponding lumped circuit scheme as a non-linear LC-resonator is shown in Fig. 8 (panels (a) and (b)). A key feature of the JJ, important for quantum informa-
<figure><img src="image_11.png" /><figcaption>Figure 8: (a) A sketch of a Josephson junction, two superconductors with a difference in the phases of their order parameters $\varphi = \varphi_L - \varphi_R$ separated by a tiny insulating barrier. The junction (labelled by a cross, X) is embedded (b) in a circuit with additional elements. In this simple version only the contribution due to the finite capacitance of the junction is indicated. (c) The Josephson junction is embedded into a stripline resonator, leading to what is known as circuit-QED element. (d) The c-QED element can be arranged in an array with the neighbouring cavities coupled by photon hopping. The gray lines indicate the arrangement of the transmission line realised in Ref. [28, 149]. The red circles indicate where transmon qubits would be placed, and the faint pink lines show the effective kagome lattice. Adapted from [148].</figcaption></figure>

tion processing, is that it is a non-linear, dissipationless element. The anharmonic spectrum of $H_{JJ}$ allows one to single out two levels that can be used to span the Hilbert space of the superconducting qubit. Different choices of the coupling parameters (i.e. the ratio $(2e^2/C)/E_J$) as well different ways to insert it in more complex structures have led to several different qubit designs (charge, flux, phase, ..., qubits). Restricting the dynamics to this computational subspace, the state of the qubit (whatever its origin is) will be labelled by the eigenstates of the Pauli matrix $\hat{\sigma}_z$.

The Josephson element illustrated in Fig. 8 (panels (a) and (b)), placed inside a stripline resonator, is at the heart of circuit-QED. The concept of circuit-QED is analogous to cavity-QED, but replacing optical cavities with microwave resonators, and replacing atoms with superconducting qubits. The microwave resonator, made of superconducting material, is shown in Fig. 8 (panel (c)) and the effective “giant atom” is realised by placing a superconducting nano-circuit inside the stripline resonator. Since the initial proposal in [148], circuit-QED has been a key enabling technology for the growth of solid-state quantum information processing. A comprehensive description of the progresses and status of the art on circuit-QED can be found in the review in [151].

Large arrays of circuit-QED cavities have been theoretically proposed and experimentally realised as simulators of quantum many-body open systems. An sketch of a design for an array of cavities, realised in Ref. [28, 149] is shown in Fig. 8 (panel (d)). Here single cavities are arranged in a regular lattice in such a way to allow photons to hop between neighbouring sites. Note that the wiggly lines in this figure are the cavities, and the “nodes” are the coupling between different cavities.

A circuit-QED array can be modelled with a Hamiltonian that describes qubit-cavity cou-
pling in each cavity, along with photon hopping between cavities. For a single cavity mode with frequency $\omega_c$, the corresponding Hamiltonian is $\hat{H}_c = \omega_c \hat{a}_i^\dagger \hat{a}_i$, where the operator $\hat{a}_i$ ($\hat{a}_i^\dagger$) annihilates (creates) a photon in the i-th cavity. The presence of Josephson elements inside each cavity leads to a strong effective non-linearity between photons. To account for it, it is enough to model the qubit as a few-level system, coupled to a cavity mode. The simplest such model is the Jaynes-Cummings model in which one mode of the cavity interacts with a two-level system:

$$ \hat{H}_i^{(0)} = \frac{\epsilon}{2} \hat{\sigma}_i^z + \omega_c \hat{a}_i^\dagger \hat{a}_i + g(\hat{\sigma}_i^+ \hat{a}_i + \hat{\sigma}_i^- \hat{a}_i^\dagger), \quad (7) $$

$\hat{\sigma}_i^\pm$ are the raising/lowering operators for the two-level system and $\epsilon$ denotes the energy difference between the two levels. In the rotating frame with respect to the uncoupled Hamiltonian the relevant quantity is the detuning $\Delta = \omega_c - \epsilon$. The spectrum of Eq. (7) is anharmonic so that, effectively, the two-level system induces a repulsion between the photons in the cavity.

If the cavities are coupled, a kinetic hopping term $-t_H \sum_{\langle i,j \rangle} (\hat{a}_i^\dagger \hat{a}_j + \text{H.c.})$ should be added to the Hamiltonian, where $t_H$ is the hopping amplitude, and the sum over $\langle i, j \rangle$ denotes neighbouring sites. In the presence of hopping the Hamiltonian of the photons is still quadratic, but now describes photon modes delocalised over the whole lattice. In general, the Hamiltonian for an array of cavities is given by

$$ \hat{H} = \sum_{i} \hat{H}_{i}^{(0)} - t_{H} \sum_{\langle i,j \rangle} (\hat{a}_{i}^{\dagger} \hat{a}_{j} + \text{H.c.}). \quad (8) $$

The interplay between the different terms in the Hamiltonian is responsible for a rich dynamics. A strong effective nonlinearity between the photons turns the cavity into a turnstile device, where no more than one photon can be present on a given site at the same time. On the other hand, photon hopping between neighbouring cavities favours delocalisation and competes with photon blockade.

So far we discussed the Hamiltonian governing the unitary dynamics of the circuit-QED array. The open system dynamics arises because of losses, either photon leakage out of each cavity and/or decay of the qubit. Both these two terms are very well accounted for by a Lindblad term. The explicit form will be described in Sec. 3.1.5, Eq. (27). In the absence of external driving, the steady state will be the one in which all the cavities have no photons and the qubits are in the ground state. In order to reach a non-trivial steady-state, cavities should be refilled by means of an external coherent or incoherent drive.

Additional details of the experimental realisation of coupled cavity arrays will be discussed in Sec. 5.3.4 and can be found in the review in [57].

## 2.6 Polaritons and photons

In most of the platforms discussed so far, light has been considered as a source of external pumping or, via photon emission, as a loss process. However—as already mentioned in Sec. 2.3 on cold atoms in optical cavities—light modes in an optical cavity can also play a role as part of the quantum system to be controlled and studied. Such an idea is key to the physics of polariton and photon condensates. This topic is particularly large, with comprehensive reviews covering a variety of aspects of polariton physics [16, 17, 152], and books e.g. [18, 19] covering key topics. As such, we will only provide a brief and partial summary of essential concepts here, to serve as context for places where we later discuss methods and ideas that were first developed for polariton and photon condensates. For a thorough review of the status of polariton physics, we refer the reader to one of these many reviews.

Microcavity polaritons are the hybrid particles that result from strong coupling between photon modes trapped in an optical microcavity and material excitations, generally excitons
<figure><img src="image_12.png" /><figcaption>Figure 9: Microcavity polaritons. (a) Cartoon showing a microcavity structure, where an active layer is confined between two distributed Bragg reflectors. (b) Polariton dispersion, arising from strong coupling between photons and excitons. (c) Illustration of parametric scattering. An external pump resonantly excites polaritons at a particular wavevector and energy, and then scatters to signal and idler modes.</figcaption></figure>

(bound electron-hole pairs). When light-matter coupling is strong, one can no longer consider excitons and photons separately as eigenstates, and instead one has new eigenstates, lower and upper polaritons, which are symmetric and anti-symmetric superpositions of photons and excitons. As a result, these particles have hybrid properties of both matter and light. This means that in contrast to bare photons they can show nonlinear interactions. Such strongly coupled polaritons have now been realised in a wide variety of material systems—such realisations have been reviewed for classes of materials including semiconductor quantum wells [16,153], organic materials [154], two-dimensional materials [155], and recently hybrid Perovskites [156]. In addition to exciton-polaritons—on which we focus in this section—there can be polaritons involving a variety of other material excitations, for an overview, see [157].

As shown in Fig. 9(a), a typical microcavity experiment involves photons confined between mirrors formed of distributed Bragg reflectors (DBRs)—layers of alternating dielectric constant [18]. Such DBRs have high reflectivity, thus leading to trapped photon modes with long lifetimes. Because the relevant photon modes are confined in a microcavity, the photons have an effective quadratic dispersion for small in-plane momentum $k_{\perp}$. That is, one can write the photon dispersion as:

$$ \omega_{k_{\perp}} = c\sqrt{k_{\perp}^2 + (n2\pi/L)^2} \approx \omega_{k_{\perp}=0} + \frac{k_{\perp}^2}{2m_{\text{eff}}} + \mathcal{O}(k^4), $$

where $L$ denotes the effective length of the cavity (see Fig. 9(a)), and $n \in \mathbb{N}$ determines the mode index. When the effective cavity length is chosen to make $\omega_{k_{\perp}=0}$ close to resonant with typical exciton energies (i.e. close to the optical bandgap for semiconductors), then the resulting effective mass photon $m_{\text{eff}}$ is small, typically $10^{-4}$ times the electron mass. This small mass is inherited by the polariton (see Fig. 9(b)).

While DBR mirrors can have high quality factors (i.e. low losses), photons will ultimately leak out of the mirrors. As discussed below, this makes the polariton system an open system. This loss is also important in allowing access to information about the polariton system. One may either image the far field of the emission (i.e. angle distribution, as suggested in Fig. 9(a)), or use lenses to image the real space distribution of polaritons in the cavity.

**Polariton condensation.** Because polaritons are a superposition of a photon and a bosonic excitation, they have (approximately) bosonic properties. In particular, one can have a state in which a single polariton mode is macroscopically occupied, analogous to the equilibrium state of Bose-Einstein condensation. Since polaritons are part photon they are inevitably subject to
losses, and so one requires external pumping to compensate this loss. In most experiments seeking condensation, pumping occurs via a non-resonant laser pump, to a highly excited state, followed by subsequent relaxation to populate excitons. This means that rather than the thermal equilibrium Bose-Einstein condensate, polaritons generally adopt a non-equilibrium steady state set by a balance of drive and dissipation [158–160]. Such experiments have indeed allowed the observation of Bose-Einstein condensation of polaritons [161]. The extent to which the system is out of equilibrium can vary, with some experiments with very long polariton lifetimes showing behaviour far closer to equilibrium [162]. In thermal equilibrium, the critical temperature required for Bose-Einstein condensation depends inversely on the mass, hence the small effective polariton mass leads to relatively high critical temperatures, which can range from tens of Kelvin to room temperature depending on the material.

The desire to understand how the properties of a non-equilibrium condensate differ from those of an equilibrium condensate—particularly questions related to superfluidity and phase coherence [163]—helped motivate the development several of the theoretical tools we discuss below. Although the underlying theory is quantum mechanical, much of the phenomenology of polariton condensates—particularly spatial pattern formation and dynamics—is classical. That is, these systems can often be understood in terms of a classical theory of nonlinear waves. This theory takes the form of a driven-dissipative version of the Gross-Pitaevskii equation. We will discuss in Secs. 4.4.3 and 4.7.1 how such a theory can arise as the the classical saddle point of a path integral, or equivalently, the mean-field approximation of a Lindblad master equation. Given our focus on many body quantum systems, we will not discuss further the varieties of behaviour that can result from this equation. For an extensive discussion see Ref. [17].

In addition to experiments with non-resonant excitation, there has also been significant work on optical parametric amplification and oscillation, with resonant drives [164, 165]. In such experiments, a resonant coherent pump directly excites polaritons, which can then scatter coherently to signal and idler modes (see Fig. 9(c)). While the pump here is resonant, this scattering process introduces a free phase, so the transition to an optical parametric oscillation state can be considered as a symmetry-breaking phase transition.

**Quantum polariton and polaritonic lattices.** As noted above, much of the physics of polaritons is essentially a theory of nonlinear classical waves. There are however some limits in which quantum effects can be seen for polaritons. Experiments on polaritons in highly confined geometries [166,167] observed weak anti-bunching of polariton emission, resulting from polariton-polariton interactions. Such experiments involved confining polaritons by using one planar mirror, and a second mirror on the tip of an optical fibre.

Polaritons can also be confined by creating micropillars [168] by etching the edges of the structure sketched in Fig. 9(a). These pillars can then also be engineered into lattices, providing another platform for the coupled cavity arrays discussed in the previous section. Experiments with such arrays have demonstrated topological non-interacting band structures, and lasing in topologically protected modes [169–171]. For reviews of this topic, see Refs. [172, 173].

In addition to improving confinement of polaritons, there are various routes that can be used to increase the intrinsic polariton-polariton interaction strength. Several ideas have been explored in this direction. We mention three of these here. *Dipolar polaritons* [174, 175]—These involve three-way hybridisation between photons, direct excitons (electrons and holes in the same quantum well), and indirect excitons (electrons and holes confined in two different quantum wells). Spatial separation of electrons and holes leads to a permanent dipole moment, increasing interactions, but reduces the coupling to light. The three-way hybridisation allows one to control the trade-off between these two effects. *Rydberg polaritons* [176, 177]—These involve hybridisation of photons with Rydberg excitons, i.e. highly excited excitonic
<figure><img src="image_13.png" /><figcaption>Figure 10: (a) Four-level lasing model. (b) Cartoon of dye laser/photon condensate experiments, with molecules in an optical cavity. Both the electronic state (ground/excited) and vibrational state of the molecules are relevant in coupling to light. (c) Effective lasing scheme using vibrational levels for a dye laser or photon condensate.</figcaption></figure>

states, analogous to Rydberg states of atoms as discussed above. *Trion polaritons* [178]—These involve polaritons in electrically doped materials, so that the optical transitions involved are transitions from an electron to a trion (two electrons and a hole).

An alternative route toward stronger quantum effects is to consider cases which involve few emitters strongly coupled to light. Such a scenario can be realised for various platforms, including coupling a quantum dot to a cavity in a photonic crystal [179,180]. For a review of similar semiconductor platforms, see Ref. [181]. More recently, single emitter strong coupling has been realised with organic molecules coupled to plasmon modes in metal nanoparticles [182].

**Photon condensation and lasing.** A class of experiments that are closely related to those studying polariton condensates are those on photon condensation. The distinction between the two is whether the coupling between light and matter is strong or weak. Photon condensation occurs in the regime of weak coupling, where one can regard the photon modes as the eigenstates. Nonetheless, photon condensation depends on coupling to matter, albeit through incoherent emission and absorption.

In discussing the case of photon condensation and weak light matter coupling, one comes very close to describing a far wider class of driven dissipative systems: photon lasers. In most typical applications of lasers, a quantum description is not really relevant or necessary. As with polariton condensate, a classical theory of nonlinear waves (with nonlinearity due to the gain medium) can often be sufficient [183]. However, one always derive an underlying quantum theory of lasing. One such theory is a model of four-level emitters coupled to a single cavity mode (see Fig. 10(a)), where coherent driving occurs on the 0–1 transition, accompanied by incoherent relaxation on the 1–2 and 3–0 transitions, while lasing occurs on the 2–3 transition, due to coupling to a cavity mode. Such a theory—described by a Lindblad master equation—can then be treated in various levels of approximation, ranging from fully classical, to semiclassical (including spontaneous as well as stimulated emission), and fully quantum—see Refs. [35–37, 184] for a detailed discussion. Quantum treatments can be relevant for understanding photon statistics. One can further make a connection between such a model and some of the generalised Dicke models discussed in the context of cold atoms in cavities. These models can show distinct lasing and superradiant phases [185–187].

For experiments on photon condensation [188], one may understand such experiments as an atypical regime of a dye laser [189]. In a dye laser, the gain medium consists of organic molecules, where electronic transitions are dressed by vibrational modes (see Fig. 10(b,c)). The photon condensate results from operating in a parameter regime where most photons emitted by the molecules can be reabsorbed, as the process of repeated absorption and re-
emission is what leads to thermalisation. In contrast, a standard dye laser is designed so that the absorption and emission are spectrally separated, so that laser light emitted by the photons will not be reabsorbed. When absorption and emission are strong enough compared to photon loss rates, the distribution of photons in a cavity can approach a quasi-equilibrium state. In contrast to true thermal equilibrium—in which the chemical potential for photons is stuck at zero—the photon condensate has a chemical potential set by the balance between external driving and loss. Combining this non-zero chemical potential with the effective quadratic dispersion of photons in a cavity leads to the possibility to realise a form of non-equilibrium Bose-Einstein condensation [188].

# 3 Theoretical frameworks

In this Section we set out the theoretical frameworks that will underpin these lecture notes. We first introduce the master equation—an equation of motion for the density matrix, and some of its various reformulations. We then discuss the Schwinger-Keldysh path integral, an alternative framework, and we show how the two are related. We also discuss the quantum trajectory approach, based on stochastic equations of motion for a pure state, such that ensemble-averaged properties recover those predicted by the master equation. We conclude with a short summary of a number of other frameworks which are used, but which are not so relevant to the rest of these lecture notes.

## 3.1 Lindblad master equation

In an open quantum system, a key requirement is to find a way to model the dynamics of the system without having to explicitly represent the dynamics of the environment to which the system is coupled. Indeed, if one includes in the description the state of the environment as well as the system, one generally has an entangled wavefunction, which implies considerable difficulties when proceeding to its resolution. Using a Schmidt decomposition [190], this state can always be written in the form

$$|\Psi\rangle = \sum_{\mu} \sqrt{p_{\mu}} |\psi_{S,\mu}\rangle |\phi_{E,\mu}\rangle, \qquad (9)$$

where $p_μ > 0$ is the real-valued and positive probability of a given state, $|\psi_{S,μ}\rangle$ are a set of orthonormal states of the system S, and $|\phi_{E,μ}\rangle$ are orthonormal states of the environment E (i.e. $\langle\psi_{S,μ}|\psi_{S,v}\rangle = \delta_{μ,v}$ and $\langle\phi_{E,μ}|\phi_{E,v}\rangle = \delta_{μ,v}$). As noted above, except in very special cases, this sum contains multiple terms, and thus the state is entangled.

The open-quantum-system paradigm is to assume we are only interested in expectations of operators on the system. As such, one may see that using the state above, one has:

$$\langle\Psi|\hat{O}_S|\Psi\rangle = \sum_{\mu} p_{\mu} \langle\psi_{S,\mu}|\hat{O}_S|\psi_{S,\mu}\rangle = \text{Tr}[\hat{O}_S \hat{\rho}_S], \qquad \hat{\rho}_S = \sum_{\mu} p_{\mu} |\psi_{S,\mu}\rangle\langle\psi_{S,\mu}|. \qquad (10)$$

This demonstrates that the expectation of system operators is determined fully by the system density matrix $\hat{\rho}_S$ as defined above, that is generically mixed: if one can find it, one can proceed to calculate all such observables. Since we are not interested in other density matrices, we will suppress the suffix S and write $\hat{\rho}$ for the system density matrix in the rest of these lecture notes.

The key problem is then to find how $\hat{\rho}$ evolves in time without needing to first find the wavefunction of the system and environment. In general this is a challenging task, however for a Markovian open quantum system, we can describe the ensemble-average state of the system by the density matrix $\hat{\rho}$, which evolves under a time-local Liouville-von Neumann equation.
For all the situations we consider in these lecture notes, the equation of motion for the density matrix takes Lindblad form [191]:

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

where the coherent dynamics generated by the system Hamiltonian $\hat{H}$ competes with the dissipative processes described by a set of jump operators $\hat{L}_{\mu}$. We give a derivation of this equation in Appendix A. Here we review some of the key properties of the Lindblad evolution and provide few examples of systems of interest, such as lattice models of quantum spins, bosonic or fermionic particles with different types of jump operators. It is often convenient to introduce the shorthand:

$$ \hat{D}[\hat{L}_{\mu}, \hat{\rho}] = \hat{L}_{\mu} \hat{\rho} \hat{L}_{\mu}^{\dagger} - \frac{1}{2} \{ \hat{L}_{\mu}^{\dagger} \hat{L}_{\mu}, \hat{\rho} \}, \quad (12) $$

to describe the dissipator corresponding to given jump operator $\hat{L}_{\mu}$. In most of the examples we will discuss, each of the dissipation terms acts on a single site, and so in general we will label jump operators by site label i as well as a label $\mu$ will be used if multiple operators act per site.

The master equation as written has various notable properties: it preserves the normalization of the density matrix $\text{Tr}[\hat{\rho}] = 1$, because $\text{Tr}[\frac{d}{dt}\hat{\rho}] = 0$. It preserves Hermiticity, i.e. the relation $\text{Tr}[\frac{d}{dt}\hat{\rho} - (\frac{d}{dt}\hat{\rho})^{\dagger}] = 0$ is correct as long as $\hat{\rho} = \hat{\rho}^{\dagger}$. The Lindblad master equation also preserves positivity of the eigenvalues of $\hat{\rho}$, as well as a stronger constraint—complete positivity—which means that not only is $\hat{\rho}$ positive, but also any combined density matrix of $\hat{\rho}$ and another system remains positive. The Lindblad form of the master equation is in fact the most general form of Markovian master equation that can guarantee all three of these properties.

The master equation as written above is linear in the system density matrix, so can be formally considered to take the form $\frac{d}{dt}\hat{\rho} = \hat{L}[\hat{\rho}]$, where $\hat{L}[\cdot]$ is the Lindbladian superoperator (such Lindbladian superoperators are a subset of “Liouvillian” superoperators, the Liouvillian refers to any superoperator describing the density matrix evolution). This superoperator describes a linear operation mapping operators into operators and thus can in principle be represented as a fourth-order tensor, namely an object with four indices. The master equation is thus amenable to a formal solution introducing the exponential of the Lindbladian superoperator:

$$ \hat{\rho}(t) = \exp(\hat{L}t)[\hat{\rho}(0)]. \quad (13) $$

### 3.1.1 Symmetries

Symmetries play an important role in understanding the collective behaviour of many-body systems, particularly in regard to spontaneous symmetry breaking as discussed in Sec. 6. In closed systems, one can simply consider the symmetries of the Hamiltonian. That is, a symmetry corresponds to a unitary transform $\hat{S}$ under which the Hamiltonian is invariant, $\hat{S}^{\dagger}\hat{H}\hat{S} = \hat{H}$, or, equivalently, which commutes with the Hamiltonian $[\hat{H}, \hat{S}] = 0$, and is thus conserved under the closed-system dynamics thanks to the Ehrenfest theorem. In the case of continuous symmetries, one may also identify the generator of the symmetry as a conserved quantity, in accordance with Noether's theorem. When conserved quantities exist, one can divide the Hilbert space into different sectors, corresponding to eigenvalues of the conserved quantity. For an initial state within a given sector, the dynamics will always remain within that sector.

For open systems described by a Lindblad master equation, the properties of symmetries are more subtle [192, 193]. One may distinguish two classes of symmetries:
**Weak symmetries.** In this case, there exists a unitary transform $\hat{S}$ such that

$$ \hat{\mathcal{L}}[\hat{S}\hat{\rho}\hat{S}^{\dagger}] = \hat{S}\hat{\mathcal{L}}[\hat{\rho}]\hat{S}^{\dagger}. $$

That is, the total Lindbladian is invariant under such a symmetry.

**Strong symmetries.** In this case one has that the Hamiltonian and each of the jump operators are individually invariant under the symmetry. Using the notation of Eq. (11) this means that $[\hat{H}, \hat{S}] = 0$ and additionally $[\hat{L}_{\mu}, \hat{S}] = 0$.

Obviously all strong symmetries are also weak symmetries, but not vice versa. The most significant distinction between a weak and strong symmetry comes when considering the question of conserved quantities, which only truly exist for a strong symmetry. Indeed one can use the Lindblad equation (11) to write an equation of motion for the expectation value of any operator $\hat{O}$, $\langle\hat{O}\rangle = \text{Tr}[\hat{O}\hat{\rho}]$ to obtain

$$ \frac{d\langle\hat{O}\rangle}{dt} = i\langle[\hat{H}, \hat{O}]\rangle + \frac{1}{2}\sum_{\mu}\langle\hat{L}_{\mu}^{\dagger}[\hat{O}, \hat{L}_{\mu}]\rangle + \frac{1}{2}\sum_{\mu}\langle[\hat{L}_{\mu}^{\dagger}, \hat{O}]\hat{L}_{\mu}\rangle, \quad (14) $$

from which it follows that in order to be associated to a conserved quantity of Lindblad dynamics an operator needs to commute with both Hamiltonian and each jump operator.

As a consequence of the above, for problems with a weak (but not strong) symmetry, one may find a unique steady state, while problems with a strong symmetry will have multiple steady states, corresponding to different values of the conserved quantity (see also Refs. [194–199] for a more in-depth discussion on the uniqueness of the steady state and Refs. [200–202] for master equations with multiple entangled steady states). Despite this distinction, a weak symmetry does still lead to a block-diagonal structure of the Lindbladian.

To illustrate the block-diagonal structure for weak symmetries, we may consider a simple example of a single bosonic mode, with incoherent pumping and dissipation:

$$ \hat{H} = \omega_c \hat{a}^{\dagger}\hat{a}, \quad \hat{L}_{-} = \sqrt{\kappa_{-}}\hat{a}, \quad \hat{L}_{+} = \sqrt{\kappa_{+}}\hat{a}^{\dagger}. \quad (15) $$

Such a model has a weak symmetry under $\hat{S} = e^{i\theta\hat{a}^{\dagger}\hat{a}}$. One may see that additionally $\hat{S}$ commutes with the Hamiltonian for this example, but $\hat{S}$ does not commute with the jump operators $\hat{L}_{-,+}$. In the absence of the dissipative terms, the symmetry would imply conservation of the number operator, $\hat{n} = \hat{a}^{\dagger}\hat{a}$. However, the pump and decay processes change the number, so the number is not conserved by the master equation, as anticipated for a weak symmetry. There is nevertheless a block-diagonal structure which follows from the weak symmetry. To appreciate this point it is useful to write the density matrix in the number basis, $\hat{n}|n\rangle = n|n\rangle$. Then, one can write $\rho_n^{(k)} = \langle n|\hat{\rho}|n+k\rangle$, and obtain that the equation of motion for $\rho_n^{(k)}$ only involves other elements of the density matrix with the same $k$. The steady state corresponds to the $k=0$ sector, and there is a unique solution within this sector, since particle number is not conserved.

### 3.1.2 Vectorisation, eigenstates and steady states

As noted above, the Lindblad master equation is a linear differential equation for the system density matrix. As such its behaviour can be understood in terms of the eigenvalues and eigenvectors of this linear operator. Such analysis will be of particular importance in Sec. 6 where we study the nature of dissipative phase transitions [203, 204], and in Sec. 7 where we study the approach to the steady state. Here we introduce these eigenvectors and values, by means of the procedure of “vectorisation” of the density matrix.
Vectorisation amounts to writing operators $\hat{\rho}$ as vectors $|\rho\rangle\rangle$, and superoperators $\hat{\mathcal{L}}[\cdot]$ as operators $\hat{\mathcal{L}}$. These objects live in a vector space with a scalar product equivalent to a trace of the product of operators, $\langle\langle a|b\rangle\rangle = \text{Tr}[\hat{a}^\dagger\hat{b}]$. After vectorisation the master equation takes the form

$$ \frac{d}{dt} |\rho\rangle\rangle = \hat{\mathcal{L}} |\rho\rangle\rangle . $$

It is important to note that the matrix $\hat{\mathcal{L}}$ appearing here is non-Hermitian and can be constructed directly from the form of the master equation, or can be done at the level of operators extending the standard second quantisation formalism. For examples, see [205–212].

In this vectorised form, we can now directly consider the eigenspectrum of $\hat{\mathcal{L}}$. As it is a non-Hermitian matrix it has distinct right-eigenvectors $|r_\mu\rangle\rangle$ and left-eigenvectors $\langle\langle l_\mu|$, with corresponding complex eigenvalues $\lambda_\mu$ (for a broad overview on non-Hermitian matrices, see Ref. [27]). These are defined as:

$$ \hat{\mathcal{L}} |r_\mu\rangle\rangle = \lambda_\mu |r_\mu\rangle\rangle, \qquad (16) $$

$$ \langle\langle l_\mu | \hat{\mathcal{L}} = \lambda_\mu \langle\langle l_\mu | , \qquad (17) $$

where $\langle\langle l_\mu | , | r_\mu \rangle\rangle$ are generally biorthogonal. The second equation can also be written in adjoint form as $\hat{\mathcal{L}}^\dagger |l_\mu\rangle\rangle = \lambda_\mu^* |l_\mu\rangle\rangle$. We choose a normalisation such that they are also biorthonormal:

$$ \langle\langle l_\mu | r_\nu \rangle\rangle = \delta_{\mu\nu}. \qquad (18) $$

To motivate the interest of this vectorised notation, we show here that within this framework it is very easy to prove that the Lindblad master equation must have at least one steady state. As noted earlier, the master equation is trace preserving, which means in superoperator form $\text{Tr}[\hat{\mathcal{L}}[\hat{\rho}]] = 0$. In vectorised form this reads $\langle\langle 1 | \hat{\mathcal{L}} \rho \rangle\rangle = 0$. As this is valid $\forall \hat{\rho}$, it follows that $\langle\langle 1 |$ is always a left-eigenvector of $\hat{\mathcal{L}}$ with zero eigenvalue, $\langle\langle 1 | \hat{\mathcal{L}} = 0$, and we identify it with the $\mu = 0$ left-eigenvector: $\langle\langle 1 | \equiv \langle\langle l_0 |$ and $\lambda_0 \equiv 0$. As the eigenvalues $\lambda_\mu$ are common for left and right-eigenvectors, there must be also at least one right-eigenvector with zero eigenvalue. We assume here that this right-eigenvector is unique and identify it with the vectorised form of the steady state density matrix, $\hat{\rho}_{ss} |r_0\rangle\rangle \equiv |\rho_{ss}\rangle\rangle$. There are interesting cases in which the Lindbladian has multiple steady states (see for example Ref. [213]), to all of them should correspond a different left eigenvector of $\hat{\mathcal{L}}$. It follows from the orthonormality condition that, with the normalisation fixed by the definition $\langle\langle 1 | \equiv \langle\langle l_0 |$, $\hat{\rho}_{ss}$ is normalised $\text{Tr}[\hat{\rho}_{ss}] = \langle\langle l_0 | \rho_{ss} \rangle\rangle = 1$.

From the orthonormality condition (18) and from the assumption of unique steady state, it follows that all the right eigenstates different from the steady state $|r_{\mu \neq 0}\rangle\rangle$ are traceless: $\langle\langle 1 | r_{\mu \neq 0} \rangle\rangle = \delta_{0,\mu \neq 0} = 0$. We can now interpret the eigenstates of the Lindbladian different from the steady state as decay modes of deviations from the stationary state. Suppose now that at $t = 0$ the system starts in some state $|\rho(0)\rangle\rangle$ that is not the stationary state. After expanding it onto the basis of right eigenvectors of $\hat{\mathcal{L}}$, using the vectorised form of Eq. (13) one can deduce that at later times, the system density matrix will be given by

$$ |\rho(t)\rangle\rangle = |\rho_{ss}\rangle\rangle + \sum_{\mu \neq 0} c_\mu e^{\lambda_\mu t} |r_\mu\rangle\rangle, \qquad (19) $$

with $c_{\mu \neq 0} = \langle\langle l_\mu | \rho(0) \rangle\rangle$ and $c_0 = \langle\langle 1 | \rho_{ss} \rangle\rangle = 1$ by orthonormality (18). We can thus interpret each Lindbladian eigenmode $\mu \neq 0$ as a possible dynamical decay mode of some initial deviation from the steady state, with a decay rate given by $-\text{Re} \lambda_\mu$. In general, a given decay mode will involve both diagonal elements of the density matrix in the energy-eigenstates basis (i.e. populations) as well as off-diagonal elements (i.e. coherences).
<figure><img src="image_14.png" /><figcaption>Figure 11: Sketch of the Lindbladian spectrum. We see the eigenvalues have negative real-part (describing decay modes, guaranteed by the complete positivity of the Lindblad map) and come in general in pairs of complex conjugate eigenvalues. The stationary state is characterised by a zero eigenvalue (zero real part and zero imaginary part) whose existence is also guaranteed by the mathematical structure of this superoperator [Based on figure from Ref. [215]].</figcaption></figure>

A special role is played by the eigenmode with the least-negative real part (excluding zero), which corresponds to the mode that decays most slowly. Thus, at long times (when all other modes have decayed) this mode typically controls the asymptotic decay to the steady state. As such, the real part of this eigenvalue is often referred to as the Asymptotic Decay Rate (ADR):

$$ \lambda_{\text{ADR}} = \min_{\mu>0} (-\text{Re} \, \lambda_{\mu}). \qquad (20) $$

It is also known as the Lindbladian gap or Liouvillian gap (as illustrated in Fig. 11). The behaviour of this eigenvalue will play a key role in the discussion of how dark states are approached at long times (see Sec. 5.2.2), and to understand the behaviour of symmetry breaking for dissipative phase transitions (see Sec. 6). We note nevertheless that not only the eigenvalues but also the expansion coefficients in Eq. (19) can play a role in controlling the long-time limit of the system, as discussed recently in the context of boundary driven systems [214].

### 3.1.3 Multi-time correlations

As discussed above, the Lindblad master equation can be used to find the steady state of the system, or to find how the system evolves toward its steady state. In addition, we will also want below to consider multi-time correlations. For example, the fluorescence spectrum (i.e. inelastic scattering of light) from a driven system is determined by knowing the two-time correlations of a system. We focus here on two-time correlations, and so the task we want to realise is to write:

$$ \langle \hat{O}_{\mu}(t+\tau) \hat{O}_{\nu}(t) \rangle = \text{Tr}_{S+E} \left[ e^{i\hat{H}_{\tau}} \hat{O}_{\mu} e^{-i\hat{H}_{\tau}} \hat{O}_{\nu} \hat{\rho}_{S+E}(t) \right], \qquad (21) $$

where $\hat{H} = \hat{H}_S + \hat{H}_E + \hat{H}_{S-E}$ as discussed above. The challenge, discussed below, is how to write $\langle \hat{O}_{\mu}(t+\tau) \hat{O}_{\nu}(t) \rangle$ without needing to explicitly reintroducing the environment.

Calculating multi-time correlations from a master equation relies on the “quantum regression theorem” [184, 191, 216]. This starts by considering that the Lindblad master equations means that one knows how to write the time-evolution of a single operator in the form $\langle \hat{O}_{\mu}(t+\tau) \rangle = \sum_{\lambda} f_{\mu\lambda}(\tau) \langle \hat{O}_{\lambda}(t) \rangle$. This condition is always true for the Lindblad master equation, as measurement of a sufficient number of operators $\hat{O}_{\lambda}$ at time $t$ (nb, $\hat{O}_{\mu}$ and $\hat{O}_{\lambda}$ are generically different operators) allows full tomography of the density matrix, and hence allows determination of $\langle \hat{O}_{\mu}(t+\tau) \rangle$. That is, by measuring enough operators one can determine the
density matrix, and then evolve that density matrix with the Lindbladian, and thus find the expectation at time $t + \tau$. We will discuss this further in Sec. 3.4.2 in terms of the adjoint equation for operators. The quantum regression theorem then states one can write two-time correlation functions by using the same coefficients $f_{\mu\lambda}(\tau)$ to give:

$$ \langle \hat{O}_{\mu}(t+\tau) \hat{O}_{\nu}(t) \rangle = \sum_{\lambda} f_{\mu\lambda}(\tau) \langle \hat{O}_{\lambda}(t) \hat{O}_{\nu}(t) \rangle. \qquad (22) $$

An equivalent—but often more practical—statement of the same idea is to write:

$$ \langle \hat{O}_{\mu}(t+\tau) \hat{O}_{\nu}(t) \rangle = \mathrm{Tr} \left[ \hat{O}_{\mu} \exp \left( \hat{\mathcal{L}} \tau \right) \left[ \hat{O}_{\nu} \hat{\rho}(t) \right] \right]. \qquad (23) $$

That is, one can calculate two-time correlations by acting with the first operator, continuing to propagate the resulting object $\hat{O}_{\nu}\hat{\rho}$ under the Lindbladian, and then measuring the second operator. Multi-time correlations follow via the obvious generalisation. Given the discussion in the previous section, one may anticipate that all such correlation functions can be expressed in terms of the eigenvalues and vectors of the Lindbladian; we will discuss this point in detail in Sec. 6.1.3.

### 3.1.4 Thermalisation

As discussed in Sec. 1.3, when a system is coupled to a single environment, one should expect that at long times it will reach equilibrium with that environment. In this section, we discuss the extent to which the Lindblad master equation—considered throughout this review—is consistent with this expectation. We note that in this section we focus on systems that are not driven, and are not coupled to multiple baths; as discussed in Sec. 1.3, driven systems are not expected to reach thermal equilibrium, whereas systems in contact with a single thermal bath are expected to.

In line with our previous discussion, we consider two aspects of this question: (1) Does the Lindblad master equation recover the correct steady state (i.e. the mean-force Gibbs state as given in Eqs. (1) and (2))? (2) Does the Lindblad master equation produce two-time correlation functions consistent with the fluctuation dissipation theorem? We discuss each of these in turn.

**Steady state properties.** For simple few-body systems—for example a single cavity mode—it is straightforward to produce a master equation that recovers the expected equilibrium Gibbs state as its steady state. If one considers the system Hamiltonian $\hat{H} = \omega_c \hat{a}^\dagger \hat{a}$, then the corresponding equilibrium Gibbs state should be a diagonal matrix in the photon number basis $|n\rangle$ with elements $\langle n|\hat{\rho}|n\rangle = e^{-\beta n\omega_c}/\mathcal{Z}$. If one considers a master equation with pumping and dissipation, i.e. with jump operators $\hat{L}_- = \sqrt{\kappa_-} \hat{a}$, and $\hat{L}_+ = \sqrt{\kappa_+} \hat{a}^\dagger$, the resulting steady state is diagonal with $\langle n|\hat{\rho}|n\rangle = (\kappa_- / \kappa_+)^n / (1 - \kappa_- / \kappa_+)$. Thus, if one takes $\kappa_- / \kappa_+ = e^{-\beta \omega_c}$ one recovers the equilibrium Gibbs state. This ratio is consistent with what one would find from a standard weak-coupling derivation of the master equation (see appendix A), which gives $\kappa_- = \kappa(n_B(\omega_c) + 1)$, $\kappa_+ = \kappa n_B(\omega_c)$ where $n_B(\omega_c)$ is the Bose-Einstein distribution, $n_B(\omega_c) = [e^{\beta \omega_c} - 1]^{-1}$. One may note that this recovered the Gibbs state for the system Hamiltonian, and not a mean-force Gibbs state: this is consistent with taking the limit of weak system-environment coupling.

While the simple example above does reproduce the equilibrium Gibbs state, this does not generally hold for more complex many-body problems. That is, the forms of master equations we discuss in this review (with local dissipation terms) does not reproduce the Gibbs state
of many body Hamiltonians.<fnref n="1" /> One key point in this is that the equilibrium state of a many body Hamiltonian will typically involve correlations between elements at different positions or different lattice sites. The number of such correlations that must be correctly determined scales with the dimension of the many-body Hilbert space, i.e. is exponential in the number of lattice sites. The steady state of a model with purely local dissipation terms will be restricted in the form of correlations that can be realised, and generally not compatible with giving the correct value for an exponential number of correlations.

In appendix A we discuss the standard Bloch-Redfield approach for deriving the master equation. For a problem on a lattice, this approach will (in general) produce a number of terms exponential in the number of lattice sites. This occurs because the system Hamiltonian generically has a non-degenerate spectrum with this number of eigenstates (i.e. a number of eigenstates equivalent to the dimension of the many-body Hilbert space). The decomposition of the system-bath coupling into eigenoperators of the system Hamiltonian introduces this large number of terms. Following such a derivation is clearly challenging, exactly due to this exponential number of terms. This challenge can be avoided by using an approximate system Hamiltonian in the derivation, neglecting coupling between sites. This will lead to an approximate local dissipation term, with only polynomial numbers of terms. However, this cannot generally recover the true thermal state of the many body Hamiltonian.

There has been extensive discussion of the question of “global” vs “local” master equations [217–222], i.e. comparing the effects of deriving the master equation with the true global system Hamiltonian and nonlocal dissipation terms, vs. using the local system Hamiltonian. There has also been significant work on the related question of how one can construct master equations to reproduce mean-force Gibbs states, even beyond the weak coupling limit, see for example Refs. [222–224].

**Multi-time correlations.** Even when the steady state of the system is correctly described by a master equation, multi-time correlations generally are not. A key point is that the quantum regression theorem in Eq. (23) is generally not consistent with the fluctuation dissipation theorem in Eq. (4). That is, if one uses the quantum regression theorem to calculate the two time correlation functions in (3), their Fourier transforms will not satisfy Eq. (4). This point, as discussed in Refs. [225, 226] holds even for the simple case of a single bosonic mode coupled to a thermal bath discussed above, i.e.

$$ \hat{H} = \omega_c \hat{a}^\dagger \hat{a}, \quad \hat{L}_- = \sqrt{\kappa_-} \hat{a}, \quad \hat{L}_+ = \sqrt{\kappa_+} \hat{a}^\dagger. \qquad (24) $$

If one uses this model and calculates the two-time correlation functions $S_x(\tau)$ and $\chi_x(\tau)$ of the operator $\hat{x} = \hat{a} + \hat{a}^\dagger$ one finds these have the ratio:

$$ \frac{S_x(\omega)}{\mathrm{Im}[\chi_x(\omega)]} = \frac{\kappa_- + \kappa_+}{\kappa_- - \kappa_+} = \coth(\beta \omega_c) \neq \coth(\beta \omega), \qquad (25) $$

where we used the forms of $\kappa_\pm$ quoted above. This expression makes clear exactly why the fluctuation dissipation theorem cannot hold: the correlation is determined by the rates $\kappa_\pm$ which in turn depend on the thermal occupation of the bath at frequency $\omega_c$. This is not sufficient to determine the occupation of the bath at all frequencies $\omega \neq \omega_c$, and so the Markovian Lindblad equation cannot correctly calculate the response of the system at an arbitrary frequency. For many-body systems, the situation can be more subtle, as the many-body Hamiltonian may lead to the bath being sampled at many frequencies, and interactions within the system might also lead to effective thermalisation, see e.g. Refs. [227, 228] for further discussion.

<footnote n="1">¹There are some cases where this still works. For example, some of the problems we consider with only loss terms have a trivial empty state as their steady state. The empty state can be the thermal equilibrium state at zero temperature, and the assumption that only loss terms exist is consistent with assuming a bath at zero temperature.</footnote>
As noted in the introduction, in this review we will focus on Markovian (i.e. weak coupling) models. Moreover, since we focus on driven-dissipative systems and other nonequilibrium states characterised by finite dissipation rate and coupling to multiple baths, we will not generically expect to reach thermal equilibrium. As such, the rest of our discussion will focus on models with local dissipation. For further discussions on open quantum systems thermalizing to a Gibbs state we refer to the literature (see for example Refs. [229,230].)

### 3.1.5 Examples of models

Here we provide some examples of Lindblad master equations for the kinds of models discussed in the rest of these lecture notes.

**Dicke model.** A heavily studied example of a driven-dissipative many-body system is the Dicke model [231]. Such a model can be used to describe many experiments on cold atoms in single-mode optical cavities, as introduced in Sec. 2.3. The Dicke Hamiltonian [108, 109] takes the form:

$$ \hat{H} = \omega_c \hat{a}^{\dagger}\hat{a} + \sum_{i} [\epsilon \hat{\sigma}_{i}^{+}\hat{\sigma}_{i}^{-} + g(\hat{a} + \hat{a}^{\dagger})\hat{\sigma}_{i}^{x}]. \quad (26) $$

Here $\omega_c$ is the cavity frequency (or sometimes cavity-pump detuning) for a cavity mode that couples to all emitters, $\epsilon$ is the energy of each emitter, and $g$ is the emitter-cavity coupling strength. Note that this differs from the Jaynes-Cummings model in Eq. (7) in that this model describes many two-level systems (labelled by $i$), and the matter-light coupling is not in the rotating-wave approximation. Other variants of this model have been studied in some cases. For example, the matter-light coupling can be split into two terms, $g(\hat{a}\hat{\sigma}_{i}^{+} + \text{H.c.}) + g'(\hat{a}\hat{\sigma}_{i}^{-} + \text{H.c.})$ where $g, g'$ denote coupling strengths for “co-rotating” and “counter-rotating” terms. The Dicke model corresponds to $g' = g$. The Tavis-Cummings model (the many-emitter version of the Jaynes-Cummings model) occurs in the rotating-wave approximation and corresponds to setting $g' = 0$, which then leads to a Hamiltonian with a symmetry $\hat{S} = \exp[i\pi(\hat{a}^{\dagger}\hat{a} + \sum_{i}\hat{\sigma}_{i}^{z}/2)]$.

Typical dissipation processes considered in such a model include:

$$ \hat{L}_{c-} = \sqrt{\kappa_-} \hat{a}, \quad \hat{L}_{i-} = \sqrt{\gamma_-} \hat{\sigma}_i^{-}, \quad \hat{L}_{i+} = \sqrt{\gamma_+} \hat{\sigma}_i^{+}, \quad \hat{L}_{i\phi} = \sqrt{\gamma_{\phi}} \hat{\sigma}_i^{z}. \quad (27) $$

The model with only the cavity loss process, $\hat{L}_{c-}$, is most-commonly considered to describe cold atoms in a single-mode cavity [110], however other variants have been studied elsewhere [232, 233].

**XYZ and Ising spin models.** Another class of models of interest will be lattice models described by short-range interactions. A class of models that has been extensively studied in this context is spin models. A fairly generic example of such models is the XYZ model:

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

Here $\hat{\sigma}_i^\alpha$ are Pauli operators, with $\alpha \in \{x, y, z\}$ denoting which operator, and $i$ denoting which site. The couplings $J_\alpha$ are allowed to be different for each vector component of the spins. As in the previous section, the label $\langle i, j \rangle$ on the sum indicates summation over all neighbouring sites. An extensively studied [234] open-system version of this model is formed by combining this Hamiltonian with jump operators:

$$ \hat{L}_{i-} = \sqrt{\gamma} \hat{\sigma}_{i}^{-}, \quad (29) $$
describing loss—taking the spin toward the spin down state—acting independently on each site. In this model one may reasonably consider the case where there is only a single class of jump operator per site, so the sum over jump operators is the sum over sites. If $J_x = J_y$, the model with this Hamiltonian and dissipation will evolve to a trivial state with all spins pointing down. A non-trivial state will exist as long as $J_x \neq J_y$, as then spin-spin interactions include process that do not conserve total $\hat{S}^z = \sum_i \hat{\sigma}_i^z/2$, thus balancing the losses.

Another particularly widely studied spin model is the Ising model, and the Ising model in a transverse field. Various different versions of this model have been studied in different contexts. The Ising model Hamiltonian takes the form:

$$ \hat{H} = J \sum_{\langle i,j \rangle} \hat{\sigma}_i^z \hat{\sigma}_j^z + h \sum_i \hat{\sigma}_i^z. \qquad (30) $$

A driven-dissipative version of this model can be realised by combining this Hamiltonian with jump operators [235]:

$$ \hat{L}_{i-} = \sqrt{\gamma_{-}} \hat{\sigma}_{i}^{-}, \quad \hat{L}_{i+} = \sqrt{\gamma_{+}} \hat{\sigma}_{i}^{+}, \quad \hat{L}_{i\phi} = \sqrt{\gamma_{\phi}} \hat{\sigma}_{i}^{z}. \qquad (31) $$

Here, because the Hamiltonian conserves the total $\hat{S}^z$, one must balance decay with incoherent pumping to achieve a steady state. This model also includes pure dephasing terms, $\hat{L}_{i\phi}$.

The transverse-field Ising model Hamiltonian includes a field transverse to the spin-spin coupling. Using the model as written above, this leads to:

$$ \hat{H} = J \sum_{\langle i,j \rangle} \hat{\sigma}_i^z \hat{\sigma}_j^z + h \sum_i \hat{\sigma}_i^x. \qquad (32) $$

As with the XYZ model, since this Hamiltonian does not conserve $\hat{S}^z$, a non-trivial steady state can exist even for a master equation with only loss terms $\hat{L}_{i-}$ [236–238]. In such a model, while the Hamiltonian has a symmetry under $\pi$-rotation around the x axis, $\hat{S} = \prod_i e^{i\pi\hat{\sigma}_i^x/2}$, the dissipation breaks this symmetry.

An alternate form of the dissipative transverse-field Ising model can be considered, where

$$ \hat{H} = J \sum_{\langle i,j \rangle} \hat{\sigma}_i^x \hat{\sigma}_j^x + h \sum_i \hat{\sigma}_i^z, \qquad (33) $$

along with $\hat{L}_{i-}$ as defined above [239]. This model has a weak symmetry under $\pi$ rotations about the z axis, $\hat{S} = \prod_i e^{i\pi\hat{\sigma}_i^z/2}$.

**Dissipative Hubbard models.** Another archetypal example of short-range lattice models is the Fermi–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}, \qquad (34) $$

written in terms of fermionic annihilation operators $\hat{c}_{i,\sigma}$ describing a fermion on site i with spin $\sigma \in \{\uparrow, \downarrow\}$, and of $\hat{n}_{i,\sigma} = \hat{c}_{i,\sigma}^{\dagger}\hat{c}_{i,\sigma}$, the corresponding number operator. Here $t_H$ defines the amplitude for hopping between neighbouring sites while preserving the spin, and $U$ defines an interaction between fermions on the same site. Because of the Pauli exclusion principle, the only way a site can be multiply occupied is with two fermions of opposite spin. When considered as a model of cold atoms, a relevant source of dissipation is two-body loss defined by the process:

$$ \hat{L}_{j-} = \sqrt{\kappa_{-}^{(2)}} \hat{c}_{j,\uparrow} \hat{c}_{j,\downarrow}. \qquad (35) $$
This describes a pair of fermions on the same site colliding and thus escaping the trap—as discussed in Sec. 2.2 such a process can also be engineered by using a laser to drive photoassociation. As with the interaction term, Pauli exclusion means this contact process can only exist between fermions with opposite spins.

The bosonic version of this model has also been extensively investigated [240], since it can describe both cold-atom experiments [90] and photonic lattices [57], both in equilibrium and in its driven-dissipative variant. The Hamiltonian of the Bose-Hubbard model reads

$$ \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 (36) $$

where each lattice site hosts a single bosonic mode with annihilation operator $\hat{a}_i$. Different kinds of dissipative processes can be considered, from incoherent pumping $\hat{L}_{i+} = \sqrt{\kappa_+} \hat{a}_i^\dagger$, heating due to spontaneous emission $\hat{L}_{i\phi} = \sqrt{\gamma_\phi} \hat{n}_i$ or m-body losses $\hat{L}_{i-(m)} = \sqrt{\kappa_-^{(m)}} (\hat{a}_i)^m$ (including $m=1$ as standard loss).

**Weakly interacting dilute Bose gas.** A model closely related to the Bose-Hubbard model arises when considering the open system of microcavity polaritons [17]. This is the weakly interacting dilute Bose gas (WIDBG). This is defined by a Hamiltonian written in terms of bosonic annihilation operators $\hat{a}_k$ for particles with momentum $k$:

$$ \hat{H} = \sum_{\mathbf{k}} \frac{k^2}{2m} \hat{a}_{\mathbf{k}}^{\dagger} \hat{a}_{\mathbf{k}} + \frac{U}{2V} \sum_{\mathbf{k},\mathbf{k}',\mathbf{q}} \hat{a}_{\mathbf{k}-\mathbf{q}}^{\dagger} \hat{a}_{\mathbf{k}'+\mathbf{q}}^{\dagger} \hat{a}_{\mathbf{k}'} \hat{a}_{\mathbf{k}}, \quad (37) $$

with $m$ being the boson mass, $U$ the (contact) interaction strength, and $V$ the quantisation volume. This has been used to describe both cold atoms [241] as well as microcavity polaritons [17]. This can be thought of as the long-wavelength (i.e. continuum) limit of the Bose-Hubbard model above.

To consider the open system one then adds processes describing gain or loss [163, 242]. Such a model is more transparent to write in terms of operators in real space, $\hat{a}(\mathbf{r}) = \sum_k \hat{a}_k e^{i\mathbf{k}\cdot\mathbf{r}} / L^{d/2}$, which leads to a Lindblad master equation involving an integral rather than a sum of dissipation terms:

$$ \frac{d}{dt} \hat{\rho} = -i[\hat{H}, \hat{\rho}] + \int d^d \mathbf{r} \hat{D}[\sqrt{\kappa_-} \hat{a}(\mathbf{r}), \hat{\rho}] + \hat{D}[\sqrt{\kappa_+} \hat{a}^\dagger(\mathbf{r}), \hat{\rho}] + \hat{D} \left[ \sqrt{\frac{\kappa_-^{(2)}}{2}} \hat{a}(\mathbf{r})^2, \hat{\rho} \right]. \quad (38) $$

Here we have included particle loss at rate $\kappa_-$, gain at rate $\kappa_+$, and two-particle loss at rate $\kappa_-^{(2)}$. The first two of these processes could equivalently be written in momentum space, leading to jump operators $\sqrt{\kappa_-} \hat{a}_k$, $\sqrt{\kappa_+} \hat{a}_k^\dagger$. The third process is more complex in momentum space as it does not lead to a diagonal Lindblad term. This model will be used to illustrate various theoretical approaches in the following sections.

**BCS model of fermion pairing.** In the same way that the WIDBG model can be thought of as the continuum analogue of the Bose-Hubbard lattice model, we next discuss the continuum analogue of the (Fermi-)Hubbard model. The archetypal example of this is the Bardeen-Cooper-Schrieffer (BCS) model that describes superconductivity [243]. Here one considers fermions described by the Hamiltonian

$$ \hat{H} = \sum_{\mathbf{k},\sigma} \frac{k^2}{2m} \hat{c}_{\mathbf{k},\sigma}^{\dagger} \hat{c}_{\mathbf{k},\sigma} - \frac{u}{L} \sum_{\mathbf{k},\mathbf{k}',\mathbf{q}} \hat{c}_{\mathbf{k}-\mathbf{q},\uparrow}^{\dagger} \hat{c}_{\mathbf{k}'+\mathbf{q},\downarrow}^{\dagger} \hat{c}_{\mathbf{k}',\downarrow} \hat{c}_{\mathbf{k},\uparrow}. \quad (39) $$
In contrast to the repulsive interaction in the Hubbard model, the standard BCS model corresponds to considering an attractive contact interaction, denoted here with strength $u$. This attractive interaction can lead to fermion pairing, and thus superconductivity for charged fermions.

The open-system dynamics of such models have been studied considering two-body loss in various works [244,245]. This is relevant to experiments on ultracold fermionic atoms: as with the examples described above, two-body losses can occur, where pairs of particles escape from the trap. In the continuum limit, this process can be described by the master equation:

$$ \frac{d}{dt} \hat{\rho} = -i[\hat{H}, \hat{\rho}] + \int d^d \mathbf{r} \hat{D} \left[ \sqrt{\kappa_-^{(2)}} \hat{c}_\uparrow(\mathbf{r}) \hat{c}_\downarrow(\mathbf{r}), \hat{\rho} \right]. \quad (40) $$

In practice, for numerical convenience, most investigation of dissipative fermion pairing models has used lattices, and so are closer to the dissipative Hubbard model above [244].

**Bose-Fermi models.** A closely related model that has also been studied in a driven dissipative regime is the Bose-Fermi model:

$$ \hat{H} = \sum_{\mathbf{k}} \left[ \omega_{\mathbf{k}} \hat{a}_{\mathbf{k}}^{\dagger} \hat{a}_{\mathbf{k}} + \sum_{i=c,v} \epsilon_{\mathbf{k}}^{i} \hat{c}_{i,\mathbf{k}}^{\dagger} \hat{c}_{i,\mathbf{k}} \right] + \frac{g}{L^{d/2}} \sum_{\mathbf{k},\mathbf{q}} \left( \hat{c}_{c,\mathbf{k}+\mathbf{q}}^{\dagger} \hat{c}_{v,\mathbf{k}} \hat{a}_{\mathbf{q}} + \text{H.c.} \right) \\ + \sum_{\mathbf{k},\mathbf{k}',\mathbf{q}} \sum_{i,i'} \frac{U_{\mathbf{q}}}{L^d} \hat{c}_{i,\mathbf{k}-\mathbf{q},\uparrow}^{\dagger} \hat{c}_{i',\mathbf{k}',\uparrow+\mathbf{q},\downarrow} \hat{c}_{i',\mathbf{k}',\downarrow} \hat{c}_{i,\mathbf{k},\uparrow}. \quad (41) $$

As written here, this describes cavity photons $\hat{a}_{\mathbf{q}}$ interconverting with electron-hole pairs. Here this is written in terms of conduction and valence bands, so creation of an electron-hole pair involves creation of a conduction electrons $\hat{c}_{c,\mathbf{k}+\mathbf{q}}^{\dagger}$ and a valence band hole $\hat{c}_{v,\mathbf{k}}$. This model also includes Coulomb repulsion between electrons. This model corresponds to describing microcavity polaritons [246,247] in terms of their constituent electrons. A similar model has also been studied in the context of pairing of ultracold atoms, when focusing on the two-channel Feshbach regime [248–250].

Driven-dissipative version of this model or closely related models [158, 251, 252] have been studied in the context of nonequilibrium polariton condensation. For such models one typically considers jump operators:

$$ \hat{L}_{a,-,\mathbf{k}} = \sqrt{\kappa} \hat{a}_{\mathbf{k}}, \quad \hat{L}_{c+,i,\mathbf{k}} = \sqrt{\gamma_{+,i,\mathbf{k}}} \hat{c}_{i,\mathbf{k}}^{\dagger}, \quad \hat{L}_{c-,i,\mathbf{k}} = \sqrt{\gamma_{-,i,\mathbf{k}}} \hat{c}_{i,\mathbf{k}}. \quad (42) $$

This describes loss of cavity photons, and exchange of electrons with some reservoir. (The rates of the fermionic process depend on the occupation of the reservoir, hence the assumption of $k$ dependence).

## 3.2 Schwinger-Keldysh path-integral formalism

An alternate—and complementary—way to study an open quantum system is through the Schwinger-Keldysh path integral formulation. A thorough introduction to the Schwinger-Keldysh path integral can be found in the book by Kamenev [54], and a comprehensive discussion of the application of this to the study of driven-dissipative systems in the reviews by Sieberer et al. [31,32] and Thompson and Kamenev [253].

The Schwinger-Keldysh approach is based on a path integral representation of a generating function for correlations. One may start by considering the time evolution of the density matrix
<figure><img src="image_15.png" /><figcaption>Figure 12: Schwinger–Keldysh closed-time-contour, illustrating labelling of two copies of time evolution labelled as forward (f) and backward (b) contours.</figcaption></figure>

in the superoperator form $\frac{d}{dt}\hat{\rho} = \hat{\mathcal{L}}[\hat{\rho}]$, and so writing the formal time evolution:

$$ \hat{\rho}(t_1) = \mathcal{T} \exp \left( \int_{t_0}^{t_1} dt' \hat{\mathcal{L}} \right) \hat{\rho}(t_0), \quad (43) $$

where $\mathcal{T}$ is a time ordering operator. One may then perform the standard path-integral derivation, by dividing the time evolution into a sequence of $N$ steps of length $\delta t$, inserting resolutions of identity at each time step, and then taking the limit $N \to \infty$, $\delta t \to 0$ with $N\delta t = t_1 - t_0$ to produce a continuous time path integral.

Because the density matrix is an operator and $\hat{\mathcal{L}}$ is a superoperator, the construction of the path integral requires inserting two resolutions of identity at each time step (or equivalently, working in a doubled Hilbert space). For purely Hamiltonian evolution, these two copies are associated with the two copies of time evolution seen in writing:

$$ \hat{\rho}(t_1) = \exp(-i\hat{H}(t_1 - t_0))\hat{\rho}(t_0)\exp(i\hat{H}(t_1 - t_0)), \quad (44) $$

thus associating the two copies with forward time evolution (for operators to the left of $\hat{\rho}(t_0)$ and backward time evolution (for operators to the right). This leads to a picture, shown in Fig. 12, where time evolution proceeds on two branches of a single closed-time-contour path. This contour is known as the Schwinger–Keldysh closed-time contour.

In the continuum limit, this yields a path integral expression:

$$ \mathcal{Z} \equiv \frac{\mathrm{Tr}[\hat{\rho}(t_1)]}{\mathrm{Tr}[\hat{\rho}(t_0)]} = \int \mathcal{D}[\Psi_f, \Psi_b] \exp(iS[\Psi_f, \Psi_b]), \quad (45) $$

where $\int \mathcal{D}[\Psi_f, \Psi_b]$ denotes an integral over all time dependent paths of fields

$$ \Psi_f(t) = \begin{pmatrix} \psi_f(t) \\ \psi_f^*(t) \end{pmatrix}, \quad \Psi_b(t) = \begin{pmatrix} \psi_b(t) \\ \psi_b^*(t) \end{pmatrix}, $$

on the forward (f) and backward (b) branches respectively. The labels used to indicate these branches vary between different texts; we use the labels $f, b$ here to make clear the role of the branches as indicating forward and backward propagation.

The action, $S[\Psi_f, \Psi_b]$ can be determined from the Lindbladian by writing:

$$ S[\Psi_f, \Psi_b] = \int_{t_0}^{t_1} dt \left( \psi_f^* i \frac{d}{dt} \psi_f - \psi_b^* i \frac{d}{dt} \psi_b - i \mathcal{L}(\Psi_f, \Psi_b) \right). \quad (46) $$

Here $\mathcal{L}(\Psi_f, \Psi_b)$ is defined by first writing the Lindbladian in normal-ordered form [54, 254] (i.e. creation operators to the left of annihilation operators), and then replacing operators acting to the left of the density matrix by fields on the forward contour, and operators to the right by fields on the backward contour.

Because density matrices are normalised to one, one finds that $\mathcal{Z} = 1$, so the value of $\mathcal{Z}$ does not give useful information about the system. However, one can calculate the expectation
of observables from the path integral, by using generating functionals in terms of classical source fields $J$. Specifically one may write:

$$ \mathcal{Z}[J_f, J_b] = \int \mathcal{D}[\Psi_f, \Psi_b] \exp \left( iS[\Psi_f, \Psi_b] + i \int dt (J_f^T \Psi_f - J_b^T \Psi_b) \right). \quad (47) $$

As above, we use $\Psi_{f,b}$ to represent the vectors $(\psi_{f,b}, \psi_{f,b}^*)$, and we use a similar notation for the source terms $J_{f,b}$. Formally [54] we can define

$$ \mathcal{Z}[J_f, J_b] = \frac{1}{\text{Tr}[\hat{\rho}(t_0)]} \text{Tr} \left[ \exp \left\{ \int_{t_0}^{t_1} \hat{\mathcal{L}}[J_f, J_b] \right\} \hat{\rho}(t_0) \right], \quad (48) $$

where we use a time evolution superoperator $\hat{\mathcal{L}}[J_f, J_b]$ which contains different source terms in the Hamiltonian $\hat{H} \rightarrow \hat{H} + J^T \hat{\Psi}$ on the time evolution to the left and right of the density matrix (i.e. forward and backward contours).

Because $\text{Tr}[\hat{\rho}(t_0)] = 1$, the prefactor in the generating function can be ignored. As such, any correlation function of $\hat{\psi}_{f,b}$ can then be found by taking derivatives with respect to the source fields $J_{f,b}$ and then ultimately setting $J_{f,b}$ to zero. While we introduced the quantities above on a finite time interval $[t_0, t_1]$, one can easily extend by taking $t_{0,1} \rightarrow \mp\infty$, and thus use this approach to find correlations at any time. In an open system, such results would then be expected to be independent of the initial state, $\hat{\rho}(t_0 \rightarrow -\infty)$.

To give a concrete example of the procedure defined above, we may consider the simple example of a bosonic mode with pumping and decay,

$$ \frac{d}{dt} \hat{\rho} = -i[\omega \hat{a}^\dagger \hat{a}, \hat{\rho}] + \hat{\mathcal{D}}[\sqrt{\kappa_-} \hat{a}, \hat{\rho}] + \hat{\mathcal{D}}[\sqrt{\kappa_+} \hat{a}^\dagger, \hat{\rho}]. \quad (49) $$

Applying the rules above leads to a Schwinger-Keldyshpath integral with the action:

$$ S = \int dt \left[ \psi_f^* \left( i \frac{d}{dt} - \omega + \frac{i}{2}(\kappa_+ + \kappa_-) \right) \psi_f - \psi_b^* \left( i \frac{d}{dt} - \omega - \frac{i}{2}(\kappa_+ + \kappa_-) \right) \psi_b - i\kappa_+ \psi_f^* \psi_b - i\kappa_- \psi_b^* \psi_f \right]. \quad (50) $$

This example can be straightforwardly generalised to dissipative many-body models such as the ones discussed in Sec. 3.1.5. The Keldysh action of a dissipative Bose-Hubbard model for example reads

$$ S = \int dt \sum_{\zeta=f,b} s_{\zeta} \left[ \sum_{i} \psi_{i\zeta}^{*} i \frac{d}{dt} \psi_{i\zeta} - H_{\zeta} \right] - i \sum_{i} \left[ L_{ib} L_{if}^{*} - \frac{1}{2} (L_{if}^{*} L_{if} + L_{ib}^{*} L_{ib}) \right], \quad (51) $$

where $s_{f,b} = \pm 1$, $H$ is the Bose-Hubbard Hamiltonian in Eq. (36) written for the boson fields $\psi_{i\zeta}$, $\psi_{i\zeta}^*$ and $L_{i\zeta}$, $L_{i\zeta}^*$ are the jump operators written in terms of the same bosonic fields. In the specific case of dephasing for example we have $L_{i\zeta} = \sqrt{\gamma_\phi} \psi_{i\zeta}^* \psi_{i\zeta}$, so that we get

$$ S = \int dt \sum_{\zeta=f,b} s_{\zeta} \left[ \sum_{i} \psi_{i\zeta}^{*} i \frac{d}{dt} \psi_{i\zeta} + t_H \sum_{\langle i,j \rangle} (\psi_{i\zeta}^{*} \psi_{j\zeta} + hc) - \frac{U}{2} \sum_{i} \psi_{i\zeta}^{*} \psi_{i\zeta} (\psi_{i\zeta}^{*} \psi_{i\zeta} - 1) \right] \\ - i \gamma_{\phi} \sum_{i} \left[ \psi_{ib}^{*} \psi_{ib} \psi_{if}^{*} \psi_{if} - \frac{1}{2} (\psi_{if}^{*} \psi_{if} \psi_{if}^{*} \psi_{if} + \psi_{ib}^{*} \psi_{ib} \psi_{ib}^{*} \psi_{ib}) \right]. \quad (52) $$
Both examples above are bosonic problems and so the fields appearing in the path integral are complex variables. For fermionic problems, one instead finds a path integral over Grassmann variables (see [54, 254] for a discussion of this case). Problems that involve other operators, such as spins, generally require rewriting in terms of bosonic or fermionic operators before constructing a Schwinger–Keldysh action.

It is often convenient to rewrite the Schwinger–Keldysh path integral in terms of symmetric and anti-symmetric combinations of the forward and backward fields; these are known as the classical and quantum fields. As we will see, this rotation simplifies the structure of the quadratic part of the action. Going back to the example of the single bosonic mode with pumping and decay these fields read

$$ \psi_{cl,q} = \frac{1}{\sqrt{2}} (\psi_f \pm \psi_b). \qquad (53) $$

In terms of these rotated fields the non-interacting action in Eq. (50) becomes:

$$ S = \int dt (\psi_{cl}^* \quad \psi_q^*) \begin{pmatrix} 0 & [D_A]^{-1} \\ [D_R]^{-1} & [D^{-1}]_K \end{pmatrix} \begin{pmatrix} \psi_{cl} \\ \psi_q \end{pmatrix}, \qquad (54) $$

where the retarded (R), advance (A), and Keldysh (K) components of the inverse Green's function are given by:

$$ [D_{R,A}]^{-1} = i \frac{d}{dt} - \omega \pm \frac{i}{2}(\kappa_{-} - \kappa_{+}), \quad [D^{-1}]_K = i(\kappa_{-} + \kappa_{+}). \qquad (55) $$

The combination $\kappa_{-} + \kappa_{+}$ appearing in the inverse Keldysh Green's function can be thought as a “noise” strength. The fact that it is the sum of these two rates appearing here describes the fact that both gain and loss lead to noise. We will see further below that this same combination appears as a noise in many different approaches.

As a second example, we may consider the driven-dissipative WIDBG model, often used to describe polariton condensates [163, 242] written in Eq. (38). Using the same rules as noted above, this yields a Keldysh action of the form:

$$ S = \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_-^{(2)}}{4} \right) (\psi_{cl}^{*2} + \psi_q^{*2}) \psi_{cl} \psi_q + \text{c.c.} \right] + i \kappa_-^{(2)} |\psi_{cl}|^2 |\psi_q|^2 \right\}, \qquad (56) $$

where the bare Green's functions (i.e. neglecting corrections due to interactions), written here in real space, take the form:

$$ [D_{R,A}^{(0)}]^{-1} = i \frac{d}{dt} + \frac{\nabla^2}{2m} \pm i(\kappa_{-} - \kappa_{+}), \quad [D^{(0)-1}]_K = i(\kappa_{-} + \kappa_{+}). \qquad (57) $$

This has a similar structure to the first model discussed above, but replacing the single-mode energy $\omega$ with the quadratic dispersion of the WIDBG model.

### 3.3 Quantum trajectories

The third framework for describing open quantum systems that we will discuss is that of quantum trajectories. This framework, also known as *unravelling*, describes the evolution of a pure quantum state due to a combination of Hamiltonian evolution and measurements which monitor the system [255–259]. Physically, the measurements can either be performed by an explicit
experimental apparatus or by the environment (in this latter case we speak of an *ungread* measurement). In quantum mechanics, a measurement is a stochastic process: the state of the quantum system changes depending on the outcome of the measure and the probability of each outcome depends on the state just before the measure is performed. When measurements are repeated multiple times, we speak of distinct *quantum trajectories*, each labelled by its own measurement record. As discussed below, the density matrix that we discussed so far is recovered by averaging over these trajectories. It is interesting to observe that different measurements can produce the same Lindblad master equation and thus the same density matrix $\rho(t)$, even if microscopically, at the level of each single quantum trajectory, the physical properties are extremely different. Indeed, there has been significant recent interest in the observation that different unravellings lead to different amounts of entanglement in the unravelled wavefunction [260–263]—this observation is closely related to the discussion of entanglement transitions which we will discuss in Sec. 8.

The unravelling method can also be seen as a theoretical technique for more efficiently studying a given Lindblad master equations: especially in many-body systems, the study of the density matrix can be numerically costly, as the size of the problem scales as the square of the total Hilbert space dimension. The cost of simulating pure states can present some advantages as it scales only as the Hilbert space dimension, even if the price to pay is that the evolution becomes stochastic and needs thus some form of repetition and averaging (which can often be parallelised) to decrease the statistical uncertainty. The previous statement implies that in these theoretical studies different forms of unravelling can be chosen, all leading to the study of the same density matrix.

We discuss here two frequently encountered ways of unravelling a Lindblad master equation. In both cases, we will produce a form which has both a stochastic noise term, and a non-unitary deterministic evolution (for example, the non-Hermitian Hamiltonian appearing in the quantum jump formalism). We emphasise that there can also be cases—where jump operators are Hermitian—where the unravelling leads to a unitary (but noisy) evolution [264–267], i.e. where no non-Hermitian Hamiltonian is required.

**Quantum jump trajectories.** “Quantum Jump” trajectories describe the case where the noise acts occasionally but abruptly on the quantum state to measure an operator $\hat{L}_\mu$ and is described by the following stochastic Schrödinger equation [256–259, 268, 269]

$$d|\psi(\mathcal{N}_t, t)\rangle = \left[-i\hat{H} - \frac{1}{2}\sum_{\mu}(\hat{L}_{\mu}^{\dagger}\hat{L}_{\mu} - \langle\hat{L}_{\mu}^{\dagger}\hat{L}_{\mu}\rangle_t)\right]dt|\psi(\mathcal{N}_t, t)\rangle + \sum_{\mu}\left(\frac{\hat{L}_{\mu}}{\sqrt{\langle\hat{L}_{\mu}^{\dagger}\hat{L}_{\mu}\rangle_t}} - 1\right)d\mathcal{N}_t^{\mu}|\psi(\mathcal{N}_t, t)\rangle, \quad (58)$$

where $\langle\circ\rangle_t \equiv \langle\psi(\mathcal{N}_t, t)|\circ|\psi(\mathcal{N}_t, t)\rangle$, and $\mathcal{N}_t = \{\mathcal{N}_t^{\mu}\}$ are a set of statistically independent Poisson processes $d\mathcal{N}_t^{\mu} \in \{0, 1\}$ with average value $d\mathcal{N}_t^{\mu} = dt \langle L_{\mu}^{\dagger}L_{\mu}\rangle_t$. This is written assuming $dt$ is sufficiently small that $d\mathcal{N}_t^{\mu} \ll 1$. As such, for most timesteps $d\mathcal{N}_t^{\mu} = 0$ and so the evolution follows the form in the first line. However, occasionally $d\mathcal{N}_t^{\mu} = 1$ for some $\mu$, in which case the second line is non-zero, and will dominate. In those cases the state will change discontinuously, such that the state at the next time $t + dt$ becomes $\hat{L}_{\mu}|\psi(\mathcal{N}_t, t)\rangle / \sqrt{\langle\hat{L}_{\mu}^{\dagger}\hat{L}_{\mu}\rangle_t}$. The density matrix can be recovered by writing:

$$\hat{\rho}(t) = \sum_{\{d\mathcal{N}_{t'} \le t\}} P(\mathcal{N}_t) |\psi(\mathcal{N}_t, t)\rangle \langle\psi(\mathcal{N}_{t'}, t)|, \quad (59)$$
summing over all measurement outcomes at previous times, with $P(N_t)$ describing the probability of a given measurement record.

This relatively complicated equation can be made more intuitively understandable by discussing explicitly how it is solved in practice, for example using a first-order Monte Carlo wave-function scheme [269]. Let us consider the initial pure state of the dynamics $|\psi(t)\rangle$. First, we need to compute the wavefunction $|\psi(t+dt)\rangle = e^{-i(\hat{H}-\frac{i}{2}\sum_{\mu}\hat{L}_{\mu}^{\dagger}\hat{L}_{\mu})dt}|\psi(t)\rangle$ obtained by applying the deterministic part of the evolution. The next step consists in the stochastic evolution. With probability $p_{\mu} = dt \langle\psi(t+dt)|\hat{L}_{\mu}^{\dagger}\hat{L}_{\mu}|\psi(t+dt)\rangle$ the quantum state is abruptly modified to $\hat{L}_{\mu}|\psi(t+dt)\rangle$ and then normalised. Otherwise, with probability $1-\sum_{\mu}p_{\mu}$, the state is unchanged and normalised. At this point, the process is then restarted applying the same procedure as before. Alternative sampling techniques are also possible, for example computing the average time till the next quantum jump and thus sampling the waiting-time distribution [269, 270]. The stochastic nature of the process creates a number of quantum trajectories, whose incoherent average gives the density matrix.

This form of stochastic unravelling often occurs for jump operators $\hat{L}_{\mu} = \sqrt{\kappa}\hat{a}$ corresponding to photon counting measurements: the measurement record consists of a set of times at which a photon was detected (i.e. times at which $N_t^{\mu} = 1$).

**Quantum state diffusion.** A different unravelling occurs when measurements act continuously but weakly. This can occur for example when light from the sample is interfered with a strong reference beam in a balanced homodyne detection scheme. In this case photon counting will occur every time step, but most of those photons come from the reference, not the system, and thus do not represent a measure the system state. One can consider this case directly by using the same equation as above but replacing $\hat{L}_{\mu} \rightarrow \hat{L}_{\mu}+\beta$ and $\hat{H} \rightarrow \hat{H}+(i/2)\sum_{\mu}(\beta\hat{L}_{\mu}^{\dagger}-\beta^{*}\hat{L}_{\mu})$, with $\beta$ denoting the strength of the reference beam. This substitution leaves the master equation unchanged, but clearly changes the probability of measurement at each time step.

In the limit of an infinitely strong reference beam, the system is described by the quantum state diffusion (QSD) equation [255]:

$$d|\psi(\xi_t, t)\rangle = \left[ -i\hat{H} - \frac{1}{2}\sum_{\mu}\left(\langle\hat{L}_{\mu}^{\dagger}\hat{L}_{\mu} + \langle\hat{L}_{\mu}^{\dagger}\rangle_{\xi_t, t}\langle\hat{L}_{\mu}\rangle_{\xi_t, t} - 2\langle\hat{L}_{\mu}^{\dagger}\rangle_{\xi_t, t}\hat{L}_{\mu}\right) \right] dt |\psi(\xi_t, t)\rangle + \sum_{\mu}\left(\hat{L}_{\mu} - \langle\hat{L}_{\mu}\rangle_{\xi_t, t}\right) d\xi_t^{\mu} |\psi(\xi_t, t)\rangle, \quad (60)$$

where $\langle 0 \rangle_{\xi_t, t} = \langle\psi(\xi_t, t)|0|\psi(\xi_t, t)\rangle$ is the average over the quantum state. The first term in Eq. (60) represents the unitary evolution, while the remaining encodes the noise effects. The $d\xi_t^{\mu}$ are $\hat{I}$to increments of a Wiener process $\xi_t = (\xi_t^1, \xi_t^2, \dots)$, responsible for the stochastic nature of the trajectories, with zero mean $d\xi_t^{\mu} = 0$ and fulfilling the exact property $d\xi_t^{\mu} d\xi_t^{\nu} = dt\delta^{\mu\nu}$. The second term in Eq. (60) describes a deterministic back-action from the measuring environment. We note the presence of a feedback mechanism, as the noise couples to the fluctuations of the measured operator $\delta\hat{L}_{\mu} = \hat{L}_{\mu} - \langle\hat{L}_{\mu}\rangle_{\xi_t, t}$. Mathematically, the role of this feedback is to preserve the norm of the state (and any of its cumulants) for any realisation of the noise.

As discussed before for the quantum jump unravelling, it is also the case that for the quantum state diffusion protocol one can reconstruct the density matrix $\hat{\rho}(t)$ by averaging over the noise

$$\hat{\rho}(t) = \int \mathcal{D}\xi_t P(\xi_t) |\psi(\xi_t, t)\rangle\langle\psi(\xi_t, t)|, \quad (61)$$

with $P(\xi_t)$ the probability distribution of the noise at all times up to $t$.
**Relating stochastic evolution to the Lindblad equation.** Monitored dynamics, as for example that governed by the stochastic Schrödinger equation in Eq. (58) or the quantum state diffusion in Eq. (60), are inherently non-linear: the evolution of the state depends on the average value of the jump operators $\langle L_\mu \rangle$, $\langle L_\mu^\dagger \rangle$. Indeed, by following the dynamics along a given trajectory, information is acquired that conditions the subsequent evolution. By averaging over the trajectories, however, one should get $\hat{\rho}(t)$ which evolves according to the (linear) Lindblad equation. Let us now see how the Lindblad equation follows from the stochastic evolution; for simplicity, we will consider this derivation in the case of the quantum state diffusion (an analogous procedure applies similarly for the case of quantum jumps).

The starting point is the evolution of the density matrix along a given trajectory, $\hat{\rho}(\xi_t, t) = |\psi(\xi_t, t)\rangle\langle\psi(\xi_t, t)|$, defining the so-called conditional state. The state at time $t+dt$ is given by

$$ \hat{\rho}(\xi_t, t+dt) = \left[ |\psi(\xi_t, t)\rangle + d|\psi(\xi_t, t)\rangle \right] \left[ \langle\psi(\xi_t, t)| + d\langle\psi(\xi_t, t)| \right]. \quad (62) $$

By substituting the expression for $d|\psi(\xi_t, t)\rangle$ given in Eq. (60) and averaging over the stochastic variable $\xi_t$, one obtains the density matrix

$$ \hat{\rho}(t) = \overline{|\psi(\xi_t, t)\rangle\langle\psi(\xi_t, t)|}, \quad (63) $$

where the average over noise is written explicitly in Eq. (61). Some care should be taken in handling the expression in Eq. (62) because the term $d|\psi(\xi_t, t)\rangle d\langle\psi(\xi_t, t)|$ cannot be simply dropped, as it contains contributions of order $dt$. This can be seen by a close inspection of the right hand side of Eq. (60). When considering $d|\psi(\xi_t, t)\rangle d\langle\psi(\xi_t, t)|$, the parts of $d|\psi(\xi_t, t)\rangle$ proportional to $dt$ can be dropped, but the last term, proportional to $d\xi_t$, should be retained since $(d\xi_t)^2 = dt$. Direct substitution of Eq. (60) in Eq. (62) and averaging over the noise leads to the Lindblad equation. Although quantum trajectories have been traditionally introduced as a computational method to solve the Lindblad master equation by averaging over different noise realisations, more recently interest has shifted towards the understanding of the many-body physics encoded in the conditional state, as we will discuss in more detail in Sec. 8.

## 3.4 Additional theoretical frameworks

The three frameworks described above are those on which we will focus in these lecture notes. However, there are a number of other frameworks that are sometimes used in describing many-body open quantum systems; here we provide a brief summary of some of them.

### 3.4.1 Phase space approaches

Phase space methods refer to representations of the density matrix in terms of distributions over continuous variables [271]. These have primarily been used to represent bosonic systems, but extensions to spin problems also exist [272, 273]. These techniques, originally developed for the study of single-atom and/or single-photon properties in the context of quantum optics, are currently being generalized and applied to the many-body context.

The earliest-studied phase space approach is the Wigner function, $W(z)$, defined in the context of a single bosonic mode $\hat{a}^{(\dagger)}$ by:

$$ W(z) = \frac{1}{\pi^2} \int d^2w e^{-wz^* + w^*z} \text{Tr} \left[ \hat{\rho} e^{w\hat{a}^\dagger - w^*\hat{a}} \right], \quad z \in \mathbb{C}, \quad (64) $$

where $d^2w$ indicates integration over the complex plane of $w$. Using such a representation, the Lindblad master equation translates into a partial differential equation for the Wigner function. The Wigner function can be interpreted as a quasi-probability density, since for
general quantum states, W can be negative. A frequently used approximation is the truncated Wigner approximation [274], which consists of neglecting third and higher-order derivatives in the partial differential equation. This then yields an equation that can be mapped to the Fokker-Planck equation for probability densities, and simulated by trajectories that sample z(t).

Other phase space distributions correspond to different mappings from the density matrix to continuous functions. Among these, it is worth mentioning the positive-P distribution [275] (for a discussion on regimes of applicability see Ref. [276]). This works by creating a distribution over two continuous variables $P(z_1, z_2)$, which can be guaranteed to remain positive. Such methods were originally applied to closed system dynamics, for which they tended to undergo numerical instabilities. For open system dynamics, dissipation damps these instabilities, and so the positive-P distribution has recently been applied to the study of driven-dissipative Bose-Hubbard lattices [277].

### 3.4.2 Heisenberg-Langevin and adjoint equations

The Lindblad master equation, as well as quantum trajectory formalisms, are based on the time evolution of the state of the system. This is the open-system analogue of the Schrödinger picture for closed quantum systems, where states evolve and operators corresponding to observables are fixed in time. The opposite picture is the Heisenberg picture, where states are constant and operators evolve. This can be extended to open quantum systems in various ways.

The most straightforward way is the use of the adjoint equation for time-evolving an operator [271, 278]. The concept is based on the notion of the adjoint operator $\bar{\mathcal{L}}$, that satisfies the following identity for the expectation value of a generic observable O:

$$ \langle O(t) \rangle = \mathrm{Tr}[\hat{O} \hat{\rho}(t)] = \mathrm{Tr}[\hat{O} \exp(\hat{\mathcal{L}}t)[\rho(0)]] = \mathrm{Tr}[\exp(\bar{\mathcal{L}}t)][\hat{O}] \hat{\rho}(0)]. \quad (65) $$

Note that in the final expression, the superoperator $\exp(\bar{\mathcal{L}}t)$ acts on $\hat{O}$ and not on $\hat{\rho}(0)$. For a generic Lindblad master equation in the form of Eq. (11), the action of the adjoint on an operator $\hat{O}$ takes the form:

$$ \bar{\mathcal{L}}[\hat{O}] = i[\hat{H}, \hat{O}] + \sum_{\mu} \left( \hat{L}_{\mu}^{\dagger} \hat{O} \hat{L}_{\mu} - \frac{1}{2} \{ \hat{L}_{\mu}^{\dagger} \hat{L}_{\mu}, \rho \} \right). \quad (66) $$

This is closely related to the quantum regression theorem [184, 191, 216], as introduced in Sec. 3.1.3 above. In terms of the adjoint representation, Eq. (23) can be rewritten as:

$$ \langle \hat{O}_{\mu}(t + \tau) \hat{O}_{\nu}(t) \rangle = \mathrm{Tr} \left[ \exp(\bar{\mathcal{L}}\tau) [\hat{O}_{\mu}] \hat{O}_{\nu} \hat{\rho}(t) \right]. \quad (67) $$

It is tempting, but wrong, to regard this as giving the time evolution of the operator $\hat{O}(t)$: solving $\frac{d}{dt}\hat{a} = \bar{\mathcal{L}}[\hat{a}]$ for the damped harmonic oscillator with Hamiltonian $\hat{H} = \hbar\omega(\hat{a}^{\dagger}\hat{a} + 1/2)$ and quantum jump $\hat{L} = \sqrt{\kappa}\hat{a}$ gives $\hat{a}(t) = e^{-i\omega t - \kappa t/2}\hat{a}(0)$. The exponential decay $e^{-\kappa t/2}$ turns $\hat{a}(t)$ into an operator that does not satisfy the canonical commutation relations $[\hat{a}(t), \hat{a}(t)^{\dagger}] \neq 1$. As we just saw, it is perfectly legitimate to use the adjoint equation to determine the time evolution of the expectation value of an operator, $\frac{d}{dt}\langle\hat{O}\rangle = \langle\bar{\mathcal{L}}[\hat{O}]\rangle$. Using this, one can calculate both $\langle\hat{a}^{\dagger}(t)\hat{a}(t)\rangle$ and $\langle\hat{a}(t)\hat{a}^{\dagger}(t)\rangle$ via the adjoint equation and thus get the correct commutator.

In general, in order to solve this issue, one can consider the Heisenberg-Langevin equation (see e.g. Refs. [184, 271]), where the presence of dissipation introduces in the equation noise
terms in addition to those describing directly the dissipation effects. As an example, for a simple model with a pumped and decaying bosonic mode $\hat{a}$, with dissipation terms

$$ \hat{L}_{-} = \sqrt{\kappa_{-}} \hat{a}, \quad \hat{L}_{+} = \sqrt{\kappa_{+}} \hat{a}^{\dagger}, \qquad (68) $$

the Heisenberg–Langevin equations for the bosonic operators take the form:

$$ \frac{d\hat{a}}{dt} = i[\hat{H}, \hat{a}(t)] - \frac{(\kappa_{-} - \kappa_{+})}{2} \hat{a}(t) + \hat{F}_{a}(t). \qquad (69) $$

$\hat{F}_a(t)$ is a quantum noise, that is, a stochastic operator which averages to zero, $\overline{\hat{F}_a(t)} = 0$ and has Gaussian correlations defined by:

$$ [\hat{F}_a(t), \hat{F}_a^\dagger(t')] = (\kappa_- - \kappa_+) \delta(t - t'), \quad \overline{\{\hat{F}_a^\dagger(t') \hat{F}_a(t)\}} = (\kappa_- + \kappa_+) \delta(t - t'). \qquad (70) $$

The role of the non-commuting noise here is crucial in preserving equal-time commutation relations. Without the noise, the bosonic operator would acquire an exponential decay $\exp(-(\kappa_- - \kappa_+)t/2)$, which is inconsistent with $[\hat{a}(t), \hat{a}(t)^\dagger] = 1$ remaining true at all times. The commutation relation of the noise term however provides a source that compensates this loss, see [184, 191] for details.

The form of Eq. (69) can be derived by considering an explicit model for the environment, writing coupled Heisenberg equations for the system and the environment, and then formally solving the equations for the environment and substituting into the system Heisenberg equation. In such an approach, the noise term arises from the initial value of environment operators. One may note a close relation between the structure of the Heisenberg–Langevin equation and the Schwinger–Keldysh Green’s functions for this problem as given in Eq. (55).

### 3.4.3 Third quantisation

Third quantisation [205, 279] refers to constructing superoperators that describe the action of applying raising and lowering operators to a density matrix. This provides a specific formalisation of the vectorised form of the Lindblad master equation. Specifically it provides a set of superoperators (which become operators in the vectorised form) from which the Lindblad superoperator can be constructed. As discussed further below, for certain problems it also provides a route to constructing the eigenstates of the Lindblad superoperator.

Third quantisation was first discussed for fermionic problems, by introducing Majorana fermionic operators [205]. It was later extended to bosonic systems [279], for which a simpler representation exists:

$$ \begin{aligned} \hat{\dot{a}}_f |\rho\rangle\rangle &= |\hat{a}\hat{\rho}\rangle\rangle, & \hat{\dot{a}}_b |\rho\rangle\rangle &= |\hat{\rho}\hat{a}\rangle\rangle, \\ \hat{\dot{a}}_f^\dagger |\rho\rangle\rangle &= |\hat{a}^\dagger\hat{\rho}\rangle\rangle, & \hat{\dot{a}}_b^\dagger |\rho\rangle\rangle &= |\hat{\rho}\hat{a}^\dagger\rangle\rangle, \end{aligned} \qquad (71) $$

where $f, b$ refer to forward (i.e. left-multiplication) and backward (right-multiplication) superoperators. Using such operators (and paying careful attention to which operators are closest to the density matrix in a given expression), one can rewrite the Lindbladian in terms of these third-quantised operators. This is particularly useful in the context of quadratic problems, as discussed further in Sec. 4.1.

An alternate form of third quantisation was introduced in Ref. [280], using quantum and classical superoperators, by direct analogy to the Keldysh action discussed above. For a bosonic problem these take the form:

$$ \begin{aligned} \hat{\dot{a}}_{cl} |\rho\rangle\rangle &= \frac{1}{\sqrt{2}} |\{\hat{a}, \hat{\rho}\}\rangle\rangle, & \hat{\dot{a}}_q |\rho\rangle\rangle &= \frac{1}{\sqrt{2}} |[\hat{a}, \hat{\rho}]\rangle\rangle, \\ \hat{\dot{a}}_{cl}^\dagger |\rho\rangle\rangle &= \frac{1}{\sqrt{2}} |\{\hat{a}^\dagger, \hat{\rho}\}\rangle\rangle, & \hat{\dot{a}}_q^\dagger |\rho\rangle\rangle &= \frac{1}{\sqrt{2}} |[\hat{a}^\dagger, \hat{\rho}]\rangle\rangle. \end{aligned} \qquad (72) $$
Using such superoperators, one can rewrite the Lindbladian into a form that is directly analogous to the Keldysh action (with fields replaced by third-quantised operators). One can also use these operators to make connections to phase space methods, by using the eigenvectors of these third-quantised operators to form an orthogonal basis for the space of density matrices [280].

### 3.4.4 Effective non-Hermitian Hamiltonian

One early approach to the study of open quantum systems involved the study of non-Hermitian Hamiltonians, attempting to incorporate decay by adding imaginary terms to the Hamiltonian (see the extensive review in Ref. [27] for examples). This approach differs from that introduced above, based on the Lindblad master equation. Here we comment on how the two approaches can be related in some limiting cases. It is important to stress that studies based solely on non-Hermitian Hamiltonians cannot reproduce the full physics described by Lindblad master equations, and can generically be employed only to gain some qualitative information on the dynamics of the setup. We also note that the simulation of a non-Hermitian Hamiltonian may be computationally cheaper than that of the master equation.

To connect our Lindblad master equation, given in Eq. (11), to a non-Hermitian Hamiltonian, we start by performing the following splitting of the Lindbladian superoperator:

$$ \hat{\mathcal{L}}[\hat{\rho}] = \hat{\mathcal{L}}_{\text{nH}}[\hat{\rho}] + \hat{\mathcal{L}}_{\text{qj}}[\hat{\rho}], \qquad (73) $$

with

$$ \hat{\mathcal{L}}_{\text{nH}}[\hat{\rho}] = -i \left( \hat{H} - \frac{i}{2} \sum_{\mu} \hat{L}_{\mu}^{\dagger} \hat{L}_{\mu} \right) \hat{\rho} + i \hat{\rho} \left( \hat{H} + \frac{i}{2} \sum_{\mu} \hat{L}_{\mu}^{\dagger} \hat{L}_{\mu} \right) = -i (\hat{H}_{\text{nH}} \hat{\rho} - \hat{\rho} \hat{H}_{\text{nH}}^{\dagger}), \quad (74a) $$

$$ \hat{\mathcal{L}}_{\text{qj}}[\hat{\rho}] = \sum_{\mu} \hat{L}_{\mu} \hat{\rho} \hat{L}_{\mu}^{\dagger}. \qquad (74b) $$

This splitting shows in a rather direct way that the Lindblad master equation can be interpreted as a dynamics with the effective non-Hermitian Hamiltonian $\hat{H}_{\text{nH}}$, described by $\hat{\mathcal{L}}_{\text{nH}}[\hat{\rho}]$, and by a quantum jump term $\hat{\mathcal{L}}_{\text{qj}}[\hat{\rho}]$. If one considers an initial pure state $\hat{\rho} = |\psi\rangle\langle\psi|$ evolving under $\hat{\mathcal{L}}_{\text{nH}}$, one sees that this evolution is equivalent to evolving $|\psi\rangle$ under the Hamiltonian $\hat{H}_{\text{nH}}$: this dynamics preserves its purity and thus we can deduce that $\hat{\mathcal{L}}_{\text{qj}}[\hat{\rho}]$ is generically responsible for turning pure states into mixed states. At the same time, we may also note that the non-Hermitian dynamics will generically reduce the norm of a state; that is, time evolution of $|\psi\rangle$ under $\hat{H}_{\text{nH}}$ does not preserve $\langle\psi|\psi\rangle$, and thus does not preserve the trace of $|\psi\rangle\langle\psi|$. The jump term $\hat{\mathcal{L}}_{\text{qj}}[\hat{\rho}]$ compensates for this, and enforces the trace preservation that is a generic property of the Lindblad master equation.

This discussion highlights that neglecting $\hat{\mathcal{L}}_{\text{qj}}[\hat{\rho}]$ compromises the possibility of describing the full time-evolution of the system. However, specific fine-tuned situations exist where some properties of the Lindblad master equation can be found using the non-Hermitian dynamics of $\hat{H}_{\text{nH}}$ alone. Let us consider, for instance, an atom in an excited state $|e\rangle$ that can undergo a spontaneous emission process to the ground state $|g\rangle$, described by the jump operator $\hat{L} = \sqrt{\gamma} |g\rangle\langle e|$. If we are interested in the probability of detecting such atom in the excited state, described by the observable $|e\rangle\langle e|$, we can invoke the adjoint equation for the dynamics introduced in Eq. (66) and perform the analogous splitting $\tilde{\mathcal{L}}[|e\rangle\langle e|] = \tilde{\mathcal{L}}_{\text{nH}}[|e\rangle\langle e|] + \tilde{\mathcal{L}}_{\text{qj}}[|e\rangle\langle e|]$. With a few steps of algebra it is possible to show that:

$$ \tilde{\mathcal{L}}_{\text{qj}}[|e\rangle\langle e|] = 0, \qquad (75) $$
by expanding formally the exponential $e^{\hat{\mathcal{L}}t}$ in its series, we can conclude that the Heisenberg dynamics of $|e\rangle\langle e|$ is entirely ruled by the effective non-Hermitian Hamiltonian $\hat{H}_{nH}$. Although this is a simple one-atom problem, similar situations exist also in the many-body context when discussing the effects of atom loss, which we discuss further below. As we will also discuss in Sec. 5, $\hat{H}_{nH}$ can additionally be used as a route to identifying the existence of “dark states” of open system dynamics.

The non-Hermitian dynamics we highlighted above is often referred to as “no-click dynamics”. If one considers the quantum trajectory interpretation of the master equation, the non-Hermitian dynamics described by $\hat{\mathcal{L}}_{nH}[\hat{\rho}]$ is what one obtains if one performs a postselection on all the possible dynamics to select only those histories where no quantum jump has ever occurred—i.e. where the detector shows no “clicks”. While these no-click configurations are in general exponentially rare, there can be situations in which they become relevant for the typical trajectory, as in the case of continuously monitored systems discussed in Sec. 8, in opposition to the averaged Lindblad dynamics.

# 4 Overview of theoretical methods

In this Section we give an overview of (analytical and numerical) methods used to study many-body open quantum systems. These include matrix-product states and tensor networks [281], the Corner-Space Renormalisation Method [282], Gutzwiller mean-field theories and their cluster extensions [283], Dynamical Mean-Field Theories [284], variational and neural network approaches [285–287] and Schwinger–Keldysh field theories [31, 32]. In addition, for small system sizes, there are software packages as QuTip [288, 289] that can be used to simulate the properties of open system. We emphasise that there are several reviews already existing on this topic, therefore we will limit ourselves to discuss some key aspects that are particularly relevant in the context of these lecture notes.

## 4.1 Exact solutions

We start by discussing classes of problems which can be solved exactly. In addition to being useful reference points against which to compare approximate methods, they can also form the starting point for perturbative expansions [290, 291]. In this regard it is worth noting that different notions of “exact solution” can exist. For some problems it is possible to find closed-form expressions for not only the steady state, but for all eigenvalues and right-eigenvectors of the Lindbladian superoperator in Eq. (16). For other problems it may be only possible to present exact solutions for the steady state.

**Quadratic (non-interacting) problems.** Quadratic problems—i.e. those in which the Hamiltonian is quadratic in bosonic or fermionic operators and the jump operators are linear—can be solved exactly. In general such problems take the form:

$$ \hat{H} = \sum_{i,j} \left[ h_{ij} \hat{a}_i^\dagger \hat{a}_j + \frac{1}{2} (k_{ij} \hat{a}_i^\dagger \hat{a}_j^\dagger + \text{H.c.}) \right], \quad \hat{L}_\mu = p_{\mu i} \hat{a}_i + q_{\mu i} \hat{a}_i^\dagger, \quad (76) $$

with $h, k, p, q$ as general matrices.

The fact that such problems can be exactly solved can be seen in various ways: in the context of the Schwinger–Keldysh path integral approach, such problems yield a quadratic action and thus involve Gaussian integrals. This means it is possible to exactly calculate any correlation functions from the Schwinger–Keldysh path integral by means of Wick’s theorem. For bosonic problems, one can also use the Heisenberg–Langevin approach discussed above.
This gives a linear equation that can be solved, and thus used to calculate any desired operator expectation value at later times.

One can also note that if the Lindbladian is quadratic, one can expect steady states to be of Gaussian form, and so one can use the well-known tools of Gaussian quantum information [292–295]. For a Gaussian state, one can express all correlation functions in terms of the covariance matrix. For our general problem, including anomalous terms, this reduces the problem to finding $\langle \hat{a}_i^\dagger \hat{a}_j \rangle$, $\langle \hat{a}_i^\dagger \hat{a}_j^\dagger \rangle$. All other correlation functions can be written in terms of these. For such Gaussian states, one may also note that phase space methods take a particularly simple form: the phase space distribution is also a Gaussian, with parameters determined by the covariance matrix. Furthermore one can use this exact solvability to classify the phases of quadratic open quantum systems [296].

Such quadratic problems can also be solved by third quantisation [205, 279, 280]. Using such an approach one can write explicit expressions for the left and right eigenvectors. These take the form of a ladder of states defined by the action of a set of third-quantised operators acting on the steady state $|\rho_{ss}\rangle$ for right eigenvectors, and a similar ladder of states starting from trace $\langle\langle 1|$. For the simplest case of a single bosonic mode with pumping and decay (i.e. $\hat{H} = \omega_c \hat{a}^\dagger \hat{a}$, $\hat{L}_- = \sqrt{\kappa_-} \hat{a}$, $\hat{L}_+ = \sqrt{\kappa_+} \hat{a}^\dagger$), and using the classical-quantum third-quantised operators given in Eq. (72), the eigenvectors and eigenvalues of the Lindbladian take the form [280]:

$$ 
\begin{aligned}
|r_{\mu, \nu}\rangle\rangle &= \frac{1}{\sqrt{\mu! \nu!}} (\hat{\hat{a}}_q^\dagger)^\mu (\hat{\hat{a}}_q)^\nu |\rho_{ss}\rangle\rangle, \\
\langle\langle l_{\mu, \nu} | &= \frac{1}{\sqrt{\mu! \nu!}} \langle\langle 1 | (\hat{\hat{a}}_{cl} + \eta \hat{\hat{a}}_q)^\mu (\hat{\hat{a}}_{cl}^\dagger - \eta \hat{\hat{a}}_q^\dagger)^\nu |\rho_{ss}\rangle\rangle, \\
\lambda_{\mu, \nu} &= \omega_c (\mu - \nu) - i \frac{\kappa_- - \kappa_+}{2} (\mu + \nu),
\end{aligned} 
\quad (77) $$

where $\eta = (\kappa_- + \kappa_+)/(\kappa_- - \kappa_+)$ gives a dimensionless noise strength.<fnref n="2" /> For a general problem, the ladder operators that appear in defining the eigenstates are linear combinations of the original operators defining the problem. The linear combination involved can be found from the eigenvectors of an effective non-Hermitian Hamiltonian [280].

All the methods discussed above can be directly related by use of the classical and quantum third-quantised approach, as described in Ref. [280]. The connection between Keldysh action and the third quantised form of the Lindbladian was already mentioned above. The connection to phase space methods can be made by considering the eigenstates of the third-quantised operators (as a generalised form of coherent state), which shows:

$$ \hat{\hat{a}}_{cl} |z_{cl}\rangle\rangle = z |w_{cl}\rangle\rangle, \quad \hat{\hat{a}}_{cl}^\dagger |z_{cl}\rangle\rangle = z^* |z_{cl}\rangle\rangle, \quad W(z) = \langle\langle z_{cl} | \rho \rangle\rangle. \quad (78) $$

Such an approach also makes clear the origin of a feature of the non-interacting system dynamics: the eigenvalues of the Lindbladian depend only on $\kappa_- - \kappa_+$, and not on the noise strength $\eta = (\kappa_- + \kappa_+)/(\kappa_- - \kappa_+)$. This independence is clear in the Schwinger–Keldysh formalism (see Eq. (55)). In the third quantised form, it can be seen that a non-unitary similarity transform on the vectorised problem:

$$ \hat{\hat{L}} \rightarrow \hat{\hat{V}}^{-1} \hat{\hat{L}} \hat{\hat{V}}, \quad \hat{\hat{V}} = \exp(-\eta \hat{\hat{a}}_q^\dagger \hat{\hat{a}}_q), \quad (79) $$

which maps the problem to an equivalent one with $\eta = 0$. This structure makes clear why the Lindbladian eigenvalues do not depend on $\eta$, but why $\eta$ does appear in the form of the eigenvectors.

<footnote n="2">²If one writes $\kappa_- = \kappa_0(n_B + 1)$, $\kappa_+ = \kappa_0 n_B$, corresponding to pumping and decay rates for a state with average population given by the Bose-Einstein occupation $n_B(\omega_c)$, then $\eta = 2n_B(\omega_c) + 1 = \coth(\beta \omega_c / 2)$.</footnote>
**Quadratic problems with dephasing.** Another class of problems for which exact results can be obtained are quadratic systems of fermions and bosons with quadratic Hermitian jump operators, e.g.

$$ \hat{H} = \sum_{i,j} h_{ij} \hat{a}_i^\dagger \hat{a}_j, \quad \hat{L}_i = \sqrt{\gamma_\phi} \hat{a}_i^\dagger \hat{a}_i, \qquad (80) $$

where $h$ is a general symmetric matrix $h_{ij} = h_{ji}$. The jump operators $\hat{L}_i$ here describe local dephasing [297–299]. Because the jump operators are quadratic, these models are described by a Lindbladian superoperator which is not of quadratic form. Still, because of the Hermitian nature of the jump operators, one can write down closed equations of motion for the correlation matrix $C_{i,j}(t) = \text{Tr}[\hat{\rho}(t)\hat{a}_i^\dagger\hat{a}_j^\dagger]$. For example, in the specific case of one-dimensional nearest-neighbour hopping, $h_{ij} = J(\delta_{i,j+1} + \delta_{i+1,j})$, we get [293, 295, 300–302]

$$ \frac{d}{dt} C_{i,j} = -iJ(C_{i-1,j}(t) + C_{i+1,j}(t)) + iJ(C_{i,j-1}(t) + C_{i,j+1}(t)) - \gamma_{\phi}(1 - \delta_{ij})C_{i,j}(t). \quad (81) $$

In addition to $C_{i,j}(t)$, single particle dynamical correlation functions, i.e. Green's functions (see Sec. 3.2), can be obtained in closed form [302]. Furthermore, one can show that in general $k$th-order correlation functions also satisfy closed equations of motion. That is, in contrast to the generic case, the equations of motion for $k$th-order correlation functions do not couple to higher-order correlation functions. The exact solvability of these types of models can be seen also by looking at the properties of the Lindbladian (see Ref. [303] and discussion below) or by unravelling the Hermitian jump operator with a classical stochastic noise and resumming the resulting diagrammatic expansion exactly [304]. Because the jump operators are Hermitian, these models describe heating and thermalisation to infinite temperature as we will discuss in Sec. 7. They can also be used to describe diffusive transport when complemented with boundary driving terms [300, 305].

**Kerr-nonlinear oscillator.** A class of interacting problems for which exact solutions also exist is the Kerr nonlinear oscillator. Various different forms of this problem have been discussed [275, 278, 280, 306–313]. A general form of this model that covers most examples studied is given by the Hamiltonian:

$$ \hat{H} = \omega \hat{a}^{\dagger} \hat{a} + \frac{U}{2} \hat{a}^{\dagger} \hat{a}^{\dagger} \hat{a} \hat{a} + (E^{(1)} \hat{a}^{\dagger} + E^{(2)} \hat{a}^{\dagger} \hat{a}^{\dagger} + \text{H.c.}), \quad (82) $$

along with dissipation terms

$$ \hat{L}_{1-} = \sqrt{\kappa_{-}^{(1)}} \hat{a}, \quad \hat{L}_{2-} = \sqrt{\kappa_{-}^{(2)}} \hat{a} \hat{a}, \quad \hat{L}_{1+} = \sqrt{\kappa_{+}^{(1)}} \hat{a}^{\dagger}. \quad (83) $$

The first example considered is the linear driving of an anharmonic oscillator [307], so consists of the above model with $E^{(2)} = 0, \kappa_{-}^{(2)} = 0$. The key observation is that the complex-P distribution [275] for this model can be exactly found, because the corresponding equation of motion has a specific relation between the drift and diffusion terms [37]. Having found the density matrix, this then allows any properties of the steady state to be extracted.

The model with $\kappa_{+}^{(n)} = 0$ (so that dissipative loss is balanced by coherent pumping) has been studied in a number of recent works [278, 311], which have found ways to re-interpret how this exact solution arises. One approach involves the quantum absorber method [311, 314]. This involves considering losses from the system as being fed into a waveguide. If one can then devises a second system that perfectly absorbs all emission from the first system, this allows an exact solution of the dynamics. That is, if one considers a cascaded quantum system [271], with the original system and the second system connected via a chiral waveguide,
the net effect is that the two systems reach a pure state, since there is no emission into the waveguide beyond the second system. This allows for an exact solution.

More recently, these results have been further re-interpreted as a signature of a hidden time-reversal symmetry [278]. This symmetry involves constructing a doubled version of the system by taking the steady state $\hat{\rho}_{ss} = \sum_{\mu} p_{\mu} |\psi_{\mu}\rangle\langle\psi_{\mu}|$ and using this to construct a state $|\Psi_{T}\rangle = \sum_{\mu} \sqrt{p_{\mu}} |\psi_{\mu}\rangle_{A} |\hat{T}\psi_{\mu}\rangle_{B}$. Here $\hat{T}$ is some generalised time-reversal operation, and A, B refer to the two copies of the system. One then constructs $\hat{\rho}_{AB}(0) = |\Psi_{T}\rangle\langle\Psi_{T}|$ and time evolves this under the doubled Lindbladian, $\hat{\mathcal{L}}_{AB} = \hat{\mathcal{L}} \otimes \hat{\mathbb{1}}$. The generalised time-reversal symmetry occurs if there exists some operator $\hat{T}$ such that the evolution of $\hat{\rho}_{AB}(t)$ under $\hat{\mathcal{L}}_{AB}$ satisfies

$$ \mathrm{Tr}[\hat{X}_A \hat{Y}_B \hat{\rho}_{AB}(t)] = \mathrm{Tr}[\hat{Y}_A \hat{X}_B \hat{\rho}_{AB}(t)], \qquad (84) $$

for any pair of operators $\hat{X}, \hat{Y}$. Reference [278] shows that this condition is equivalent to demanding a duality condition of the Lindbladian under a class of mappings. This identification provides an explicit construction of the required quantum absorber system, connecting to the method above. This approach also shows that the hidden time reversal symmetry is broken by thermal noise, so only survives for $\kappa_{+}^{(n)} = 0$.

A different special case is that without coherent pumping, i.e. where $E^{(n)} = 0$ [280, 309, 315]. In this case, the model has a weak symmetry corresponding to the phase of the photon mode, $\hat{S} = e^{i\theta\hat{a}^{\dagger}\hat{a}}$; this symmetry can be used to solve this model. As noted above, weak symmetries are sufficient to ensure the Lindbladian has a block diagonal structure; the key observation of Ref. [315] is that for some specific models (including the incoherently pumped Kerr model when $\kappa_{-}^{(2)} = 0$), each sector can be mapped to an equivalent non-interacting model. This allows one to find the complete Lindbladian eigenspectrum. The model can also be solved by transforming the third-quantised Lindbladian [280]. This has the effect of mapping the problem to a quadratic (and thus solvable model), but with a non-trivial time dependent transform for all operators. This nonetheless allows direct calculation of the time evolution of any initial state.

**Dissipative all-to-all Ising model.** Another example of exactly solvable models are certain versions of the Ising model, as written in Eqs. (30) and (31) [235,315]. One way to understand the solution of this model is by noting it has many weak symmetries [315] under the operations $\hat{S}_i = e^{i\theta\sigma_i^z}$. As with the example of the Kerr oscillator above, this means the Lindbladian has a block diagonal structure, and further elements of the problem make it possible to solve the model in each of these sectors.

There also exists an exact solution for a version of the transverse-field Ising model [316], but here with all-to-all coupling of the spins:

$$ \hat{H} = J \sum_{i \neq j} \hat{\sigma}_i^z \hat{\sigma}_j^z + \sum_i \sum_\alpha h_i^\alpha \hat{\sigma}_i^\alpha, \qquad \hat{L}_0 = \sqrt{\Gamma} \sum_i \sigma_i^-, \qquad \hat{L}_i = \sqrt{\gamma} \sigma_i^-, \qquad (85) $$

with both collective and individual decay, and with both transverse and longitudinal fields $h_i^\alpha$ on each spin. For this model it is possible in certain cases to construct an operator $\hat{\Psi}$ such that in terms of the non-Hermitian Hamiltonian $\hat{H}_{nH} = \hat{H} - \frac{i}{2} \sum_{\mu} L_{\mu}^{\dagger} L_{\mu}$, one has $\hat{H}_{nH}^{\dagger} \hat{\Psi} = \hat{\Psi} \hat{H}_{nH}$ and also $\hat{L}_{\mu} \hat{\Psi} = \hat{\Psi} \hat{L}_{\mu}^{\dagger}$. This is a sufficient condition to mean that $\hat{\rho}_{ss} = \hat{\Psi}\hat{\Psi}^{\dagger}$ is a steady state solution. The operator $\hat{\Psi}$ can be found when either $\Gamma = 0$ or when $h_i^\alpha$ is uniform (independent of site).

**Integrable Lindbladians.** Finally, we close by mentioning recent developments in solving Lindblad problems with integrability techniques, most notably the Bethe Ansatz. We focus
in particular on the case of bulk dissipation [317, 758] (see for example Ref. [318] for the boundary-driven case). The idea behind these approaches is that the Lindblad superopera- tor, when written in vectorised form (as discussed in Sec. 3.1.2), takes the form of a non- Hermitian Hamiltonian which can in certain cases be solved by the Bethe Ansatz. One ex- ample is again the problem of free fermions with dephasing discussed above. This can be mapped to a non-Hermitian Hubbard model with a purely imaginary interaction and solved by Bethe Ansatz [303, 319, 758]. A similar approach can be explored for models of two-level systems with collective dissipation, whose vectorised Lindbladian maps to the non-Hermitian Gaudin-Richardson model, which is again Bethe-Ansatz solvable [320, 321]. Other examples involve mapping to non-Hermitian XXZ spin chains with different types of boundary condi- tions [322, 323]. In all these works one can obtain the full spectrum of the Lindbladian and sometime even obtain matrix elements of local operators from which the dissipative dynamics can be reconstructed.

## 4.2 Cluster expansions

Cluster expansion methods—as we will define below—can be relevant in cases where cer- tain properties of many-body open systems can be extracted through a detailed knowledge of short-range correlations. This means that in some cases, one can extract useful results from exact numerical simulations of finite-sized systems by direct simulation of the master equation, which can be considered as the open-system analogue of exact diagonalisation.

To access behaviour of the infinite system from finite-sized simulations, one can use var- ious scaling methods. A convenient technique is the linked-cluster expansion that has been successfully applied over several decades in classical and quantum statistical mechanics [324]. Powerful series expansions and resummation techniques [325, 326], developed over the years, allow to further determine long-range behaviour close to critical points from the scaling of a given observable for different cluster sizes, as discussed next.

The key idea behind the linked-cluster expansion is related to the fact that extensive quan- tities of lattice systems can be computed, in the thermodynamic limit, by means of series expansions whose different terms associated to clusters of increasing numbers of sites. The way the various terms appear in the series is dictated by the topology of the lattice and by the connectivity of the generator of the thermodynamics (the Hamiltonian, for the study of equi- librium properties) or of the dynamics (the Lindbladian, for the study of stationary properties); here, we will be interested in the latter situation. We briefly outline the method of numeri- cal linked-cluster expansions (NLCEs), where the contributions to the series are obtained by means of exact diagonalisation techniques on finite-size clusters sites [325, 326].

We start with a presentation of the formalism following Ref. [327], to make evident its natural applicability to the study of driven-dissipative quantum systems governed by a Lindblad dynamics. In our discussion, we focus on the steady-state properties, but the approach is immediately adapted to the real-time dynamics as well. The goal is to compute expectation values of generic extensive observables $\hat{O}$ over the stationary density matrix $\hat{\rho}_{ss}$.

The key to the method relies on the fact that the Lindbladian can be written as a sum of local terms; for the sake of simplicity let us assume that only couplings between at most neighbouring sites matter:

$$ \hat{\mathcal{L}} = \sum_{i,j} \alpha_{ij} \hat{\mathcal{L}}_{ij}, \quad \text{with } \alpha_{ij} = 0, \text{ if } i, j \text{ are not neighbours.} \qquad (86) $$

The observable $\langle \hat{O} \rangle = \text{Tr}[\hat{O} \hat{\rho}_{ss}]$ can be always formally arranged in an expansion in powers of the $\alpha_{ij}$:

$$ \langle \hat{O} \rangle (\{\alpha_{ij}\}) = \sum_{\{n_{ij}\}} o_{\{n_{ij}\}} \prod_{k} \alpha_{ij}^{n_{ij}}, \qquad (87) $$
where $n_j$ runs over all non-negative integers for each $(i, j)$, such that any possible monomial in the $\alpha_{ij}$ is included; the $O_{\{n_{ij}\}}$ are the coefficients of the resulting polynomial. The expansion (87) can be then organised in clusters:

$$ \langle \hat{O} \rangle (\{\alpha_{ij}\}) = \sum_c W_{[O]}(c), \qquad (88) $$

where each $c$ represents a non-empty set of links between sites, which define the given cluster. Specifically, the so-called cluster weight $W_{[O]}(c)$ contains all terms of the expansion (87), which have at least one power of $\alpha_k$ if $k \in c$, and no powers of $\alpha_k$ if $k \notin c$. Vice-versa, all terms in Eq. (87) can be included in one of these clusters. Note that this writing is fully general, and can be used on finite systems as well as on infinite ones.

Using the inclusion-exclusion principle, we can define $\langle O(c) \rangle = \mathrm{Tr}[\hat{O}\hat{\rho}_{ss}(c)]$ as the steady-state expectation value of the observable $\hat{O}$ calculated for the finite cluster $c$ and using the definition given by the sum in Eq. (88) we obtain the recurrence relation:

$$ W_{[O]}(c) = \langle O(c) \rangle - \sum_{s \subset c} W_{[O]}(s), \qquad (89) $$

where the sum runs over all the sub-clusters $s$ contained in $c$, and $\hat{\rho}_{ss}(c)$ is the steady state of the restricted Lindbladian $\hat{\mathcal{L}}(c)$ defined only over the cluster $c$. An important property of Eq. (89) is that, if $c$ is formed out of two disconnected clusters $c_1$ and $c_2$, its weight $W_{[O]}(c)$ is zero. This follows from the fact that $\langle O \rangle$ is an extensive property $\langle O(c) \rangle = \langle O(c_1) \rangle + \langle O(c_2) \rangle$ and $c = c_1 + c_2$. This means that in the expansion (88) one can simply focus on connected clusters, hence the name of *linked* cluster expansion. Obviously, this implies a significant calculation advantage.

The symmetries of the Lindbladian $\hat{\mathcal{L}}$ may drastically simplify the summation (88), since it is typically not needed to compute all the contributions coming from each connected cluster. This can be immediately seen, e.g., for situations where the interaction term $\alpha$ between different pairs of sites is homogeneous throughout the lattice. In such cases, it is possible to identify the topologically-distinct linked clusters, so that a representative $c_n$ for each class can be chosen and counted according to its multiplicity $\ell(c_n)$ per lattice site (the lattice constant of the graph $c_n$). Here the subscript $n$ denotes the number of $k$-spatial indexes that are grouped in the cluster, that is, its size. The property $\langle O \rangle$ per lattice site can be thus written directly in the thermodynamic limit $N \to \infty$ as:

$$ \frac{\langle \hat{O} \rangle}{N} = \sum_{n=1}^{\infty} \left[ \sum_{\{c_n\}} \ell(c_n) W_{[O]}(c_n) \right]. \qquad (90) $$

The outer sum runs over all possible cluster sizes, while the inner one accounts for all topologically distinct clusters $\{c_n\}$ of a given size $n$. Let us emphasise that, if the series expansion (90) is truncated up to order $n = R$, only clusters $c$ at most of size $R$ have to be considered. Indeed each of them should include at least one power of $\alpha_{i,j}$, $\forall (i, j) \in c$. Therefore a cluster of size $R+1$ or larger does not contribute to the expansion, up to order $\alpha^R$.

As a matter of fact, dealing with open many-body systems significantly reduces our ability to compute large orders in the expansion, with respect to the *closed-system* scenario. The size of the space onto which the linear Lindblad superoperator $\hat{\mathcal{L}}$ acts scales as $d_H^{2n}$, where $d_H$ is the dimension of the local Hilbert space and $n$ is the number of sites of a given cluster. In isolated systems, one would need to evaluate the ground state of the cluster Hamiltonian, of size $d_H^n$. Therefore, for the case of spin-1/2 systems ($d_H = 2$), for clusters up to $n = 8$, the study of $\hat{\mathcal{L}}$ requires to work with linear spaces of dimension $2^{2 \times 8} = 65536$. In addition to the calculation
<figure><img src="image_16.png" /><figcaption>Figure 13: Numerical linked-cluster expansion (NLCE) calculation of the averaged magnetic susceptibility of the model in Eqs. (28), (29). Magnetic susceptibility is averaged over the x-y plane, and plotted as a function of J<sub>y</sub> with other parameters set to J<sub>x</sub> = 0.9, J<sub>z</sub> = 1, γ = 1. Solid lines show the Euler resummed data of the NLCE up to the highest-achievable expansion order, R = 8. The dashed lines show the corresponding bare expansions. The shaded region, 0 < J<sub>y</sub> - 0.9 < 0.02 is where the NLCE converges. Adapted From Ref. [327] [Copyright (2018) by the American Physical Society].</figcaption></figure>

of the observable on a given cluster, algorithms exist to compute all topologically distinct clusters, for a given size and lattice geometry. A remarkable advantage of this method is that it enables a direct access to the thermodynamic limit by only counting the cluster contributions of sizes equal or smaller than a certain size, thus using a limited amount of resources). There is no small parameter which controls this expansion, the actual control parameter for the expansion is given by the amount of correlations that are present in the system. Needless to say, optimised algorithms for cluster generation are crucial to boost the power of this method, see [328] and references therein for additional details on the method with applications to Hamiltonian systems.

Numerical linked-cluster expansions have been used in driven-dissipative systems, for example, in [327, 329]. As an example of the method, we show how it performs for two-dimensional dissipative quantum lattice models of interacting spin-1/2 particles, specifically an anisotropic Heisenberg spin-1/2 model with local (single-site) incoherent spin relaxation, as given by Eq. (28) and (29). We set $J_x = 0.9$, $J_z = 1$, and the dissipation $\gamma = 1$. For $J_y = 0.9$ the Hamiltonian conserves the magnetisation along the z direction, and the steady state is the pure state with all the spins pointing down in the z direction. Away from this singular point, for a certain $J_y > J_{y,c} > 0.9$ the system undergoes a second-order phase transition associated to the spontaneous breaking of the $Z_2$. That is, a transition from a paramagnetic phase for $J_y < J_{y,c}$ to a ferromagnetic (FM) phase for $J_y > J_{y,c}$. In the FM phase, a finite magnetisation in the x-y plane develops: $\langle \hat{\sigma}^x \rangle, \langle \hat{\sigma}^y \rangle \neq 0$, which also defines the order parameter of the transition. This phase transition has recently received a lot of attention, see for example [234, 283, 330, 331] and for this choice of parameters the critical point has been estimated to be $J_{y,c} \sim 1.0 - 1.1$. As an example of linked-cluster expansion at work, we show in Fig. 13 the angular-averaged susceptibility to an external field in the x-y plane as a function of the coupling parameter.
## 4.3 Perturbative expansion

If the Lindblad superoperator can be written in the form

$$ \hat{\mathcal{L}} = \hat{\mathcal{L}}_0 + \epsilon \hat{\mathcal{L}}_1, \qquad (91) $$

where $\hat{\mathcal{L}}_0$ can be exactly solved, and $\epsilon \ll 1$, then the second term on the right hand side of Eq. (91) can be treated as a small perturbation. Perturbation theory for open quantum systems can be formulated in different frameworks. Here we discuss the Lindbladian superoperator approach, while for an approach based on Keldysh field theory we refer to Sec. 4.7 or to the literature [253].

The procedure to make a perturbative expansion of the Lindblad superoperator is identical in spirit to that for the eigenstates and eigenvalues of a Hamiltonian system. The key difference is that the Lindblad superoperator is non-Hermitian, thus requiring non-Hermitian perturbation theory.

A detailed description on how to proceed with perturbation theory in Markovian quantum systems is provided in Refs. [290, 291, 332, 333]. The perturbative approach using the vectorisation described in Section 3.1.2 has been used in particular in Refs. [290, 291]. Here we briefly sketch the idea without dwelling on all the subtleties involved in the procedure; we follow the discussion in [291]. We focus on the steady-state of the full problem, which is given by $|\rho_{ss}\rangle\rangle = \sum_i |\rho_{ss}^{(i)}\rangle\rangle$ (using vectorised notation), and we assume the non-degenerate case for which the unperturbed Lindbladian has a unique steady-state solution. Using the fact that the steady state corresponds to a solution with eigenvalue 0 (at all orders), the order-by-order contributes to the steady state obeys

$$ |\rho_{ss}^{(i)}\rangle\rangle = -\hat{\mathcal{L}}_0^{-1} \hat{\mathcal{L}}_1 |\rho_{ss}^{(i-1)}\rangle\rangle. \qquad (92) $$

Here the inverse $\hat{\mathcal{L}}_0^{-1}$ has to be interpreted as the Moore-Penrose pseudoinverse: using right and left eigenvectors $|r_{\mu}^{(0)}\rangle\rangle$, $\langle\langle l_{\mu}^{(0)}|$ of $\hat{\mathcal{L}}_0$ as defined in Eq. (16) we write

$$ \hat{\mathcal{L}}_0^{-1} = \sum_{\mu>0} \frac{1}{\lambda_{\mu}^{(0)}} |r_{\mu}^{(0)}\rangle\rangle \langle\langle l_{\mu}^{(0)}|, \qquad (93) $$

excluding the steady state contribution $\lambda_0^{(0)} = 0$ from the inverse. Similar expressions can be found not only for the steady state, but for all eigenvectors of the Lindbladian, however the general case is more complicated as it involves non-zero eigenvalues.

As noted earlier, the above result assumes that $\hat{\mathcal{L}}_0 \rho_{ss}^{(0)} = 0$ has a unique solution. There are however important cases where this may not be true. For example, weakly driven systems in the presence of approximate conservation laws as studied in Ref. [333]. Here a particular cases is considered in which the zeroth order term is purely Hamiltonian evolution, $\hat{\mathcal{L}}_0[\hat{\rho}] = -i[\hat{H}_0, \hat{\rho}]$, and the Hamiltonian $\hat{H}_0$ has several conserved quantities $\hat{C}_i$ such that $[\hat{C}_i, \hat{H}_0] = 0$. This means the zeroth order evolution has strong symmetries, leading to a degenerate problem. For cases where the perturbation then breaks these symmetries, a unique solution should exist once including the perturbation. Ref. [333] presents an extension of perturbation theory to these problems, by using two key ingredients. Firstly, a generalised Gibbs ensemble (GGE) for the zeroth order solution, i.e. $\rho_{ss}^{(0)} = \exp(-\sum_i \lambda_i \hat{C}_i)/Z$, where $\hat{C}_i$ are the conserved quantities, $\lambda_i$ are generalised chemical potentials. The values $\lambda_i$ are not determined by $\hat{\mathcal{L}}_0$, reflecting the degeneracy. The values of these parameters are fixed either by requiring $\partial_t \langle\hat{C}_i\rangle = \text{Tr}[\hat{C}_i \epsilon \hat{\mathcal{L}}_1 \hat{\rho}_{ss}^{(0)}] = 0$, or in some cases a more complex condition (see [333] for details). Secondly, to overcome the larger degenerate subspace of $\hat{\mathcal{L}}_0$ when inverting this