Keldysh field theory for driven open quantum systems

Cavity quantum electrodynamics (cavity QED), with its focus on strong light-matter interactions, is a growing field of research, which has experienced several groundbreaking advances in the past few years. Historically, these systems were developed as few or single atom experiments, detecting the radiation properties of atoms that are strongly coupled to a quantized light field. The focus has recently been shifted towards loading more and more atoms inside a cavity. Thereby, not only single particle dynamics in strong radiation fields can be probed, but also collective, macroscopic phenomena, which are driven by light-matter interactions. For an excellent review on this topic, covering important experimental and theoretical developments, see [10]. In the experiments, cold atoms are loaded inside an optical or microwave cavity, for which the coherent interaction between the atomic internal states and a single cavity mode dominates over dissipative processes [117]. The atoms absorb and emit cavity photons, thereby changing their internal states. Due to this process, the spatial modulation of the intra-cavity light field induces a coupling of the cavity photons to the atomic internal state as well as their motional degree of freedom [118]. In this way, cavity photons mediate an effective atom–atom interaction, which leads to a back-action of a single atom on the motion of other atoms inside the cavity. This cavity mediated interaction is long-ranged in space and represents the source of collective effects in cold atomic clouds within cavity QED experiments. One hallmark of collective dynamics in cavity QED has been the observation of self-organization of a Bose–Einstein condensate in an optical cavity. This is accompanied by a Dicke phase transition via breaking of a discrete Z₂-symmetry of the underlying model [9], [119–121].

Although the dominant dynamics in these systems is coherent, there are dissipative effects which cannot be discarded. These are the loss of cavity photons due to imperfections in the cavity mirrors, or spontaneous emission processes of atoms, which emit photons into transverse modes. These effects modify the dynamics of the system on the longest time scales. Therefore, they become relevant for macroscopic phenomena, such as phase transitions and collective dynamics. For instance, dissipative effects have been shown theoretically [122, 123] and experimentally [124] to modify the critical exponent of the Dicke transition compared to its zero temperature value. This illustrates that for the analysis of collective phenomena in cavity QED experiments, the dissipative nature of the system has to be taken into account properly [112, 125].

Typically, for a cavity field which has a very narrow spectrum, the atomic internal degrees of freedom can be reduced to two internal states, whose transitions are nearly resonant to the photon frequency. The operators acting on these two internal states can be represented by Pauli matrices, making them equivalent to a spin-1/2 degree of freedom. A very important model in the framework of cavity QED is the Dicke model [126–128], which describes N atoms (i.e. two-level systems) coupled to a single quantized photon mode. This is expressed by the Hamiltonian (here and in the following we set $\hslash =1$)

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{H}_{\text{D}}}={{\omega}_{0}}\,{{a}^{\dagger}}a+\frac{{{\omega}_{z}}}{2}\underset{i=1}{\overset{N}{\sum}}\,\sigma _{i}^{z}+\frac{g}{\sqrt{N}}\underset{i=1}{\overset{N}{\sum}}\,\sigma _{i}^{x}\left({{a}^{\dagger}}+a\right). \end{array}\end{eqnarray} \tag{ 1 }$$

Here, ${{\omega}_{0}}$ is the photon frequency, ${{\sigma}^{z}}=\,|1\rangle \langle 1|\,-\,|0\rangle \langle 0|$ represents the splitting of the two atomic levels with energy difference ${{\omega}_{z}}$. ${{\sigma}^{x}}=\,|0\rangle \langle 1|\,+\,|1\rangle \langle 0|$ describes the coherent excitation and de-excitation of the atomic state proportional to the atom–photon coupling strength g. The Dicke model features a discrete Z₂ Ising symmetry: it is invariant under the transformation $\left({{a}^{\dagger}},a,\sigma _{i}^{x}\right)\mapsto \left(-{{a}^{\dagger}},-a,-\sigma _{i}^{x}\right)$. In the thermodynamic limit, for $N\to \infty$, it features a phase transition, which spontaneously breaks the Ising symmetry. Crossing the transition, the system enters a superradiant phase, characterized by finite expectation values $\langle a\rangle \ne 0,\langle \sigma _{i}^{x}\rangle \ne 0$. This describes condensation of the cavity photons, i.e. the formation of a macroscopically occupied, coherent intra-cavity field, and a 'ferromagnetic' ordering of the atoms in the x-direction.

Although the Dicke model is a standard model for cavity QED experiments in the ultra-strong coupling limit $\sqrt{N}g>{{\omega}_{z}}{{\omega}_{0}}$, it has been realized only very recently in cold atom experiments, where an entire BEC was placed inside an optical cavity [9]. It has been shown that this setup maps to a Dicke model, with a 'collective' spin degree of freedom [120]. Here, the detuning of the pump laser was chosen such that the atoms effectively remain in the internal ground state, but acquire a characteristic recoil momentum when scattering with a cavity photon. This scattering creates a collective, motionally excited state, which replaces the role of an individual, internally excited atom. The experimental realization of a superradiance transition in the Dicke model is usually inhibited, since the required coupling strength by far exceeds the available value of the atomic dipole coupling. However, for the BEC in the cavity, the energy scales of the excited modes are much lower than the optical scale of the atomic modes. In this way, the superradiance transition indeed became experimentally accessible. This was inspired by a theoretical proposal using two balanced Raman channels between different internal atomic states inside an optical cavity, which reduced the effective level splitting of the internal states to much lower energy scales [129].

In addition to the unitary dynamics represented by the Dicke model, the cavity is subject to permanent photon loss due to imperfections in the cavity mirrors. For high finesse cavities, the coupling of the intra cavity photons to the surrounding vacuum radiation field is very weak, and the latter can be eliminated in a Born-Markov approximation [81]. This results in a Markovian quantum master equation for the system's density matrix

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{\partial}_{t}}\rho =-\text{i}\left[{{H}_{\text{D}}},\rho \right]+{{\mathcal{L}}_{d}}\rho, \end{array}\end{eqnarray} \tag{ 2 }$$

where ρ is the density matrix for the intra cavity system and ${{\mathcal{L}}_{d}}$ adds dissipative dynamics to the coherent evolution of the Dicke model. For a vacuum radiation field, it is given by

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{\mathcal{L}}_{d}}\rho =\kappa \left(L\rho {{L}^{\dagger}}-\frac{1}{2}\left\{{{L}^{\dagger}}L,\rho \right\}\right). \end{array}\end{eqnarray} \tag{ 3 }$$

The Lindblad operator L acting on the density matrix describes pure photon loss (L  =  a) with an effective loss rate κ; the latter depends on system specific parameters, but is typically the smallest scale in the master equation (2) [10].

In generic cavity experiments, there are also atomic spontaneous emission processes. The atoms scatter a laser or cavity photon out of the cavity, and this represents a source of decoherence. This process has been considered in [112], and leads to an effective decay rate of the atomic excited state and therefore to a dephasing of the atoms, described by additional Lindblad operators ${{L}_{i}}=\sigma _{i}^{z}$ [130]. However, these losses are at least three orders of magnitude smaller than the cavity decay rate and typically not considered [9, 125]. The basic model for cavity QED with cold atoms is therefore represented by the Dicke model with dissipation, formulated in terms of the master equation (2). Its dynamics is discussed in section 3, including the Dicke superradiance transition.

The Dicke model Hamiltonian takes the form ${{H}_{\text{D}}}={{H}_{c}}+{{H}_{s}}+{{H}_{cs}}$, where H_c,s,cs represent the bosonic cavity, spin, and spin-boson sectors, respectively. There are many directions to go beyond the Dicke model with single collective spin, still keeping the basic feature of coupling spin to boson (cavity photon) degrees of freedom. One direction—relevant to future cold atom experiments—is to consider multimode cavities instead of a single one. In particular, in conjunction with quenched disorder, intriguing analogies to the physics of quantum glasses can be established in this way [113, 131]. Here, the global coupling of all spins to a single mode ${{g}_{i}}\equiv g$ is replaced by random couplings ${{g}_{i,\ell}}$, where the index $\ell$ now refers to a collection of cavity modes. The many-body physics of such an open system is discussed in section 3.

The basic building block of the Dicke model is the spin-boson term of the form ${{H}_{sc}}\sim \left(a+{{a}^{\dagger}}\right){{\sigma}^{x}}$, i.e. a Rabi type non-linearity that preserves the Z₂ symmetry of H_c,s. In the context of circuit quantum electrodynamics, a natural many-body generalization of Hamiltonians with a spin-boson interaction is to consider entire arrays of microcavities (instead of considering many modes within a single cavity). These cavities can be coupled to each other by single photon tunnelling processes between adjacent cavities, giving rise to Hubbard-type hopping terms $\sim Ja_{i}^{\dagger}{{a}_{j}}+\text{h}.\text{c}.$, where i, j now label the spatial index of the cavities. In cirquit QED, strong non-linearities can be generated, e.g. by coupling to adjacent qubits made of Cooper pair boxes [132, 133]. This gives rise to many-body variants of the Rabi model [134], whose phase diagrams have been studied recently [135–137]. Furthermore, for the implementation of lattice Dicke models with large collective spins, the use of hybrid quantum systems consisting of superconducting cavity arrays coupled to solid-state spin ensembles have been proposed [138]. Spontaneous collective coherence in driven–dissipative cavity arrays has been studied in [139].

In many physical situations (away from the ultra-strong coupling limit), a form of the spin-cavity interaction alternative to the Rabi term is more appropriate, with non-linear building block ${{H}_{cs}}\sim a{{\sigma}^{+}}+{{a}^{\dagger}}{{\sigma}^{-}}$. In fact, this form results naturally from the weak coupling rotating wave approximation of a driven spin-cavity problem. For a single cavity mode and spin, the resulting model is the Jaynes–Cummings model, which in contrast to the Dicke model possesses a continuous U(1) phase rotation symmetry under $a\to {{\text{e}}^{\text{i}\theta}}a,{{\sigma}^{-}}\to {{\text{e}}^{\text{i}\theta}}{{\sigma}^{-}}$.

Clearly, when such systems are driven coherently via a Hamiltonian ${{H}_{d}}=\Omega \left(a+{{a}^{\dagger}}\right)$ or suitable multimode generalizations thereof, both the U(1) and even the Z₂ symmetries of the above models are broken explicitly. Coherent drive is usually the simplest way to compensate for unavoidable losses due to cavity leakage, although incoherent pumping schemes are conceivable using multiple qubits [140]. An advantage of such schemes is that the symmetries of the underlying dynamics are preserved or less severely corrupted in this way (see the discussion in section 2.4). In other platforms, such as exciton-polariton systems (see the subsequent section), incoherent pumping is more natural from the outset.

All the systems discussed here represent genuine instances of driven open many-body systems. Besides the coherent drive, they undergo dissipative processes, which have to be taken into account for a proper understanding of their time evolution and stationary states. A generic feature of these processes is their locality. For the effective spin degrees of freedom, the typical processes are qubit decay (local Lindblad operators ${{L}_{i}}=\sigma _{i}^{-}$) and dephasing (${{L}_{i}}=\sigma _{i}^{z}$). For the bosonic component, local single-photon loss is dominant, ${{L}_{i}}={{a}_{i}}$.

Under various circumstances, such as a low population of the excited spin states, the latter degrees of freedom can be integrated out. In this limit, their physical effect is to generate Kerr-type bosonic non-linearities, giving rise to driven open variants of the celebrated Bose–Hubbard model. These models can even be brought into the correlation dominated regime [24, 132, 133, 141]. Oftentimes, these approximations on the spin sector actually apply, and it is both useful and interesting to study these effective low frequency bosonic theories instead of the full many-body spin-boson problems, see [142–146] for recent work in this direction.

Exciton-polaritons are an extremely versatile experimental platform, which is documented by the richness of physical phenomena that have been studied in these systems both in theory and experiment. For a comprehensive account of the subject, we refer to a number of excellent review articles [13, 147, 148]. A Keldysh functional integral approach is discussed in [149, 150], which provides both a microscopic derivation of an exciton-polariton model and a mean field analysis including Gaussian fluctuations. At this point, we content ourselves with a short introduction, with the aim of showing that in a suitable parameter regime, exciton-polaritons very naturally provide a test-bed to study Bose condensation phenomena out of thermal equilibrium. Similar physics can also arise in a variety of other systems, including condensates of photons [151], magnons [152], and potentially excitons [153]. Remarkably, even cold atoms could be brought to condense in a non- equilibrium regime, where continuous loading of atoms balances three-body losses [154], or in atom laser setups [155–157].

A basic experimental setup for exciton-polaritons consists of a planar semiconductor microcavity embedding a quantum well (see figure 1(a)). This setting allows for a strong coupling of cavity light and matter in the quantum well, as originally proposed in [158]. The free dynamics of the elementary excitations of this system—i.e. of cavity photons and Wannier-Mott excitons—is described by the quadratic Hamiltonian [13]

$$\begin{eqnarray}&&{{H}_{0}}={{H}_{C}}+{{H}_{X}}+{{H}_{C-X}},\end{eqnarray} \tag{ 4 }$$

where the parts of the Hamiltonian involving only photons and excitons, respectively, take the same form, which is given by (here the index α labels cavity photons, $\alpha =C$, and excitons, $\alpha =X$, respectively)⁶

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{H}_{\alpha}}={\int}^{}\frac{\text{d}\mathbf{q}}{{{(2\pi )}^{2}}}\underset{\sigma}{\sum}\,{{\omega}_{\alpha}}(q)a_{\alpha,\sigma}^{\dagger}\left(\mathbf{q}\right){{a}_{\alpha,\sigma}}\left(\mathbf{q}\right). \end{array}\end{eqnarray} \tag{ 5 }$$

Zoom In Zoom Out Reset image size

Field operators $a_{\alpha,\sigma}^{\dagger}\left(\mathbf{q}\right)$ and ${{a}_{\alpha,\sigma}}\left(\mathbf{q}\right)$ create or destroy a photon or exciton (note that both are bosonic excitations) with in-plane momentum $\mathbf{q}$ and polarization σ (there are two polarization states of the exciton which are coupled to the cavity mode [13]). For simplicity, we neglect polarization effects leading to an effective spin–orbit coupling [13]. Due to the confinement in the transverse (z) direction, i.e. along the cavity axis, the motion of photons in this direction is quantized as ${{q}_{z,n}}=\pi n/{{l}_{z}}$, where n is a positive integer, and l_z is the length of the cavity. In writing the Hamiltonian (5), we are assuming that only the lowest transverse mode is populated, which leads to a quadratic dispersion as a function of the in-plane momentum $q=\,|\mathbf{q}|\,=\sqrt{q_{x}^{2}+q_{y}^{2}}$:

$$\begin{eqnarray}&&{{\omega}_{C}}(q)=c\sqrt{q_{z,1}^{2}+{{q}^{2}}}=\omega _{C}^{0}+\frac{{{q}^{2}}}{2{{m}_{C}}}+O\left({{q}^{4}}\right).\end{eqnarray} \tag{ 6 }$$

Here, c is the speed of sound, $\omega _{C}^{0}=c{{q}_{z,1}}$, and the effective mass of the photon is given by ${{m}_{C}}={{q}_{z,1}}/c$. Typically, the value of the photon mass is orders of magnitude smaller than the mass of the exciton, so that the dispersion of the latter appears to be flat on the scale of figure 1(b).

Upon absorption of a photon by the semiconductor, an exciton is generated. This process (and the reverse process of the emission of a photon upon radiative decay of an exciton) is described by

$$\begin{eqnarray}&&{{H}_{C-X}}={{\Omega }_{R}}{\int}^{}\frac{\text{d}\mathbf{q}}{{{(2\pi )}^{2}}}\underset{\sigma}{\sum}\,\left(a_{X,\sigma}^{\dagger}\left(\mathbf{q}\right){{a}_{C,\sigma}}\left(\mathbf{q}\right)+\text{H}.\text{c}.\right),\end{eqnarray} \tag{ 7 }$$

where ${{\Omega }_{R}}$ is the rate of the coherent interconversion of photons into excitons and vice versa. The quadratic Hamiltonian (4) can be diagonalized by introducing new modes—the lower and upper exciton-polaritons, ${{\psi}_{\text{LP},\sigma}}\left(\mathbf{q}\right)$ and ${{\psi}_{\text{UP},\sigma}}\left(\mathbf{q}\right)$ respectively, which are linear combinations of photon and exciton modes. The dispersion of lower and upper polaritons is depicted in figure 1(b). In the regime of strong light-matter coupling, which is reached when ${{\Omega }_{R}}$ is larger than both the rate at which photons are lost from the cavity due to mirror imperfections and the non-radiative decay rate of excitons, it is appropriate to think of exciton-polaritons as the elementary excitations of the system.

In experiments, it is often sufficient to consider only lower polaritons in a specific spin state, and to approximate the dispersion as parabolic [13]. Interactions between exciton-polaritons originate from various physical mechanisms, with a dominant contribution stemming from the screened Coulomb interactions between electrons and holes forming the excitons. Again, in the low-energy scattering regime, this leads to an effective contact interaction between lower polaritons. As a result, the low-energy description of lower polaritons takes the form (in the following we drop the subscript indices in ${{\psi}_{\text{LP},\sigma}}$) [13]

$$\begin{eqnarray}&&{{H}_{\text{LP}}}={\int}^{}\text{d}\mathbf{x}\left[{{\psi}^{\dagger}}\left(\mathbf{x}\right)\left(\omega _{\text{LP}}^{0}-\frac{{{\nabla}^{2}}}{2{{m}_{\text{LP}}}}\right)\psi \left(\mathbf{x}\right)+{{u}_{c}}{{\psi}^{\dagger}}{{\left(\mathbf{x}\right)}^{2}}\psi {{\left(\mathbf{x}\right)}^{2}}\right].\end{eqnarray} \tag{ 8 }$$

While this Hamiltonian is quite generic and arises also, e.g. in cold bosonic atoms in the absence of an external potential, the peculiarity of exciton-polaritons is that they are excitations with relatively short lifetime. In turn, this necessitates continuous replenishment of energy in the form of laser driving in order to maintain a steady population. In figure 1(b), we consider the case in which the excitation laser is tuned to energies well above the lower polariton band. The high-energy excitations thus created are deprived of their excess energy via phonon–polariton and stimulated polariton–polariton scattering. Eventually, they accumulate at the bottom of the lower polariton band. As a consequence of multiple scattering processes, the coherence of the incident laser field is quickly lost, and the effective pumping of lower polaritons is incoherent.

A phenomenological model for the dynamics of the lower polariton field, which accounts for both the coherent dynamics generated by the Hamiltonian (8) and the driven–dissipative one described above, was introduced in [159]. It involves a dissipative Gross–Pitaevskii equation for the lower polariton field that is coupled to a rate equation for the reservoir of high-energy excitations. However, for the study of universal long-wavelength behavior in section 4, this degree of microscopic modeling is actually not required: indeed, any (possibly simpler) model that possesses the relevant symmetries (see the discussion at the beginning of section 4) will yield the same universal physics. Such a model can be obtained by describing incoherent pumping and losses of lower polaritons by means of a Markovian master equation⁷:

$$\begin{eqnarray}&&{{\partial}_{t}}\rho =-\text{i}\left[{{H}_{\text{LP}}},\rho \right]+{{\mathcal{L}}_{d}}\rho,\end{eqnarray} \tag{ 9 }$$

where ${{\mathcal{L}}_{d}}\rho$ encodes incoherent single-particle pumping and losses, as well as two-body losses:

$$\begin{eqnarray}&&{{\mathcal{L}}_{d}}\rho ={\int}^{}\text{d}\mathbf{x}\left({{\gamma}_{p}}\mathcal{D}\left[{{\psi}^{\dagger}}\left(\mathbf{x}\right)\right]\rho +{{\gamma}_{l}}\mathcal{D}\left[\psi \left(\mathbf{x}\right)\right]\rho +2{{u}_{d}}\mathcal{D}\left[\psi {{\left(\mathbf{x}\right)}^{2}}\right]\rho \right),\end{eqnarray} \tag{ 10 }$$

where

$$\begin{eqnarray}&&\mathcal{D}[L]\rho =L\rho {{L}^{\dagger}}-\frac{1}{2}\left\{{{L}^{\dagger}}L,\rho \right\}\end{eqnarray} \tag{ 11 }$$

reflects the Lindblad form, and ${{\gamma}_{p}}$, ${{\gamma}_{l}}$, 2u_d are the rates of single-particle pumping, single-particle loss, and two-body loss, respectively. The inclusion of the non-linear loss term ensures saturation of the pumping. An analogous mechanism is implemented in the above-mentioned Gross–Pitaevskii description. More precisely, in the spirit of universality, the above quantum master equation (9) and the abovementioned phenomenological model reduce to precisely the same low frequency model for bosonic degrees of freedom upon taking the semiclassical limit in the Keldysh path integral associated to equation (9) (see section 2.3 for its implementation), and integrating out the upper polariton reservoir in the phenomenological model.

In recent years, experiments with cold atoms in optical lattices have shown remarkable progress in the simulation of many-body model systems both in and out of equilibrium. A particular strength of cold atom experiments is the unprecedented tuneability of model parameters, such as the local interaction strength and the lattice hopping amplitude. This becomes possible by, e.g. manipulation of the lattice laser and external magnetic fields. It comes along with a very weak coupling of the system to the environment, such that the dynamics can often be seen as isolated on relevant time scales for typical measurements of static, equilibrium correlations. However, more and more experiments start to investigate the realm of non-equilibrium phenomena with cold atoms, e.g. by letting systems prepared in a non-equilibrium initial condition relax in time towards a steady state [37, 39], [163–165]. With these experiments, time scales are reached for which the dissipative coupling to the environment becomes visible in experimental observables. Such dissipation may even hinder the system from relaxation towards a well defined steady state.

A relevant example of a dissipative coupling is decoherence of an atomic cloud induced by spontaneous emission of atoms in the lattice [166, 167]. In this way, the many-body system is heated up, and therefore driven away from the low entropy state in which it was prepared initially. A detailed discussion of the microscopic physics and its long time dynamics can be found in [45–47], see also the review article [48]. For bosonic atoms in optical lattices, the coherent dynamics is described by the Bose–Hubbard Hamiltonian [168, 169]

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{H}_{\text{BH}}}=-J\underset{\langle l,m\rangle}{\sum}\,b_{l}^{\dagger}{{b}_{m}}+\frac{U}{2}\underset{l}{\sum}\,{{n}_{l}}\left({{n}_{l}}-1\right), \end{array}\end{eqnarray} \tag{ 12 }$$

which models bosonic atoms in terms of the creation and annihilation operators $\left[{{b}_{l}},b_{m}^{\dagger}\right]={{\delta}_{lm}}$ in the lowest band of a lattice with site indices l, m. The atoms hop between neighboring lattice sites with an amplitude J, and experience an on-site repulsion U. The lattice potential V(x), which leads to the second quantized form of the Bose–Hubbard model, is created by the superposition of counter-propagating laser beams in each spatial dimension. The coherent laser field couples two internal atomic states via stimulated absorption and emission, which leads to the single particle Hamiltonian

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{H}_{\text{atom}}}=\frac{{{{\mathbf{\hat{p}}}}^{2}}}{2m}-\frac{\Delta }{2}{{\sigma}^{z}}-{{\sigma}^{x}}\frac{\Omega \left(\mathbf{\hat{x}}\right)}{2}, \end{array}\end{eqnarray} \tag{ 13 }$$

where $\mathbf{\hat{p}}$ is the atomic momentum operator, $\Delta$ is the detuning of the laser from the atomic transition frequency, and $\Omega \left(\mathbf{\hat{x}}\right)$ is the laser field at the atomic position. For large detuning, the excited state of the atom can be traced out, which leads to the lattice Hamiltonian

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} H_{\text{atom}}^{\text{eff}}=\frac{{{{\mathbf{\hat{p}}}}^{2}}}{2m}+\frac{|\Omega \left(\mathbf{\hat{x}}\right){{|}^{2}}}{4\Delta }. \end{array}\end{eqnarray} \tag{ 14 }$$

This describes a lattice potential for the ground state atoms generated by the spatially modulated AC Stark shift. Adding an atomic interaction potential and expanding both the single particle Hamiltonian (14) and the interaction in terms of Wannier states, the leading order Hamiltonian is the Bose–Hubbard model (12) [169].

In the semi-classical treatment of the atom-laser interaction, spontaneous emission events are neglected, as their probability is typically very small (see below for a more precise statement). They can be taken into account on the basis of optical Bloch equations [166], which leads, after elimination of the excited atomic state, to an additional, driven–dissipative term in the atomic dynamics. It describes the decoherence of the atomic state due to spontaneous emission, i.e. position dependent random light scattering. The leading order contribution to the dynamics in the basis of lowest band Wannier states is captured by the master equation

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{\partial}_{t}}\rho =-\text{i}\left[{{H}_{\text{BH}}},\rho \right]-\gamma \underset{l}{\sum}\,\left[{{n}_{l}},\left[{{n}_{l}},\rho \right]\right], \end{array}\end{eqnarray} \tag{ 15 }$$

where ${{n}_{l}}=b_{l}^{\dagger}{{b}_{l}}$ is the local atomic density, and ρ is the many-body density matrix. For a red detuned laser the rate $\gamma =\Gamma\frac{|\Omega {{|}^{2}}}{4{{\Delta }^{2}}}$ is proportional to the microscopic spontaneous emission rate $\Gamma$ and the laser amplitude $|\Omega |$. Note the suppression of the scale γ by a factor $\Gamma/\Delta \ll 1$ for large detuning, compared to the strength of $H_{\text{atom}}^{\text{eff}}$. The coherent and incoherent contribution of the atom–laser coupling to the dynamics is illustrated in figure 2.

Zoom In Zoom Out Reset image size

The dissipative term in the master equation (15) leads to an energy increase $\langle {{H}_{\text{BH}}}\rangle (t)\sim t$ linear in time, and thus to heating. Furthermore, it introduces decoherence in the number state basis, i.e. it projects the local density matrix on its diagonal in Fock space and leads to a decrease of the coherences in time [45–48]. Starting from a low entropy state at t  =  0, the heating leads to a crossover from coherence dominated dynamics at short and intermediate times [43] to a decoherence dominated dynamics at long times [45–47]. In one dimension, both regimes have been analyzed extensively both numerically (with a focus on decoherence dominated dynamics) and analytically, and display several aspects of non-equilibrium universality, see section 5. Therefore, heating in interacting lattice systems represents a crucial example for universality in out-of-equilibrium dynamics, which can be probed by cold atom experiments.

The master equation describing the decay of photons in a single-mode cavity takes the general form of equation (16), with $H={{\omega}_{0}}{{a}^{\dagger}}a$, where ${{a}^{\dagger}}$ and a are creation and annihilation operators of photons, and ${{\omega}_{0}}$ is the frequency of the cavity mode. Assuming the external electromagnetic field to be in the vacuum state, there is only a single term in the sum over α with L  =  a, describing the decay of the cavity field at a rate $2\kappa$ (the factor of 2 is chosen for convenience). The corresponding Keldysh action is given by [112]

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} S={{{\int}^{}}_{t}}\left\{a_{+}^{\ast}\left(\text{i}{{\partial}_{t}}-{{\omega}_{0}}\right){{a}_{+}}-a_{-}^{\ast}\left(\text{i}{{\partial}_{t}}-{{\omega}_{0}}\right){{a}_{-}}\right. \\ \left. -\text{i}\kappa \left[2{{a}_{+}}a_{-}^{\ast}-\left(a_{+}^{\ast}{{a}_{+}}+a_{-}^{\ast}{{a}_{-}}\right)\right]\right\}, \end{array}\end{eqnarray} \tag{ 42 }$$

where ${{a}_{\pm }},a_{\pm }^{\ast}$ represent the complex photon field. Performing the basis rotation to classical and quantum fields as in equation (32) in section 2.1, and going to Fourier space, the action becomes

$$\begin{eqnarray}S={{{\int}^{}}_{\omega}}\left(a_{c}^{\ast}(\omega ),a_{q}^{\ast}(\omega )\right)\left(\begin{array}{*{35}{c}} 0 & {{P}^{A}}(\omega ) \\ {{P}^{R}}(\omega ) & {{P}^{K}} \end{array}\right)\left(\begin{array}{*{35}{l}} {{a}_{c}}(\omega ) \\ {{a}_{q}}(\omega ) \end{array}\right),\end{eqnarray} \tag{ 43 }$$

where we used the shorthand ${{{\int}^{}}_{\omega}}\equiv {\int}_{-\infty}^{\infty}\frac{\text{d}\omega}{2\pi}$. Furthermore, we have

$$\begin{eqnarray}&&{{P}^{R}}(\omega )={{P}^{A}}{{(\omega )}^{\ast}}=\omega -{{\omega}_{0}}+\text{i}\kappa,\quad {{P}^{K}}=2\text{i}\kappa.\end{eqnarray} \tag{ 44 }$$

P^R/A are the inverse retarded and advanced Green's functions, and P^K is the Keldysh component of the inverse Green's function. To see this, we evaluate the generating functional (33) by Gaussian integration (see appendix B). Then, the generating functional (35) for connected correlation functions is given by

$$\begin{eqnarray}W\left[{{J}_{c}},{{J}_{q}}\right]=-{{{\int}^{}}_{\omega}}\left(\,j_{q}^{\ast}(\omega ),j_{c}^{\ast}(\omega )\right)\left(\begin{array}{*{35}{c}} {{G}^{K}}(\omega ) & {{G}^{R}}(\omega ) \\ {{G}^{A}}(\omega ) & 0 \end{array}\right)\left(\begin{array}{*{35}{l}} {{j}_{q}}(\omega ) \\ {{j}_{c}}(\omega ) \end{array}\right).\end{eqnarray} \tag{ 45 }$$

According to equation (36), the second variation (i.e. the matrix in the above equation, see appendix A) represents the Green's function. It is obtained by inversion of the matrix in the action in equation (43),

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{G}^{R}}(\omega ) & ={{G}^{A}}{{(\omega )}^{\ast}}=\frac{1}{{{P}^{R}}(\omega )}=\frac{1}{\omega -{{\omega}_{0}}+\text{i}\kappa}, \end{array}\end{eqnarray} \tag{ 46 }$$

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{G}^{K}}(\omega ) & =-{{G}^{R}}(\omega ){{P}^{K}}{{G}^{A}}(\omega )=-\frac{i2\kappa}{{{\left(\omega -{{\omega}_{0}}\right)}^{2}}+{{\kappa}^{2}}}. \end{array}\end{eqnarray} \tag{ 47 }$$

We now summarize a few key structural properties that can be gleaned from this explicit discussion. Indeed, as we argue below, these properties are valid in general.

Conservation of probability—There is a zero matrix entry in equation (43), or, equivalently, in equation (45). Technically, as anticipated above, this property reflects a redundancy in the  ±  basis and eliminates it. This simplifies practical calculations in the Keldysh basis. Physically, this property ensures the normalization of the partition function ($Z=\text{tr}\rho (t)=1$), and can thus be interpreted as manifestation of the conservation of probability (${{\partial}_{t}}\text{tr}\rho (t)=0$), which is an exact property of physical problems. This can be seen as follows: consider the more general property, which implies the vanishing matrix element in the quadratic sector, $S\left[{{a}_{c}},a_{c}^{\ast},{{a}_{q}}=0,a_{q}^{\ast}=0\right]=0$. Any Keldysh action associated to the Liouville operator equation (16) has this property, as can be seen by setting ${{\psi}_{+}}={{\psi}_{-}}$ (i.e. ${{\phi}_{q}}=0$) in equation (28). Indeed, this operation on the Keldysh action may be interpreted as taking the trace in the operator based master equation (16): in this way, using the cyclic property of the trace allows us to shift all operators to one side of the density matrix, leading to the cancellation of terms such that ${{\partial}_{t}}\text{tr}\rho (t)=0$.

We still need to argue that the above property of the classical action also holds for the full theory, i.e. the effective action,

$$\begin{eqnarray}&&\Gamma\left[{{\bar{a}}_{c}},\bar{a}_{c}^{\ast},{{\bar{a}}_{q}}=0,\bar{a}_{q}^{\ast}=0\right]=0,\end{eqnarray} \tag{ 48 }$$

or, more schematically in the notation of section 2.1, $\Gamma\left[{{\rm{\bar \Phi}}_{c}},{{\rm{\bar \Phi}}_{q}}=0\right]=0$. To this end, we note that $W\left[{{J}_{c}},{{J}_{q}}=0\right]=0$ ($Z\left[{{J}_{c}},{{J}_{q}}=0\right]=1$) holds actually for arbitrary classical sources J_c (whereas in section 2.1 we worked additionally with J_c  =  0 for conceptual clarity): any term $\sim {{J}_{c}}$ can be absorbed into the underlying Hamiltonian, describing nothing but a Hamiltonian contribution in the presence of a classical external potential, and thus it cannot affect the normalization property of the theory. The above properties are equivalent:

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \Gamma\left[{{{\rm{\bar \Phi}}}_{c}},{{{\rm{\bar \Phi}}}_{q}}=0\right]=0\Leftrightarrow \,\,W\left[{{J}_{c}},{{J}_{q}}=0\right]=0. \end{array}\end{eqnarray} \tag{ 49 }$$

This is seen using the definition of the Legendre transform equation (38) and the definition of the quantum field in terms of equation (37) for one direction of the mutual implication. The other direction results from the involutory property of the Legendre transform¹⁴.

Sometimes, the property of conservation of probability—expressed in the effective action formalism as $\Gamma\left[{{\rm{\bar \Phi}}_{c}},{{\rm{\bar \Phi}}_{q}}=0\right]=0$—is referred to as 'causality structure' in the literature.

Hermeticity properties of the Green's functions—As can be read off from equations (46) and (47), ${{G}^{R}}(\omega )$ and ${{G}^{A}}(\omega )$ are Hermitian conjugates, and ${{G}^{K}}(\omega )$ is anti-Hermitian. These properties are exact, as can be seen from the definition of the Green's function in terms of functional derivatives in equation (34). Equivalently, these properties hold for the corresponding components of the inverse Green's function (see equation (41)).

Analytic structure—The poles of ${{G}^{R}}(\omega )$ are located at $\omega ={{\omega}_{0}}-\text{i}\kappa$, i.e. in the lower half of the complex plane (accordingly, the poles of ${{G}^{A}}(\omega )$ are in the upper half). In the real time domain, this implies that G^R describes the retarded response (and accordingly, G^A the advanced): indeed, taking the inverse Fourier transform, we obtain

$$\begin{eqnarray}&&{{G}^{R}}(t)=-\text{i}\theta (t){{\text{e}}^{-\left(\text{i}{{\omega}_{0}}+\kappa \right)t}},\end{eqnarray} \tag{ 50 }$$

where $\theta (t)$ is the Heaviside step function. Hence, the response of the system, if it is perturbed at t  =  0, is retarded (non-zero only for t  >  0) and decays.

The analytic structure of retarded and advanced Green's functions is a general property of Keldysh actions, too; one may then think of these Green's functions as renormalized, full single particle Green's functions connected to the bare, microscopic ones via renormalization, and thus preserving the analyticity properties. We note, however, that the pole of the retarded bosonic Green's function can approach the real axis from below via tuning of microscopic parameters. The touching point typically signals a physical instability: beyond that point, the description of the system must be modified qualitatively. Such a scenario is discussed in the next subsection.

Connection to the operator formalism—It is sometimes useful to restore the precise relation between the operator formalism and the functional integral description at the level of the Green's functions. At the single particle level, they read

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{G}^{R}}\left(t,{{t}^{\prime}}\right) & =-\text{i}\theta \left(t-{{t}^{\prime}}\right)\langle \left[a(t),{{a}^{\dagger}}\left({{t}^{\prime}}\right)\right]\rangle, \end{array}\end{eqnarray} \tag{ 51 }$$

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{G}^{K}}\left(t,{{t}^{\prime}}\right) & =-\text{i}\langle \left\{a(t),{{a}^{\dagger}}\left({{t}^{\prime}}\right)\right\}\rangle, \end{array}\end{eqnarray} \tag{ 52 }$$

and we note again that in stationary state ${{G}^{R/A/K}}\left(t,{{t}^{\prime}}\right)=$ ${{G}^{R/A/K}}\left(t-{{t}^{\prime}}\right)$. These relations are exact and can be obtained by going back to the  ±  basis: the retarded Green's function is

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{G}^{R}}\left(t,{{t}^{\prime}}\right) & =-\text{i}\langle {{a}_{c}}(t)a_{q}^{\ast}\left({{t}^{\prime}}\right)\rangle \\ {} & =-\frac{\text{i}}{2}\langle \left({{a}_{+}}(t)+{{a}_{-}}(t)\right)\left(a_{+}^{\ast}\left({{t}^{\prime}}\right)-a_{-}^{\ast}\left({{t}^{\prime}}\right)\right)\rangle \\ {} & =-\frac{\text{i}}{2}\left(\langle Ta(t){{a}^{\dagger}}\left({{t}^{\prime}}\right)\rangle +\langle \left[a(t),{{a}^{\dagger}}\left({{t}^{\prime}}\right)\right]\rangle -\langle \tilde{T}a(t){{a}^{\dagger}}\left({{t}^{\prime}}\right)\rangle \right) \\ {} & =-\text{i}\theta \left(t-{{t}^{\prime}}\right)\langle \left[a(t),{{a}^{\dagger}}\left({{t}^{\prime}}\right)\right]\rangle, \end{array}\end{eqnarray} \tag{ 53 }$$

where $T,\tilde{T}$ are the time-ordering, anti-time-ordering operators, which lead to a cancellation of the commutator for ${{t}^{\prime}}>t$ (see [102] for a more detailed discussion of time ordering on the Keldysh contour). The Keldysh Green's function in the operator formalism is obtained the same way,

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{G}^{K}}\left(t,{{t}^{\prime}}\right) & =-\text{i}\langle {{a}_{c}}(t)a_{c}^{\ast}\left({{t}^{\prime}}\right)\rangle \\ {} & =-\frac{\text{i}}{2}\langle \left({{a}_{+}}(t)+{{a}_{-}}(t)\right)\left(a_{+}^{\ast}\left({{t}^{\prime}}\right)+a_{-}^{\ast}\left({{t}^{\prime}}\right)\right)\rangle \\ {} & =-\frac{\text{i}}{2}\left(\langle Ta(t){{a}^{\dagger}}\left({{t}^{\prime}}\right)\rangle +\langle \left\{a(t),{{a}^{\dagger}}\left({{t}^{\prime}}\right)\right\}\rangle +\langle \tilde{{T}}\,a(t){{a}^{\dagger}}\left({{t}^{\prime}}\right)\rangle \right) \\ {} & =-\text{i}\langle \left\{a(t),{{a}^{\dagger}}\left({{t}^{\prime}}\right)\right\}\rangle. \end{array}\end{eqnarray} \tag{ 54 }$$

Response versus correlation functions—In order to generate response and correlation functions directly from the partition function, it is convenient to introduce source fields ${{j}_{+}},\,{{j}_{-}}$ and express the partition function as the average (see equation (31))

$$\begin{eqnarray}&&Z\left[\,{{j}_{+}},{{j}_{-}}\right]=\langle {{\text{e}}^{-\text{i}{{S}_{j}}}}\rangle ={\int}^{}\mathcal{D}\left[a_{+}^{\ast},a_{-}^{\ast},{{a}_{+}},{{a}_{-}}\right]{{\text{e}}^{\text{i}S-\text{i}{{S}_{j}}}},\end{eqnarray} \tag{ 55 }$$

where the source action is defined as

$$\begin{eqnarray}&&{{S}_{j}}={\int}_{t}\left(\,j_{+}^{\ast}{{a}_{+}}-j_{-}^{\ast}{{a}_{-}}+\text{c}.\text{c}.\right)={\int}_{t}\left(\,j_{c}^{\ast}{{a}_{q}}+j_{q}^{\ast}{{a}_{c}}+\text{c}.\text{c}.\right).\end{eqnarray} \tag{ 56 }$$

In the second step, we have performed a Keldysh rotation. Due to the normalization of the Keldysh path integral $Z\left[\,{{j}_{c}}=0,\,{{j}_{q}}=0\right]=1$, and as a consequence, expectation values of n-point functions of the fields a^*, a can be expressed via nth order functional derivatives of the partition function with respect to the source fields (see appendix A). For example

$$\begin{eqnarray}\begin{array}{*{35}{l}} \langle {{a}_{c}}(t)\rangle & =\text{i}\frac{\delta Z\left(\,{{j}_{c}},{{j}_{q}}\right)}{\delta j_{q}^{\ast}(t)}{{|}_{{{j}_{c}}={{j}_{q}}=0}}, \\ {{G}^{K}}\left(t,{{t}^{\prime}}\right) & =-\text{i}\langle {{a}_{c}}(t)a_{c}^{\ast}\left({{t}^{\prime}}\right)\rangle =\text{i}\frac{{{\delta}^{2}}Z\left(\,{{j}_{c}},{{j}_{q}}\right)}{\delta j_{q}^{\ast}(t)\delta {{j}_{q}}\left({{t}^{\prime}}\right)}{{|}_{{{j}_{c}}={{j}_{q}}=0}}, \\ {{G}^{R}}\left(t,{{t}^{\prime}}\right) & =-\text{i}\langle {{a}_{c}}(t)a_{q}^{\ast}\left({{t}^{\prime}}\right)\rangle =\text{i}\frac{{{\delta}^{2}}Z\left(\,{{j}_{c}},{{j}_{q}}\right)}{\delta j_{q}^{\ast}(t)\delta {{j}_{c}}\left({{t}^{\prime}}\right)}{{|}_{{{j}_{c}}={{j}_{q}}=0}}. \end{array}\end{eqnarray} \tag{ 57 }$$

Due to causality, the field j_q has to be zero in any physical setup and the introduction of this field only serves as a technical tool to compute expectation values via derivatives. On the other hand, the field j_c can, in principle, be different from zero and we will see in the following what the physical meaning of this classical source field is, and in which respect it generates the response function.

Responses: The retarded Green's function, often called synonymously response function, describes the linear response of a system which is perturbed by a weak external source field. As an illustrative example, consider a driven cavity system consisting of atoms and photons, which is considered to be in a stationary state. Due to imperfections in the cavity mirrors, there is a finite rate with which a photon is escaping the cavity, or an external photon is entering the cavity. This process is expressed by the Hamiltonian

$$\begin{eqnarray}&&{{H}_{p}}=\sqrt{2\kappa}\left({{b}^{\dagger}}a+{{a}^{\dagger}}b\right),\end{eqnarray} \tag{ 58 }$$

where the operators $b,{{b}^{\dagger}}$ represent photons outside the cavity. Shining a laser through the cavity mirror, the operators $b,{{b}^{\dagger}}$ can be replaced by the coherent laser field j(t), j^*(t), which oscillates with the laser frequency ${{\omega}_{j}}$. The corresponding Hamiltonian, describing the complete system, is (for the present purposes we can leave the Hamiltonian H of photons and atoms in the cavity unspecified)

$$\begin{eqnarray}&&{{H}_{j}}=H+\sqrt{2\kappa}\left(\,{{j}^{\ast}}(t)a+j(t){{a}^{\dagger}}\right).\end{eqnarray} \tag{ 59 }$$

Since the fields j, j^* are classical external fields, they are equal on the plus and minus contour and the Keldysh action in the presence of the laser is

$$\begin{eqnarray}&&{{S}_{j}}=S-{\int}_{t}\left(\,{{j}^{\ast}}(t){{a}_{q}}(t)+j(t)a_{q}^{\ast}(t)\right).\end{eqnarray} \tag{ 60 }$$

This action is very similar to the action in equation (56), with a linear source term, which is however ${{j}_{+}}={{j}_{-}}=j$. In the present example, the meaning of the source term is physically very transparent: it is nothing but the coherent laser field that is coupled to the cavity photons. For a weak source field, one is interested in the first order correction of observables induced by the coupling to the source. In the present case, this is the coherent light field inside the cavity $\langle {{a}_{c}}(t)\rangle$. Up to first order in j(t) it is given by

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{\langle {{a}_{c}}(t)\rangle}_{j}} & ={{\langle {{a}_{c}}(t)\rangle}_{j=0}}-\text{i}\sqrt{\kappa}{{{\int}^{}}_{{{t}^{\prime}}}}{{\langle {{a}_{c}}(t)a_{q}^{\ast}\left({{t}^{\prime}}\right)\rangle}_{j=0\,}}j\left({{t}^{\prime}}\right) \\ {} & ={{\langle {{a}_{c}}(t)\rangle}_{j=0}}+\sqrt{\kappa}{{{\int}^{}}_{{{t}^{\prime}}}}{{G}^{R}}\left(t-{{t}^{\prime}}\right)\,j\left({{t}^{\prime}}\right). \end{array}\end{eqnarray}$$

The retarded Green's function is therefore a measure of the system's response to an external perturbation. An experimental setup, with which one can measure the coherent light field and therefore the response function of the cavity via so-called homodyne detection, is illustrated in figure 4. A closely related function of interest is the spectral function

$$\begin{eqnarray}&&A(\omega )=-2\text{Im}\,{{G}^{R}}(\omega ).\end{eqnarray} \tag{ 61 }$$

Zoom In Zoom Out Reset image size

It is the distribution of excitation levels of the system, i.e. when adding a single photon with frequency ω to the system, $A(\omega )$ is the probability to hit the system at resonance. Indeed, it can be shown [178] that the spectral function is positive, $A(\omega )>0$, and fulfills the sum rule

$$\begin{eqnarray}&&{{{\int}^{}}_{\omega}}A(\omega )=\langle \left[a,{{a}^{\dagger}}\right]\rangle =1.\end{eqnarray} \tag{ 62 }$$

For our example of a single-mode cavity, the spectral function is given by

$$\begin{eqnarray}&&A(\omega )=\frac{2\kappa}{{{\left(\omega -{{\omega}_{0}}\right)}^{2}}+{{\kappa}^{2}}},\end{eqnarray} \tag{ 63 }$$

i.e. it is a Lorentzian, which is centered at the cavity frequency ${{\omega}_{0}}$ and has a half-width at half-maximum given by κ. Note that for $\kappa \to 0$, the photon number states become exact eigenstates and the spectral density reduces to a δ-function peaked at ${{\omega}_{0}}$, $A(\omega )=2\pi \delta \left(\omega -{{\omega}_{0}}\right).$ As these considerations illustrate, the retarded Green's function contains essential information on the system's response towards an external perturbation, and on the spectral properties.

Correlations: The Keldysh Green's function contains elementary information on the system's correlations and the occupation of the individual quantum mechanical modes. A prominent example of a correlation function in quantum optics is the photonic g⁽²⁾ correlation function. It is defined as the four-point correlator

$$\begin{eqnarray}&&{{g}^{(2)}}(t,\tau )=\frac{\langle {{a}^{\dagger}}(t){{a}^{\dagger}}(t+\tau )a(t+\tau )a(t)\rangle}{|\langle {{a}^{\dagger}}(t)a(t)\rangle {{|}^{2}}},\end{eqnarray} \tag{ 64 }$$

and it is proportional to the intensity fluctuations of the intracavity radiation field ${{g}^{(2)}}(t,\tau )\propto \langle I(t)I(t+\tau )\rangle$. In the limit of $\tau \to 0$, the g⁽²⁾ correlation function reveals the statistics of the cavity photons, i.e. it demonstrates super-Poissonian (g⁽²⁾(0)  >  1) or sub-Poissonian statistics (g⁽²⁾(0)  <  1) as an effect of the light-matter interactions. In the absence of interactions (as in our example), it is straightforward to show that

$$\begin{eqnarray}&&{{g}^{(2)}}(t,\tau )=1+\frac{|{{G}^{K}}(t,\tau +t)+{{G}^{R}}(t,\tau +t)-{{G}^{A}}(t,\tau +t){{|}^{2}}}{|\text{i}{{G}^{K}}(t,t)-1{{|}^{2}}}.\end{eqnarray} \tag{ 65 }$$

For equal times, the retarded and advanced Green's functions satisfy ${{G}^{R}}(t,t)-{{G}^{A}}(t,t)=-i$, indicating the expected photon-bunching at $\tau \to 0$. On the other hand, in the presence of interactions, equation (64) is modified and contains additional higher order terms, as well as off-diagonal contributions. However, it remains a measure of the photonic statistics in the cavity due to the physical relation to the intensity fluctuations.

A particularly relevant and instructive limit of the two-time Keldysh Green's function equation (52) concerns equal times $t={{t}^{\prime}}$, which in general describes static correlation functions, or covariances in the quantum optics language. In the cavity example, it yields the mode occupation number,

$$\begin{eqnarray}&&\text{i}{{G}^{K}}(t,t)=2\langle {{a}^{\dagger}}a\rangle +1.\end{eqnarray} \tag{ 66 }$$

The appearance of the combination $2\langle {{a}^{\dagger}}a\rangle +1$ is rather intuitive when taking into account the relation to the operatorial formalism. Indeed, in operator language, $2{{a}^{\dagger}}a+1=\left\{a,{{a}^{\dagger}}\right\}=\left\{{{a}^{\dagger}},a\right\}$ is invariant under permutation of the operators $a,{{a}^{\dagger}}$, and therefore lends itself to a direct functional integral representation, which in fact carries no information on operator ordering¹⁵. At the same time, the anticommutator carries all the physical information on state occupation, while the commutator carries none, $\left[a,{{a}^{\dagger}}\right]=1$ (note ${{a}^{\dagger}}a=\frac{1}{2}\left(\left\{a,{{a}^{\dagger}}\right\}-\left[a,{{a}^{\dagger}}\right]\right)$).

Then, the explicit form of the Keldysh Green's function stated in equation (47) leads to the result

$$\begin{eqnarray}&&\langle {{a}^{\dagger}}a\rangle =\frac{1}{2}\left(\text{i}{{{\int}^{}}_{\omega}}{{G}^{K}}(\omega )-1\right)=0,\end{eqnarray} \tag{ 67 }$$

showing that the cavity is empty in the steady state—which is not surprising in the absence of external pumping.

We note that in the master equation formalism, information on the static or spatial correlations is most easily accessible—only the state (density matrix) has to be known, but not the dynamics acting on it. Temporal correlation functions can be extracted using the quantum regression theorem [186]. In the Keldysh formalism, spatial and temporal correlation functions are treated on an equal footing.

We can concisely summarize the above discussion of the response and Keldysh Green's functions of the system, at the risk of oversimplifying:

Relation to thermodynamic equilibrium—Before we move on, let us briefly discuss how we have to modify the formalism in order to describe a cavity in thermodynamic equilibrium. In general, a system is in thermodynamic equilibrium at a temperature $T=1/\beta$, if it is in a Gibbs state $\rho ={{\text{e}}^{-\beta H}}/Z$ where $Z=\text{tr}{{\text{e}}^{-\beta H}}$, and its dynamics is coherent and generated by the Hamiltonian H. (In particular, dissipative dynamics described by a term in the Lindblad form in equation (16) is not compatible with equilibrium conditions: see [187–189], and the discussion in section 2.4.1.) Specifying the thermal density matrix at t₀ in the Keldysh functional integral in equation (25) explicitly in terms of its matrix elements is rather inconvenient, especially if one is interested in steady state properties and wants to take the limit ${{t}_{0}}\to -\infty.$ An alternative is suggested by the observation that thermodynamic equilibrium can be established in the system if it is weakly coupled to a thermal bath. Then, the finite decay rate κ in the retarded Green's function in equation (46) should be replaced by an infinitesimally weak system–bath coupling $\delta \to 0$. Additionally, the Keldysh Green's function has to be modified in such a way that the integral in equation (67) yields the Bose distribution function $n\left({{\omega}_{0}}\right)$ implying thermal occupations of the cavity mode,

$$\begin{eqnarray}&&\langle {{a}^{\dagger}}a\rangle =n\left({{\omega}_{0}}\right)=\frac{1}{{{\text{e}}^{\beta {{\omega}_{0}}}}-1}.\end{eqnarray} \tag{ 68 }$$

This can be achieved by replacing the Keldysh component of the inverse Green's function in equation (43) by ${{P}^{K}}(\omega )=\text{i}2\delta \left(2n(\omega )+1\right)$. We emphasize the key structural difference of the thermal Keldysh component, which is strongly frequency dependent, to the Markovian case discussed previously, where this entry is frequency independent—this gives a strong hint that Markovian systems can behave quite differently from systems in thermal equilibrium. With ${{P}^{R}}(\omega )={{P}^{A}}{{(\omega )}^{\ast}}=\omega -{{\omega}_{0}}-\text{i}\delta$ we obtain the equilibrium Green's functions

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{G}^{R}}(\omega ) & ={{G}^{A}}{{(\omega )}^{\ast}}=\frac{1}{\omega -{{\omega}_{0}}+\text{i}\delta}, \end{array}\end{eqnarray} \tag{ 69 }$$

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{G}^{K}}(\omega ) & =-\text{i}2\pi \delta \left(\omega -{{\omega}_{0}}\right)\left(2n(\omega )+1\right). \end{array}\end{eqnarray} \tag{ 70 }$$

Note that the Green's functions obey a thermal fluctuation-dissipation relation (FDR), which for the present example of a single bosonic mode reads

$$\begin{eqnarray}&&{{G}^{K}}(\omega )=\left(2n(\omega )+1\right)\left({{G}^{R}}(\omega )-{{G}^{A}}(\omega )\right).\end{eqnarray} \tag{ 71 }$$

This is discussed in more detail in section 2.4.1. For the time being, let us mention that this construction to describe thermodynamic equilibrium by adding infinitesimal dissipative terms not only works in the present case of a quadratic action but can also be applied in the interacting case. Then, the construction ensures that the free Green's functions (i.e. the ones obtained by ignoring the interactions) obey an FDR. If the non-linear terms in the action obey the equilibrium symmetry discussed in section 2.4.1 (which is the case for generic interaction terms), this property is shared by the full Green's functions of the non-linear system.

General Fluctuation–Dissipation Relations—While equation (71) is valid only in thermodynamic equilibrium, it is always possible to parameterize the (anti-Hermitian) Keldysh Green's function in terms of the retarded and advanced Green's functions (which, as discussed above, are Hermitian conjugates of each other) and a Hermitian matrix $F={{F}^{\dagger}}$ in the form [102, 178]

$$\begin{eqnarray}&&{{G}^{K}}={{G}^{R}}\circ F-F\circ {{G}^{A}},\end{eqnarray} \tag{ 72 }$$

where ∘ denotes convolution. In this parametrization, F is the distribution function, which describes the distribution of (quasi-) particles over the modes of the system. For a non-equilibrium steady state, F is time-translational invariant and its Fourier transform in frequency space $F(\omega )$ represents the energy resolved occupation of (quasi-) particle modes. On the other hand, as we discuss in detail in section 5.2, for the case of a time-evolving system, for which time translation invariance is absent, the Wigner transform of F corresponds to the instantaneous local distribution function. For the important case of thermodynamic equilibrium (of a bosonic system), it is $F(\omega )=\coth (\beta \omega /2)=2n(\omega )+1$, and equation (72) reduces to (71). For the case where the bosonic Green's function has a matrix structure in Nambu space, a subtlety concering the preservation of the symplectic structure of that space arises, see section 5.2. There, also, a time-dependent variant of the non-equilibrium fluctuation-dissipation relations is discussed.

In the previous section, we discussed the simple case of a quadratic Keldysh action, which allowed us to perform the Keldysh functional integral explicitly. An additional simplification resulted from the fact that we were considering only a single bosonic mode. Let us now consider a genuine many-body problem, which is non-linear and in which the system consists of a continuum of modes. To be specific, in this section we discuss the model introduced in section 1.3.2 for a bosonic many-body system with interactions and non-linear loss processes (i.e. with a loss rate that is proportional to the density) in addition to the linear dissipative terms which were already present in the example of the single-mode cavity. Then, for a specific value of the mean-field density ${{\rho}_{0}}$, non-linear loss and linear pump exactly balance each other and the system reaches a stationary state. If ${{\rho}_{0}}$ is different from zero, this signals the presence of a condensate, which is accompanied by the breaking of a specific phase rotation symmetry as will be discussed in detail in section 2.4, and the establishment of long-range order. In section 4, we give a detailed account of the influence of fluctuations on this driven–dissipative condensation transition. Here, we content ourselves with a mean-field analysis, which serves to illustrate some of the field theoretical concepts introduced in section 2.1—in particular, the effective action and field equations—in a simple setting.

In the basis of classical and quantum fields, the Keldysh action associated with the quantum master equation (9) reads

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} S={{{\int}^{}}_{t,\mathbf{x}}}\left\{\phi _{q}^{\ast}\left(\text{i}{{\partial}_{t}}+{{K}_{c}}{{\nabla}^{2}}-{{r}_{c}}\,+\,\text{i}{{r}_{d}}\right){{\phi}_{c}}+\text{c}.\text{c}.\right. \\ -\left[\left({{u}_{c}}-\text{i}{{u}_{d}}\right)\left(\phi _{q}^{\ast}\phi _{c}^{\ast}\phi _{c}^{2}+\phi _{q}^{\ast}\phi _{c}^{\ast}\phi _{q}^{2}\right)+\text{c}.\text{c}.\right] \\ \left. +\,\text{i}2\left(\gamma +2{{u}_{d}}\phi _{c}^{\ast}{{\phi}_{c}}\right)\phi _{q}^{\ast}{{\phi}_{q}}\right\}, \end{array}\end{eqnarray} \tag{ 73 }$$

where ${{K}_{c}}=1/\left(2{{m}_{\text{LP}}}\right)$ and ${{r}_{c}}=\omega _{\text{LP}}^{0}$; as additional parameters, we introduced the noise level $\gamma =\left({{\gamma}_{l}}+{{\gamma}_{p}}\right)/2$ and the spectral mass or gap ${{r}_{d}}=\left({{\gamma}_{l}}-{{\gamma}_{p}}\right)/2$. Hence, the rates of losses and pumping add up to the total noise level; in contrast, the difference of these rates enters in the spectral gap r_d, which becomes negative when the rate of incoherent pumping exceeds the single-particle loss rate, signaling the physical instability against condensation. In a mean-field analysis of the condensation transition, we perform a saddle-point approximation of the functional integral in equation (39). To leading order, fluctuations around the field expectation values are completely neglected. The expectation values are then obtained as spatially homogenous and stationary solutions to the classical field equations (here, 'classical' refers to the fact that these field equations are derived from the classical (or: bare, microscopic) action S discarding fluctuations, in contrast to the field equations in equation (40), which involve the effective action)

$$\begin{eqnarray}&&\frac{\delta S}{\delta \phi _{c}^{\ast}}=0,\quad \frac{\delta S}{\delta \phi _{q}^{\ast}}=0.\end{eqnarray} \tag{ 74 }$$

As already mentioned above in section 2.2.1, there is no term in the action equation (38) with zero power of either $\phi _{q}^{\ast}$ or ${{\phi}_{q}}$, and the same is clearly true for $\delta S/\delta \phi _{c}^{\ast}$. Therefore, the first equation in (74) is solved by ${{\phi}_{q}}=0$. Inserting this condition into the second equation, we have

$$\begin{eqnarray}&&\left[-{{r}_{c}}+\text{i}{{r}_{d}}-\left({{u}_{c}}-\text{i}{{u}_{d}}\right)|{{\phi}_{0}}{{|}^{2}}\right]{{\phi}_{0}}=0.\end{eqnarray} \tag{ 75 }$$

The solution ${{\phi}_{c}}={{\phi}_{0}}$ is determined by the imaginary part of equation (75): for ${{r}_{d}}\geqslant 0$, in the so-called symmetric phase, the classical field expectation value is zero, ${{\rho}_{0}}=|{{\phi}_{0}}{{|}^{2}}=0$, whereas for r_d  <  0 we have a finite condensate density ${{\rho}_{0}}=-{{r}_{d}}/{{u}_{d}}$. Taking the real part of equation (75), we obtain for the parameter r_c the relation ${{r}_{c}}={{u}_{c}}{{r}_{d}}/{{u}_{d}}$. This condition can always be satisfied by proper choice of a rotating frame, i.e. by performing a gauge transformation ${{\phi}_{c}}\mapsto {{\phi}_{c}}{{\text{e}}^{-\text{i}\omega t}}$ such that ${{r}_{c}}\mapsto {{r}_{c}}-\omega$. In the original (laboratory) frame this simply means that the condensate amplitude oscillates at a finite frequency.

In a first step beyond mean field theory, quadratic fluctuations around the mean-field order parameter can be investigated within a Bogoliubov or tree-level expansion: we set ${{\phi}_{c}}={{\phi}_{0}}+\delta {{\phi}_{c}},{{\phi}_{q}}=\delta {{\phi}_{q}}$ in the action equation (73) and expand the resulting expression to second order in the fluctuations $\delta {{\phi}_{c,q}}$. The inverse retarded, advanced and Keldysh Green's functions now become $2\times 2$ matrices in the space of Nambu spinors $\delta {{\Phi}_{\nu}}=\left(\delta {{\phi}_{\nu}},\delta \phi _{\nu}^{\ast}\right)$. In particular, we have in the frequency and momentum domain

$$\begin{eqnarray}{{P}^{R}}\left(\omega,\mathbf{q}\right)=\left(\begin{array}{*{35}{c}} \omega -{{K}_{c}}{{q}^{2}}-\left({{u}_{c}}-\text{i}{{u}_{d}}\right){{\rho}_{0}} & -\left({{u}_{c}}-\text{i}{{u}_{d}}\right){{\rho}_{0}} \\ -\left({{u}_{c}}+\text{i}{{u}_{d}}\right){{\rho}_{0}} & -\omega -{{K}_{c}}{{q}^{2}}-\left({{u}_{c}}+\text{i}{{u}_{d}}\right) \end{array}\right),\\ &&\left. {{P}^{A}}\left(\omega,\mathbf{q}\right)\right)={{P}^{R}}{{\left(\omega,\mathbf{q}\right)}^{\dagger}},\\ &&{{P}^{K}}=\text{i}\gamma {\mathbb 1}.\end{eqnarray} \tag{ 76 }$$

The excitation spectrum is obtained from the condition $\det {{P}^{R}}\left(\omega,\mathbf{q}\right)=0$. Indeed, this is the condition for the field equation of the fluctuations ${{P}^{R}}\left(\omega,\mathbf{q}\right)\delta {{\Phi}_{c}}\left(\omega,\mathbf{q}\right)=0$ to have nontrivial solutions. This yields [159]

$$\begin{eqnarray}&&\omega _{1,2}^{R}=-\text{i}{{u}_{d}}{{\rho}_{0}}\pm \sqrt{{{K}_{c}}{{q}^{2}}\left({{K}_{c}}{{q}^{2}}+2{{u}_{c}}{{\rho}_{0}}\right)-{{\left({{u}_{d}}{{\rho}_{0}}\right)}^{2}}}.\end{eqnarray} \tag{ 77 }$$

We note that due to the tree-level shifts $\propto {{\rho}_{0}}$ the above-described instability for r_d  <  0 is lifted: both poles are consistently located in the lower complex half-plane, indicating a physically stable situation with decaying single-particle excitations. For u_d  =  0, equation (77) reduces to the standard Bogoliubov result [190], where for $q\to 0$ the dispersion is phononic, $\omega _{1,2}^{R}=\pm cq$ ($c=\sqrt{2{{K}_{c}}{{u}_{c}}{{\rho}_{0}}}$ the speed of sound), whereas particle-like behavior $\omega _{1,2}^{R}\sim {{K}_{c}}{{q}^{2}}$ is obtained at high momenta. Here, due to the presence of two-body loss ${{u}_{d}}\ne 0$, the dispersion is qualitatively modified: while at high momenta the dominant behavior is still $\omega _{1,2}^{R}\sim {{K}_{c}}{{q}^{2}}$, at low momenta we find purely incoherent diffusive, non-propagating modes $\omega _{1}^{R}\sim -\text{i}\frac{{{K}_{c}}{{u}_{c}}}{{{u}_{d}}}{{q}^{2}}$ and $\omega _{2}^{R}\sim -\text{i}2{{u}_{d}}{{\rho}_{0}}$. In particular, for q  =  0 we have $\omega _{1}^{R}=0$: this is a dissipative Goldstone mode [159, 191, 192], associated with the spontaneous breaking of the global U(1) symmetry in the ordered phase. The existence of such a mode is not bound to the mean-field approximation, but rather is an exact property guaranteed by the U(1) invariance of the effective action, even in the present case of a driven–dissipative condensate. We will come back to this point in section 2.4.6.

A system is in thermodynamic equilibrium at a temperature $T=1/\beta$, if (i) the density matrix is given by $\rho ={{\text{e}}^{-\beta H}}/Z$ where $Z=\text{tr}{{\text{e}}^{-\beta H}}$, and (ii) the very same Hamiltonian operator H appearing in ρ generates the unitary time evolution of the system, $U={{\text{e}}^{-\text{i}Ht}}$. Condition (ii) implies that static correlations are in general not sufficient to prove that a system is in thermal equilibrium; however, dynamical correlations are: this can be inferred from the fact that the static correlations of a physical system can always be encoded in a density matrix ρ, and the latter can always be parameterized formally as an equilibrium density matrix, $\rho ={{\text{e}}^{-\beta {{H}^{\prime}}}}/{{Z}^{\prime}}$, with a Hermitian operator ${{H}^{\prime}}$. (In other words, any state can be thought of as a thermal state with respect to some Hamiltonian.) On the other hand, static (i.e. purely momentum- or space-dependent) properties do not allow us to discriminate whether the generator of dynamics coincides with ${{H}^{\prime}}$. In sharp contrast, any dynamical (i.e. frequency- or time-dependent) observable is manifestly governed by this generator, as is easily seen in the Heisenberg picture (or a suitable generalization to open systems). Response functions at finite frequency are such dynamical observables, and the equilibrium conditions formulated above are reflected in the fact that in thermodynamic equilibrium the response of the system to an external perturbation at a frequency ω is related to thermal fluctuations within the system at the same frequency by what is known as a fluctuation-dissipation relation [200] (FDR; for a discussion in the context of non-equilibrium Bose–Einstein condensation see [201]). Such a relation, which is equivalent to the combination of the Kubo–Martin–Schwinger (KMS) condition [202, 203] with time reversal [204], is valid for any pair of operators. In particular, for the case that these are the basic field operators of a single-component Bose system at the temperature $T=1/\beta$, the FDR reads (note that this is a generalization of equation (71) to the case of a spatial continuum of degrees of freedom)

$$\begin{eqnarray}&&{{G}^{K}}\left(\omega,\mathbf{q}\right)=\coth \,\left(\frac{\omega}{2T}\right)\left({{G}^{R}}\left(\omega,\mathbf{q}\right)-{{G}^{A}}\left(\omega,\mathbf{q}\right)\right)\end{eqnarray} \tag{ 84 }$$

(in the presence of a chemical potential μ, in the argument of the trigonometric function we have to shift $\omega \to \omega -\mu$). In quantum field theory, relations among Green's functions often follow as consequences of a symmetry of the action. Then, they are known as Ward–Takahashi identities associated with the symmetry [86, 205]¹⁶. This raises the question, whether FDRs and hence the presence of thermodynamic equilibrium conditions are also connected to a symmetry of the Keldysh action.

To make the point clear, let us rephrase the question: above we made the observation that the defining property of (canonical) thermal equilibrium is the validity of FDRs. Importantly, these FDRs hold for correlation functions of arbitrary order, i.e. not just the two-point Keldysh Green's function in equation (84) can be expressed through response functions, but also any higher order correlation function is determined by corresponding higher order response functions. Then, the question—which we answer in the affirmative below—is, whether this infinite hierarchy of FDRs expressing thermal equilibrium is related to a structural property of the theory. Indeed, this structural property is a symmetry of the Keldysh action, i.e. a transformation of the fields $\Psi\mapsto {{\mathcal{T}}_{\beta}}\Psi$ such that

$$\begin{eqnarray}&&S[\Psi]=\tilde{S}\left[{{\mathcal{T}}_{\beta}}\Psi\right].\end{eqnarray} \tag{ 85 }$$

Here, we denote $\Psi=\left({{\psi}_{+}},\psi _{+}^{\ast},{{\psi}_{-}},\psi _{-}^{\ast}\right)$, and we specify the precise form of ${{\mathcal{T}}_{\beta}}$ below in equation (88). The tilde on the RHS of the equation indicates that all external fields appearing in S have to be replaced by their corresponding time-reversed values. To name an example, the signs of magnetic fields have to be inverted [189]. Evidently, discussing a single symmetry instead of an infinite hierarchy of equations is much more elegant and practical: checking whether a given Keldysh action obeys equation (85) is straightforward and can be accomplished in finite time, in contrast to verifying the full set of FDRs.

It is easily seen how these FDRs can be deduced as a consequence of the symmetry of the Keldysh action (85). To this end, we note that by a Keldysh rotation (32) and after Fourier transformation, the Green's functions appearing in equation (84) can be expressed as a sum of averages of the form $\langle {{\psi}_{\sigma}}\left(t,\mathbf{x}\right)\psi _{{{\sigma}^{\prime}}}^{\ast}\left({{t}^{\prime}},{{\mathbf{x}}^{\prime}}\right)\rangle$. Writing these explicitly as field integrals, and anticipating that a change of integration variables $\Psi\to {{\mathcal{T}}_{\beta}}\Psi$ leaves the functional measure $\mathcal{D}[\Psi]$ invariant [189], we find

$$\begin{eqnarray}\begin{array}{*{35}{l}} \langle {{\psi}_{\sigma}}\left(t,\mathbf{x}\right)\psi _{{{\sigma}^{\prime}}}^{\ast}\left({{t}^{\prime}},{{\mathbf{x}}^{\prime}}\right)\rangle & ={\int}^{}\mathcal{D}[\Psi]\,{{\psi}_{\sigma}}\left(t,\mathbf{x}\right)\psi _{{{\sigma}^{\prime}}}^{\ast}\left({{t}^{\prime}},{{\mathbf{x}}^{\prime}}\right){{\text{e}}^{\text{i}S[\Psi]}} \\ {} & ={\int}^{}\mathcal{D}\left[{{\mathcal{T}}_{\beta}}\Psi\right]\,{{\mathcal{T}}_{\beta}}{{\psi}_{\sigma}}\left(t,\mathbf{x}\right){{\mathcal{T}}_{\beta}}\psi _{{{\sigma}^{\prime}}}^{\ast}\left({{t}^{\prime}},{{\mathbf{x}}^{\prime}}\right){{\text{e}}^{\text{i}S\left[{{\mathcal{T}}_{\beta}}\Psi\right]}} \\ {} & =\langle {{\mathcal{T}}_{\beta}}{{\psi}_{\sigma}}\left(t,\mathbf{x}\right){{\mathcal{T}}_{\beta}}\psi _{{{\sigma}^{\prime}}}^{\ast}\left({{t}^{\prime}},{{\mathbf{x}}^{\prime}}\right)\rangle. \end{array}\end{eqnarray} \tag{ 86 }$$

In the last equality, we used the symmetry of the Keldysh action (85) (assuming for simplicity that there are no external fields). Inserting here the explicit form of the symmetry transformation specified in equation (88) below where $\beta =1/T$ is the inverse temperature, it can be seen that equation (86) is in fact equivalent to the FDR (84) [189].

By generalizing this argument to arbitrary field averages, one can establish the full equivalence between the infinite hierarchy of FDRs and the symmetry property of the Keldysh action (85) [189]. Thus, the symmetry of the Keldysh action under this transformation is a direct proof of the presence of thermodynamic equilibrium conditions. The existence of a symmetry that is related to thermal equilibrium has first been realized in the context of classical stochastic models [197, 198], [206–208], and these considerations have been extended to the realm of quantum systems in [189, 209].

What is the explicit form of the transformation ${{\mathcal{T}}_{\beta}}$? We can guess its essential ingredients by reminding ourselves that, as stated above equation (84), FDRs can be obtained from the KMS condition [202, 203]. The latter reads, for operators $A(t)={{\text{e}}^{\text{i}Ht}}A(0){{\text{e}}^{-\text{i}Ht}}$ and $B\left({{t}^{\prime}}\right)$ in the Heisenberg representation,

$$\begin{eqnarray}&&\langle A(t)B\left({{t}^{\prime}}\right)\rangle =\langle B\left({{t}^{\prime}}\right)A\left(t+\text{i}\beta \right)\rangle.\end{eqnarray} \tag{ 87 }$$

Comparing this with equation (86) indicates that ${{\mathcal{T}}_{\beta}}$ involves a translation of t into the complex plane by an amount proportional to the inverse temperature β. Moreover, the order of operators on the RHS of the KMS condition is reversed as compared to the LHS. The original time order can be restored by means of a time reversal transformation [210]—this step is necessary in order to obtain a time-ordered expression which can be expressed as a Keldysh field integral (by construction, field integrals yield time-ordered averages [180]). A careful analysis [189] leads to the precise form of the thermal symmetry transformation ($\sigma =+(-)$ for fields on the forward (backward) branch):

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{\mathcal{T}}_{\beta}}{{\psi}_{\sigma}}\left(t,\mathbf{x}\right) & =\psi _{\sigma}^{\ast}\left(-t+\text{i}\sigma \beta /2,\mathbf{x}\right), \\ {{\mathcal{T}}_{\beta}}\psi _{\sigma}^{\ast}\left(t,\mathbf{x}\right) & ={{\psi}_{\sigma}}\left(-t+\text{i}\sigma \beta /2,\mathbf{x}\right) \end{array}\end{eqnarray} \tag{ 88 }$$

(in the presence of a chemical potential μ, we have to multiply the RHS of the first (second) line by ${{\text{e}}^{\sigma \beta \mu /2}}\left({{\text{e}}^{-\sigma \beta \mu /2}}\right)$). The transformation ${{\mathcal{T}}_{\beta}}$ is a composition of complex conjugation of the fields ${{\psi}_{\sigma}}$ and inversion of the sign of the time t—both originating from the time reversal transformation—and translation of the value of t by an amount $i\sigma \beta /2$. The latter is induced by the KMS condition equation (87). We note that the translation of t in equation (88) takes opposite signs depending on whether a field on the forward or on the backward branch is being transformed. As we show in section 2.4.4 below, a similar form of time translations is connected to conservation of energy: indeed, if the Keldysh action is invariant under time translations ${{\psi}_{\sigma}}\left(t,\mathbf{x}\right)\mapsto {{\psi}_{\sigma}}\left(t+\sigma s,\mathbf{x}\right)$ with $s\in \mathbb{R}$, the total energy in the system is conserved. A crucial difference from the time translation that is part of ${{\mathcal{T}}_{\beta}}$ is that energy conservation requires invariance under shifts by an arbitrary real value s, whereas ${{\mathcal{T}}_{\beta}}$ involves a shift by the purely imaginary value $i\beta /2,$ where $\beta =1/T$ is fixed and determined by the temperature.

Under which conditions does the Keldysh action (30) with $\mathcal{L}$ defined in equation (28) have the thermal symmetry, i.e. under which conditions does it describe a system in thermal equilibrium? Let us first consider the parts of the action corresponding to unitary time evolution, i.e. the first two terms in equation (30) and the first line in equation (28). It is straightforward to check [189] that these terms are symmetric if the Hamiltonian densities ${{H}_{\pm }}$ do not explicitly depend on time. On the other hand, adding an external classical driving field such as a laser, and thus breaking time translational invariance by making ${{H}_{\pm }}$ time dependent, the thermal symmetry is also broken¹⁷. The violation of the thermal symmetry by classical driving fields on the level of a microscopic Hamiltonian description indicates that quantum master equations correspond to genuine non-equilibrium conditions [187–189]. Physically, this is due to the fact that a system for which such a description is appropriate is necessarily driven.

In an effective description of a driven–dissipative system in terms of a Markovian master equation in a rotating frame, there is often no explicit time dependence—indeed, our derivation of the dissipative Keldysh action in section 2.1 started from equation (16) with time-independent Hamiltonian. Even then, the thermal symmetry can be used to diagnose non-equilibrium conditions [189]. Again, it is sufficient to study only the time-translation part of ${{\mathcal{T}}_{\beta}}$: for time-independent ${{H}_{\pm }}$ and ${{L}_{\alpha,\pm }}$ in equation (28), the Keldysh action is still invariant under time translations of the form ${{\psi}_{\sigma}}\left(t,\mathbf{x}\right)\mapsto {{\psi}_{\sigma}}\left(t+s,\mathbf{x}\right)$. Importantly, this differs from the time translation that occurs in ${{\mathcal{T}}_{\beta}}$ by the absence of a factor of σ, meaning that t is shifted by the same amount on the forward and backward branches. Then, by means of a simple shift of the integration variable t in equation (30) the original form of the action can be restored. This strategy fails for the dissipative contributions in equation (28) that couple forward and backward branches, when the transformation involves the contour index σ¹⁸. Again we reach the conclusion that quantum master equations describe genuine non-equilibrium conditions, even though by a slightly different argument than in the driven but purely Hamiltonian setting.

In some cases, the Markov and rotating-wave approximations leading to a description in terms of a quantum master equation or equivalent Keldysh action are applied in the absence of external driving fields. Then, these approximations might still be justified to study the behavior of specific observables, even though they explicitly break the symmetry [189]. To give an example, if one attempts to study thermalization of a system due to the coupling with a heat bath, and one integrates out the bath using the abovementioned approximations, the resulting dynamics of the system will still lead to a thermal stationary state. Correspondingly, all static properties of the system will appear thermal. However, as discussed at the beginning of the present section, to unambiguously prove the presence of thermal equilibrium conditions, one has to consider dynamical signatures such as FDRs. Then, as a consequence of the explicit violation of the thermal symmetry by the Markov and rotating-wave approximations, the description of the system dynamics in terms of a quantum master equation will lead to the wrong prediction that fluctuations in the system do not obey an FDR.

A simple example that illustrates the above discussion is given by a single bosonic mode a with energy ${{\omega}_{0}}$, driven coherently at a frequency ${{\omega}_{l}}$, and coupled to a bath of harmonic oscillators ${{b}_{\mu}}$ in thermal equilibrium. The associated Keldysh action can be decomposed as $S={{S}_{s}}+{{S}_{sb}}+{{S}_{b}}$, where the action for the system consists of two parts, ${{S}_{s}}={{S}_{0}}+{{S}_{l}}$ which read

$$\begin{eqnarray}&&{{S}_{0}}=\underset{\sigma}{\sum}\,\sigma {\int}_{\omega}a_{\sigma}^{\ast}(\omega )\left(\omega -{{\omega}_{0}}\right){{a}_{\sigma}}(\omega ),\end{eqnarray} \tag{ 89 }$$

and

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{S}_{l}} & =\Omega \underset{\sigma}{\sum}\,\sigma {\int}^{}\text{d}t\left({{a}_{\sigma}}(t){{\text{e}}^{\text{i}{{\omega}_{l}}t}}+a_{\sigma}^{\ast}(t){{\text{e}}^{-\text{i}{{\omega}_{l}}t}}\right) \\ {} & =\Omega \underset{\sigma}{\sum}\,\sigma \left({{a}_{\sigma}}\left({{\omega}_{l}}\right)+a_{\sigma}^{\ast}\left({{\omega}_{l}}\right)\right). \end{array}\end{eqnarray} \tag{ 90 }$$

$\Omega$ is the amplitude of the driving field, and due to the harmonic time dependence of the drive it affects only the component $a\left({{\omega}_{l}}\right)$ in frequency space. The action of the bath is expressed most conveniently on the basis of classical and quantum fields (see equation (32)). In addition to the coherent part stemming from the oscillator frequencies, it involves infinitesimal dissipative regularization terms specifying the thermal equilibrium state of the bath at inverse temperature $\beta =1/T$ as discussed in section 2.2.1,

$$\begin{eqnarray}&&{{S}_{b}}=\underset{\mu}{\sum}\,{{{\int}^{}}_{\omega}}\left(\begin{array}{*{35}{c}} b_{\mu,c}^{\ast}(\omega ),b_{\mu,q}^{\ast}(\omega ) \end{array}\right)\\ \,\,\,\,\times \,\left(\begin{array}{*{35}{l}} 0 & \omega -{{\omega}_{\mu}}-\text{i}\delta \\ \omega -{{\omega}_{\mu}}+\text{i}\delta & \text{i}2\delta \coth (\beta \omega /2) \end{array}\right)\,\left(\begin{array}{*{35}{l}} {{b}_{\mu,c}}(\omega ) \\ {{b}_{\mu,q}}(\omega ) \end{array}\right).\end{eqnarray} \tag{ 91 }$$

Finally, the system–bath interaction with coupling strength λ corresponds to the following contribution to the action:

$$\begin{eqnarray}&&{{S}_{sb}}=\lambda \underset{\sigma}{\sum}\,\sigma \underset{\mu}{\sum}\,{\int}_{\omega}\left(a_{\sigma}^{\ast}(\omega ){{b}_{\mu,\sigma}}(\omega )+{{a}_{\sigma}}(\omega )b_{\mu,\sigma}^{\ast}(\omega )\right).\end{eqnarray} \tag{ 92 }$$

To check whether the transformation (88) is a symmetry of the action even in the presence of the driving term in equation (90), it is most convenient to rewrite the transformation in frequency space. Moreover, for future reference, we give the form of the transformation including a chemical potential:

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{\mathcal{T}}_{\beta,\mu}}{{a}_{\sigma}}(\omega ) & ={{\text{e}}^{-\sigma \beta (\omega -\mu )/2}}a_{\sigma}^{\ast}(\omega ), \\ {{\mathcal{T}}_{\beta,\mu}}a_{\sigma}^{\ast}(\omega ) & ={{\text{e}}^{\sigma \beta (\omega -\mu )/2}}{{a}_{\sigma}}(\omega ). \end{array}\end{eqnarray} \tag{ 93 }$$

It is straightforward to check that both S_b and S_sb are invariant under this transformation with $\mu =0$ [189]; the same holds true for the contribution S₀ to the action of the system. However, the driving part (90) becomes after the transformation

$$\begin{eqnarray}&&{{S}_{l}}=\Omega \underset{\sigma}{\sum}\,\sigma \left({{a}_{\sigma}}\left({{\omega}_{l}}\right){{\text{e}}^{\sigma \beta {{\omega}_{l}}/2}}+a_{\sigma}^{\ast}\left({{\omega}_{l}}\right){{\text{e}}^{-\sigma \beta {{\omega}_{l}}/2}}\right),\end{eqnarray} \tag{ 94 }$$

where the appearance of the exponentials shows that the symmetry is violated. One might wonder whether it is possible to restore the symmetry in a rotating frame in which the explicit time dependence of the action is eliminated, i.e. by introducing new variables ${{\tilde{a}}_{\sigma}}(t)$ and ${{\tilde{b}}_{\sigma}}(t)$ rotating at the frequency of the driving field, ${{\tilde{a}}_{\sigma}}(t)={{a}_{\sigma}}(t){{\text{e}}^{\text{i}{{\omega}_{l}}t}}$ and analogously for ${{\tilde{b}}_{\sigma}}(t)$. In terms of these variables, the driving part of the action becomes

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{S}_{l}} & =\Omega \underset{\sigma}{\sum}\,\sigma {\int}^{}\text{d}t\left({{{\tilde{a}}}_{\sigma}}(t)+\tilde{a}_{\sigma}^{\ast}(t)\right) \\ {} & =\Omega \underset{\sigma}{\sum}\,\sigma \left({{{\tilde{a}}}_{\sigma}}(0)+\tilde{a}_{\sigma}^{\ast}(0)\right), \end{array}\end{eqnarray} \tag{ 95 }$$

i.e. the drive couples to the zero-frequency component of ${{\tilde{a}}_{\sigma}}(\omega )$. Clearly, applying again the same transformation equation (93) with $\mu =0$ to the new variables, the driving term, as well as S₀ in equation (89) and the system–bath coupling S_sb in equation (92) are invariant. However, in the action for the bath (91), as a consequence of the transformation to the rotating frame the distribution function acquires an effective chemical potential and becomes $\coth \left(\beta \left(\omega -{{\omega}_{l}}\right)/2\right)$. Hence, to leave this part of the action invariant, the transformation equation (93) has to be applied with $\mu ={{\omega}_{l}}$, and again the full action is not invariant under the equilibrium transformation. No frame exists in which the reference to the driving scale ${{\omega}_{l}}$ were eliminated.

While this simple example demonstrates explicitly that generically external driving takes a system out of thermal equilibrium, and how this is manifest in the symmetry properties of the Keldysh action, it is interesting to note that there are surprising exceptions to this rule, emerging as limiting cases. One of them has been identified in [211]. Along the lines of this reference, we discard the driving term S_l and instead consider a parametric system–bath coupling of the form¹⁹

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{S}_{sb}} & =\lambda \underset{\sigma}{\sum}\,\sigma \underset{\mu}{\sum}\,{\int}^{}\text{d}t\left(a_{\sigma}^{\ast}(t){{b}_{\mu,\sigma}}(t){{\text{e}}^{-\text{i}{{\omega}_{\text{p}}}t}}+{{a}_{\sigma}}(t)b_{\mu,\sigma}^{\ast}(t){{\text{e}}^{\text{i}{{\omega}_{\text{p}}}t}}\right) \\ {} & =\lambda \underset{\mu}{\sum}\,\underset{\sigma}{\sum}\,\sigma {\int}_{\omega}\left(a_{\sigma}^{\ast}\left(\omega +{{\omega}_{\text{p}}}\right){{b}_{\mu,\sigma}}(\omega )+{{a}_{\sigma}}\left(\omega +{{\omega}_{\text{p}}}\right)b_{\mu,\sigma}^{\ast}(\omega )\right), \end{array}\end{eqnarray} \tag{ 96 }$$

in the limit $\lambda \to 0$. Then, the full action $S={{S}_{0}}+{{S}_{sb}}+{{S}_{b}}$, where S₀ and S_b are as above (equations (89) and (91)), is invariant if the system fields are transformed with ${{\mathcal{T}}_{\beta,\mu ={{\omega}_{\text{p}}}}}$ and the bath oscillators with ${{\mathcal{T}}_{\beta,\mu =0}}$. In other words, the parametric coupling acts to thermalize the system at the temperature of the bath while at the same time shifting the chemical potential by ${{\omega}_{\text{p}}}$ with respect to the bath chemical potential. Concomitantly, the FDR for the system variables takes the form of equation (84) with $\omega \to \omega -{{\omega}_{\text{p}}}$, while for the bath it is equation (84) without modification, and only FDRs for cross-correlations between system and bath observables would reveal that there is no true equilibrium in the sense of a single global temperature and chemical potential. Such cross-correlations, however, are suppressed in the limit $\lambda \to 0$. In this limit, both subsystems decouple, and each of them exhibits thermal behavior.

For many applications, as discussed in section 2.3, the Keldysh action in the semiclassical limit (79) is appropriate. Correspondingly we should consider the semiclassical limit of the transformation (88). In the limit of high temperatures $\beta =1/T\to 0$, we can perform an expansion of the transformed fields, with their arguments shifted by $i\sigma \beta /2$, in terms of derivatives. To leading order²⁰, this yields

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{\mathcal{T}}_{\beta}}{{\Phi}_{c}}\left(t,\mathbf{x}\right) & ={{\sigma}_{x}}{{\Phi}_{c}}\left(-t,\mathbf{x}\right), \\ {{\mathcal{T}}_{\beta}}{{\Phi}_{q}}\left(t,\mathbf{x}\right) & ={{\sigma}_{x}}\left({{\Phi}_{q}}\left(-t,\mathbf{x}\right)+\frac{\text{i}}{2T}{{\partial}_{t}}{{\Phi}_{c}}\left(-t,\mathbf{x}\right)\right), \end{array}\end{eqnarray} \tag{ 97 }$$

where ${{\sigma}_{x}}$ is the usual Pauli matrix. 'High temperatures' thus means that the typical frequency scale of the field is much smaller than temperature; indeed this recovers the intuition of a semiclassical limit. If we replace the quantum field ${{\Phi}_{q}}$ by the response field $\rm{\tilde \Phi}=-\text{i}{{\sigma}_{z}}{{\Phi}_{q}}$, equation (97) takes the form of the classical symmetry introduced in [207]. The FDR in the semiclassical limit, which may be derived as a Ward–Takahashi of the symmetry equation (97), simplifies to the Raleigh–Jeans form

$$\begin{eqnarray}&&{{G}^{K}}\left(\omega,\mathbf{q}\right)=\frac{2T}{\omega}\left({{G}^{R}}\left(\omega,\mathbf{q}\right)-{{G}^{A}}\left(\omega,\mathbf{q}\right)\right).\end{eqnarray} \tag{ 98 }$$

The thermal FDR is just one consequence of the presence of the symmetry equation (97). A second one concerns the possible values of the coupling constants defining the action in the semiclassical limit. This allows us to state precisely in which sense the driven–dissipative systems represent a genuine non-equilibrium situation. To this end, it is most convenient to discuss the equivalent Langevin equation (80). We rewrite it by splitting the deterministic parts on the RHS into reversible (coherent) and irreversible (dissipative) contributions according to (here we replace ${{\phi}_{c}}=\psi$)

$$\begin{eqnarray}&&\text{i}{{\partial}_{t}}\psi =\frac{\delta {{H}_{c}}}{\delta {{\psi}^{\ast}}}-\text{i}\frac{\delta {{H}_{d}}}{\delta {{\psi}^{\ast}}}+\xi \end{eqnarray} \tag{ 99 }$$

with effective coherent and dissipative Hamiltonians ($\alpha =c,d$)

$$\begin{eqnarray}&&{{H}_{\alpha}}={{{\int}^{}}_{t,\mathbf{x}}}\left({{K}_{\alpha}}|\nabla \psi {{|}^{2}}\,+\,{{r}_{\alpha}}|\psi {{|}^{2}}\,+\,\frac{{{u}_{\alpha}}}{2}|\psi {{|}^{4}}\right).\end{eqnarray} \tag{ 100 }$$

It can be shown [181, 212] that the presence of the symmetry (97), or, in physical terms, relaxation of a system to thermodynamic equilibrium (a state with global detailed balance, where arbitrary subparts are in equilibrium with each other) requires the condition

$$\begin{eqnarray}&&{{H}_{c}}=r{{H}_{d}}\quad \Leftrightarrow \quad r=\frac{{{K}_{c}}}{{{K}_{d}}}=\frac{{{u}_{c}}}{{{u}_{d}}}.\end{eqnarray} \tag{ 101 }$$

(Note that there is no condition on the ratio ${{r}_{c}}/{{r}_{d}}$ since the effective chemical potential r_c can always be adjusted by a gauge transformation $\psi \mapsto \psi {{\text{e}}^{-\text{i}\omega t}}$ such that ${{r}_{c}}\mapsto {{r}_{c}}-\omega$ without changing the physics.) In equilibrium dynamics, the ratio of real (reversible) and imaginary (dissipative) parts is thus locked to one common value for all couplings. This is illustrated in figure 6(a). In the complex plane spanned by real and imaginary parts of the couplings $K={{K}_{c}}+\text{i}{{K}_{d}}$ and $u={{u}_{c}}+\text{i}{{u}_{d}}$, they lie on one single ray. The intuition behind this seeming fine-tuning is the following: a microscopically reversible dynamics starts from a Hamiltonian functional H_c alone, i.e. all couplings are located on the real axis. Coarse graining the system from the microscopic to the macroscopic scales introduces irreversible dynamics in the form of finite imaginary parts, however preserving their location on a single ray: the ray just rotates under coarse graining, but does not spread out. The geometric constraint is thus not due to fine tuning, but results from the microscopic 'initial condition' for the RG flow, in combination with the presence of a symmetry. In stark contrast, in a driven non-equilibrium system, the microscopic origins of reversible and irreversible dynamics are independent, as illustrated in figure 6(b). For example, in the microscopic description of equation (73) the rates can be tuned fully independently from the Hamiltonian parameters—they have completely different physical origins.

Zoom In Zoom Out Reset image size

Finally, we note that, while an explicit violation of the symmetry is present on the microscopic scales on which the quantum master equation (16) or the corresponding Keldysh action (27) provide a suitable description of the system, in several cases it has been found that driven–dissipative systems appear as approximately thermal at low frequencies [25, 33, 54, 110, 112, 122, 191, 213, 214]. In driven–dissipative condensates in three spatial dimensions, this thermalization at low frequencies is particularly sharply reflected via the emergence of the thermal symmetry in the RG flow in this regime [26, 181], see figure 6(c); this is discussed in more detail in section 4. On the other hand, there are also cases where the opposite behavior occurs, and the non-equilibrium character becomes more pronounced as one coarse-grains to the macroscale. This occurs in driven two dimensional systems, as explained in section 4.

A global continuous symmetry of the Keldysh action is a transformation ${{T}_{\alpha}}$ of the fields $\Psi={{\left({{\psi}_{+}},\psi _{+}^{\ast},{{\psi}_{-}},\psi _{-}^{\ast}\right)}^{T}}$ that leaves the value of the action invariant, i.e.

$$\begin{eqnarray}&&S\left[{{T}_{\alpha}}\Psi\right]=S[\Psi],\end{eqnarray} \tag{ 102 }$$

where α is a real time- and space independent parameter, and for $\alpha =0$ the transformation is the identity, ${{T}_{0}}={\mathbb 1}$. The Noether theorem states that any such global continuous symmetry entails the existence of a current j with components ${{j}^{\mu}},$ which obeys a continuity equation on average, i.e. $\langle {{\partial}_{\mu}}\,{{j}^{\mu}}\rangle =0$ (here and in the following summation over repeated indices is implied), where ${{\partial}_{0}}={{\partial}_{t}}$, and ${{\partial}_{1,2,\ldots d}}$ are derivatives with respect to spatial coordinates. Then, the integral over space $\mathcal{Q}={{{\int}^{}}_{\mathbf{x}}}{{j}^{0}}$—the Noether charge—is an integral of motion, and we have $\langle \text{d}\,\mathcal{Q}/\text{d}t\rangle =0$. In the following, we prove this relation, which states that $\mathcal{Q}$ is conserved on average, in the framework of the Keldysh formalism. We note, however, that a global continuous symmetry implies the even stronger statement $\text{d}\,\mathcal{Q}/\text{d}t=0$ of conservation of $\mathcal{Q}$ on the operator level [86].

In order to prove the Noether theorem, it is sufficient to consider infinitesimal transformations. Then, for $\alpha \ll 1$, we expand the transformation as

$$\begin{eqnarray}&&{{T}_{\alpha}}\Psi=\Psi+\alpha \mathcal{G}\Psi+O\left({{\alpha}^{2}}\right),\end{eqnarray} \tag{ 103 }$$

where $\mathcal{G}$ is called the generator of the transformation. In general, $\mathcal{G}$ is a $4\times 4$ matrix with derivative operators as entries. We consider specific examples below in section 2.4.4. In the following we assume that as in equation (30) the action S can be written in terms of a Lagrangian density $\mathcal{L}$ as $S={{{\int}^{}}_{t,\mathbf{x}}}\mathcal{L}$, and that the Lagrangian density is a local function of the fields and their first derivatives with respect to time and space. This assumption is appropriate for most practical purposes.

In an expansion of the LHS of equation (102) in powers of α, each coefficient has to vanish individually. The first-order contribution yields the relation

$$\begin{eqnarray}&&{{{\int}^{}}_{t,\mathbf{x}}}\left(\frac{\partial \mathcal{L}}{\partial {{\Psi}^{T}}}\mathcal{G}\Psi+\frac{\partial \mathcal{L}}{\partial {{\partial}_{\mu}}{{\Psi}^{T}}}{{\partial}_{\mu}}\mathcal{G}\Psi\right)=0,\end{eqnarray} \tag{ 104 }$$

where $\partial \mathcal{L}/\partial \Psi\equiv {{\left(\partial \mathcal{L}/\partial {{\psi}_{+}},\partial \mathcal{L}/\partial \psi _{+}^{\ast},\partial \mathcal{L}/\partial {{\psi}_{-}},\partial \mathcal{L}/\partial \psi _{-}^{\ast}\right)}^{T}}.$ In the cases we consider, equation (104) holds true because the integrand can be written as the divergence of a vector field ${{f}^{\mu}}$, i.e. we have

$$\begin{eqnarray}&&{{\partial}_{\mu}}\,{{f}^{\mu}}=\frac{\partial \mathcal{L}}{\partial {{\Psi}^{T}}}\mathcal{G}\Psi+\frac{\partial \mathcal{L}}{\partial {{\partial}_{\mu}}{{\Psi}^{T}}}{{\partial}_{\mu}}\mathcal{G}\Psi.\end{eqnarray} \tag{ 105 }$$

To proceed we consider local transformations, i.e. we consider ${{T}_{\alpha}}$ with $\alpha =\alpha \left(t,\mathbf{x}\right)$. We perform a change of integration variables $\Psi\to {{T}_{\alpha}}\Psi$ in the partition function. Then, assuming that the functional measure is invariant with respect to the local transformation, we have

$$\begin{eqnarray}&&Z={\int}^{}\mathcal{D}[\Psi]{{\text{e}}^{\text{i}S[\Psi]}}={\int}^{}\mathcal{D}[\Psi]{{\text{e}}^{\text{i}S\left[{{T}_{\alpha}}\Psi\right]}}.\end{eqnarray} \tag{ 106 }$$

As before, we expand the RHS of this equality in a power series in α. Since by assumption the latter is a function of $\left(t,\mathbf{x}\right)$, equation (102) does not hold true any more. Instead, we find

$$\begin{eqnarray}&&S\left[{{T}_{\alpha}}\Psi\right]=S[\Psi]+{{{\int}^{}}_{t,\mathbf{x}}}\alpha {{\partial}_{\mu}}\left(\,{{f}^{\mu}}-\frac{\partial \mathcal{L}}{\partial {{\partial}_{\mu}}{{\Psi}^{T}}}\mathcal{G}\Psi\right)+O\left({{\alpha}^{2}}\right),\end{eqnarray} \tag{ 107 }$$

where we used equation (105) and integration by parts to write the RHS in a form that does not contain derivatives of α explicitly. Inserting equation (107) in the exponential on the RHS of equation (106), and expanding the latter to first order in α, we obtain the condition

$$\begin{eqnarray}&&{{{\int}^{}}_{t,\mathbf{x}}}\alpha {{\partial}_{\mu}}\langle \,{{f}^{\mu}}-\frac{\partial \mathcal{L}}{\partial {{\partial}_{\mu}}{{\Psi}^{T}}}\mathcal{G}\Psi\rangle =0.\end{eqnarray} \tag{ 108 }$$

Since there are no restrictions on the choice of α, this equation implies that the divergence of the expectation value vanishes. The latter is just the quantum Noether current,

$$\begin{eqnarray}&&{{j}^{\mu}}=\frac{\partial \mathcal{L}}{\partial {{\partial}_{\mu}}{{\Psi}^{T}}}\mathcal{G}\Psi-{{f}^{\mu}}.\end{eqnarray} \tag{ 109 }$$

Note that ${{j}^{\mu}}$ is not unique: for constant $a,{{b}^{\mu}}$, the combination $a{{j}^{\mu}}+{{b}^{\mu}}$ is also a conserved current. The associated Noether charge $\mathcal{Q}$ corresponds to the integral over space of the zeroth component of the Noether current,

$$\begin{eqnarray}&&\mathcal{Q}={{{\int}^{}}_{\mathbf{x}}}{{j}^{0}}={{{\int}^{}}_{\mathbf{x}}}\left(\frac{\partial \mathcal{L}}{\partial {{\partial}_{t}}{{\Psi}^{T}}}\mathcal{G}\Psi-{{f}^{0}}\right).\end{eqnarray} \tag{ 110 }$$

While this derivation is quite general, a more concrete expression for the Noether current can be obtained by restricting the form of the generator $\mathcal{G}$ in equation (103). In particular, in the next section we consider the specific cases of classical and quantum transformations.

Energy and momentum conservation—This property is expected for systems whose classical action is determined by a translationally invariant Hamiltonian alone. Indeed, the construction of the Keldysh action for a Markovian master equation in section 2.1 shows that terms which couple the forward and backward branches are only contained in the dissipative parts of the action (30). In other words, they arise upon coupling the system to a bath and integrating out the latter. On the other hand, the Keldysh action for a closed system with unitary dynamics does not contain such terms and can be written as

$$\begin{eqnarray}&&S=\underset{\sigma}{\sum}\,\sigma {\int}_{t,\mathbf{x}}{{\mathcal{L}}_{\sigma}},\end{eqnarray} \tag{ 111 }$$

where ${{\mathcal{L}}_{+}}$ and ${{\mathcal{L}}_{-}}$ are the Lagrangian densities evaluated with fields on the forward and backward branches respectively. Given this structure of the action, let us consider a transformation ${{T}_{\alpha}}$ which does not mix fields on the forward and backward branches. Then, the generator $\mathcal{G}$ has the block-diagonal structure

$$\begin{eqnarray}\mathcal{G}=\left(\begin{array}{*{35}{l}} {{g}_{+}} & 0 \\ 0 & {{g}_{-}} \end{array}\right).\end{eqnarray} \tag{ 112 }$$

The cases of physical interest are the classical and quantum transformations mentioned above, corresponding to the choices ${{g}_{+}}={{g}_{-}}$ and ${{g}_{+}}=-{{g}_{-}}$ respectively. In both cases, the Noether current can be represented as superposition of currents on the forward and backward branches, which we define in terms of $g\equiv {{g}_{+}}$ as

$$\begin{eqnarray}&&j_{\sigma}^{\mu}=\frac{\partial {{\mathcal{L}}_{\sigma}}}{\partial {{\partial}_{\mu}}\Psi_{\sigma}^{T}}g{{\Psi}_{\sigma}}-f_{\sigma}^{\mu},\end{eqnarray} \tag{ 113 }$$

where

$$\begin{eqnarray}&&{{\partial}_{\mu}}f_{\sigma}^{\mu}=\frac{\partial {{\mathcal{L}}_{\sigma}}}{\partial \Psi_{\sigma}^{T}}g{{\Psi}_{\sigma}}+\frac{\partial {{\mathcal{L}}_{\sigma}}}{\partial {{\partial}_{\mu}}\Psi_{\sigma}^{T}}{{\partial}_{\mu}}g{{\Psi}_{\sigma}}.\end{eqnarray} \tag{ 114 }$$

With this definition, j₋ can be obtained from j₊ simply by replacing all instances of ${{\Psi}_{+}}$ appearing in j₊ by ${{\Psi}_{-}}$. Then, the symmetry of the action under a quantum transformation yields a classical Noether current and vice versa. These currents are given by (as noted above we are free to choose a convenient multiplicative normalization of the currents)

$$\begin{eqnarray}&&{{j}_{c}}=\frac{1}{2}\left(\,{{j}_{+}}+{{j}_{-}}\right),\quad {{j}_{q}}={{j}_{+}}-{{j}_{-}}.\end{eqnarray} \tag{ 115 }$$

As pointed out in section 2.2.2, due to causality only the classical component of the field can acquire a finite expectation value. The same is true for the currents defined in equation (115): in the presence of an external field that induces a current, we have $\langle \,{{j}_{+}}\rangle =\langle \,{{j}_{-}}\rangle =\langle \,{{j}_{c}}\rangle \ne 0,$ whereas $\langle \,{{j}_{q}}\rangle =0.$ Hence, as anticipated above, the continuity equation $\langle {{\partial}_{\mu}}\,j_{q}^{\mu}\rangle =0$, which follows from the symmetry under a classical transformation, does not entail a non-trivial Noether charge. In the next section, we show that such a symmetry nevertheless does have physical consequences, as it can be broken spontaneously in a condensation transition.

First, let us derive the Noether charge associated with the symmetry of a Hamiltonian Keldysh action under contour- dependent—i.e. quantum—translations in space-time. To be specific, we consider the transformation ${{\Psi}_{\sigma}}(X)\mapsto {{\Psi}_{\sigma}}\left(X+\sigma \alpha {{e}_{\nu}}\right)$, where $X=\left(t,\mathbf{x}\right)$, and ${{e}_{\nu}}$ is the basis vector in the direction of the νth coordinate of d  +  1-dimensional space-time. This is a symmetry of the action (111): the shift of the coordinates can be absorbed by performing a change of integration variables $X\to X-\sigma \alpha {{e}_{\nu}}$ in the integral over the Lagrangian density ${{\mathcal{L}}_{\sigma}}$. Clearly, the symmetry is violated in the presence of dissipative terms that couple the forward and backward branches.

For space-time translations, the generators ${{g}_{\pm }}$ in equation (112) acquire an additional index ν and are given by ${{g}_{\pm,\nu}}=\pm {{\partial}_{\nu}}$. It follows that the vector field $f_{\sigma,\nu}^{\mu}$ in equation (114) takes the form $f_{\sigma,\nu}^{\mu}=\delta _{\nu}^{\mu}{{\mathcal{L}}_{\sigma}}$, and the classical Noether current, which we denote as ${{T}_{c,\nu}}$, is given by

$$\begin{eqnarray}&&T_{c,\nu}^{\mu}=\frac{1}{2}\underset{\sigma}{\sum}\,T_{\sigma,\nu}^{\mu},\quad T_{\sigma,\nu}^{\mu}=\frac{\partial {{\mathcal{L}}_{\sigma}}}{\partial {{\partial}_{\mu}}\Psi_{\sigma}^{T}}{{\partial}_{\nu}}{{\Psi}_{\sigma}}-\delta _{\nu}^{\mu}{{\mathcal{L}}_{\sigma}},\end{eqnarray} \tag{ 116 }$$

where $T_{\pm,\nu}^{\mu}$ are the components of the energy-momentum tensor [86, 205] on the forward and backward branches. In particular, the component $T_{\pm,0}^{0}$ is the energy density, whereas $T_{\pm,i}^{0}$ is the density of momentum in the spatial direction $i=1,2,\ldots,d$. The spatial integrals over these densities yield the associated Noether charges, e.g. if we express the Lagrangian density ${{\mathcal{L}}_{\sigma}}$ in terms of a Hamiltonian density (as in equation (30), where $\mathcal{L}$ is given by equation (28) with ${{\gamma}_{\alpha}}=0$) as

$$\begin{eqnarray}&&{{\mathcal{L}}_{\sigma}}=\frac{1}{2}\Psi_{\sigma}^{\dagger}\text{i}{{\sigma}_{z}}{{\partial}_{t}}{{\Psi}_{\sigma}}-{{H}_{\sigma}},\end{eqnarray} \tag{ 117 }$$

we obtain the conserved energy density $\mathcal{E}$, which is as expected given by the sum of the Hamiltonian densities on the forward and backward branches,

$$\begin{eqnarray}&&\mathcal{E}=\frac{1}{2}\underset{\sigma}{\sum}\,{{{\int}^{}}_{\mathbf{x}}}T_{\sigma,0}^{0}=\frac{1}{2}{\int}_{\mathbf{x}}\left({{H}_{+}}+{{H}_{-}}\right).\end{eqnarray} \tag{ 118 }$$

Along the same lines, conservation of angular momentum follows from the quantum symmetry of the Keldysh action with respect to rotations of the spatial coordinates.

Number conservation—As another example, we consider phase rotations and their relation to particle number conservation. In this case the classical transformation reads ${{\psi}_{\pm }}\mapsto {{U}_{c}}{{\psi}_{\pm }}={{\text{e}}^{\text{i}\alpha}}{{\psi}_{\pm }}$ and the quantum transformation is ${{\psi}_{\pm }}\mapsto {{U}_{q}}{{\psi}_{\pm }}={{\text{e}}^{\pm \text{i}\alpha}}{{\psi}_{\pm }}$. Classical phase rotations are a symmetry of the Keldysh action, if the fields appear only in the combinations $\psi _{\sigma}^{\ast}{{\psi}_{{{\sigma}^{\prime}}}}$. The quantum transformation is more restrictive: in this case only products with $\sigma ={{\sigma}^{\prime}}$ are allowed. Hence, the symmetry with regard to quantum phase rotations implies the classical symmetry. The generators of U_q are given by ${{g}_{\pm }}=\pm \text{i}{{\sigma}_{z}}$. Then, considering particles with quadratic dispersion relation and inserting the representation of the Lagrangian density in equation (117) in the expression for the Noether current (113), we find that the latter can be written as ${{j}_{\sigma}}={{\left({{\rho}_{\sigma}},{{\mathbf{j}}_{\sigma}}\right)}^{T}},$ where ${{\rho}_{\sigma}}=|{{\psi}_{\sigma}}{{|}^{2}}$ is the density on the contour with index σ and the spatial components of the current are given by

$$\begin{eqnarray}&&{{\mathbf{j}}_{\sigma}}=\frac{1}{i2m}\left(\psi _{\sigma}^{\ast}\nabla {{\psi}_{\sigma}}-{{\psi}_{\sigma}}\nabla \psi _{\sigma}^{\ast}\right),\end{eqnarray} \tag{ 119 }$$

which is just the ordinary quantum mechanical current, evaluated on the closed time path. What is the physical meaning of the classical and quantum components, which are defined in equation (115), of the current given in equation (119)? If we introduce in the Hamiltonian an externally imposed gauge field, in the Keldysh action it would couple to the quantum current, by analogy to the source terms in the generating functional equation (33). However, the observable effect of such a gauge field is that it can induce a classical current $\langle \,{{j}_{c}}\rangle \ne 0$.

We close this section with a comment on the status of the U_q(1) symmetry. It will be present on the microscopic level for an action corresponding to a Hamiltonian which commutes with the number operator. Since the symmetry considerations are properly applied to the microscopic action (and the functional measure), the conservation law ensues. However, this does not imply that the effective action must manifestly preserve a U_q(1) symmetry. In fact, classical models of number conserving dynamics [1] do not exhibit this symmetry. This leads to the picture that this symmetry is broken spontaneously under RG. The precise workings of such a mechanism are, to the best of our knowledge, not settled so far.

The conservation of particle number in a closed system follows from the continuity equation for the classical Noether current that is associated with the symmetry under quantum phase rotations. In an open system, this symmetry is absent, and the continuity equation has to be extended to account for the exchange of particles with the environment, as we discuss in the following.

To be specific, we consider the Keldysh action (73) for a system with single-particle pump as well as single-particle and two-body loss. In order to derive the extended continuity equation, as in equation (106) we perform a change of integration variables $\Psi\to {{U}_{q}}\Psi$, where U_q is a local quantum phase rotation, in the Keldysh partition function. Since the partition function is invariant under this transformation, the coefficients in an expansion of the partition function in powers of the phase shift must vanish. In linear order, we find the condition

$$\begin{eqnarray}&&\langle {{\partial}_{\mu}}\,j_{c}^{\mu}\rangle -{{\gamma}_{p}}\langle \psi _{+}^{\ast}{{\psi}_{-}}\rangle +{{\gamma}_{l}}\langle \psi _{-}^{\ast}{{\psi}_{+}}\rangle +4{{u}_{d}}\langle {{\left(\psi _{-}^{\ast}{{\psi}_{+}}\right)}^{2}}\rangle =0.\end{eqnarray} \tag{ 120 }$$

This should be compared to equation (108): if U_q were a symmetry of the action, we would have found an ordinary continuity equation. This would have been the case in the absence of the pumping and loss terms proportional to ${{\gamma}_{p}},{{\gamma}_{l}},$ and u_d. The interpretation of these terms is as follows: in a system that is perfectly isolated from its environment, the temporal change of the density at a given point is due to the motion of particles toward or away from this point, i.e. the rate of change of density is a source or sink for the mass flow as measured by the divergence of the current $\langle \nabla \centerdot {{\mathbf{j}}_{c}}\rangle$. In an open system, on the other hand, there is in addition a dissipative current [215, 216] due to the exchange of particles between the system and its environment. In particular, in equation (120) this dissipative current has contributions from the pumping and loss terms.

It is interesting to note that equation (120) can also be obtained from the equation of motion of the local density $n\left(\mathbf{x}\right)={{\psi}^{\dagger}}\left(\mathbf{x}\right)\psi \left(\mathbf{x}\right)$ in the operator formalism. From the master equation ${{\partial}_{t}}\rho =\mathcal{L}\rho$ with Liouvillian $\mathcal{L}$ in Lindblad form given in equation (9), it follows that the time evolution of $n\left(t,\mathbf{x}\right)$ in the Heisenberg representation is given by ${{\partial}_{t}}n\left(t,\mathbf{x}\right)={{\mathcal{L}}^{\ast}}n\left(t,\mathbf{x}\right)$ with the adjoint Liouvillian defined by $\text{tr}\left(A\mathcal{L}B\right)=\text{tr}\left(B{{\mathcal{L}}^{\ast}}A\right)$. We find

$$\begin{eqnarray}&&{{\partial}_{t}}n\left(t,\mathbf{x}\right)=-\nabla \centerdot \mathbf{j}\left(t,\mathbf{x}\right)+{{\gamma}_{p}}\psi \left(t,\mathbf{x}\right){{\psi}^{\dagger}}\left(t,\mathbf{x}\right)\\ &&-{{\gamma}_{l}}n\left(t,\mathbf{x}\right)-4{{u}_{d}}{{\psi}^{\dagger}}{{\left(t,\mathbf{x}\right)}^{2}}\psi {{\left(t,\mathbf{x}\right)}^{2}},\end{eqnarray} \tag{ 121 }$$

where the first term encodes coherent dynamics and corresponds to the Heisenberg commutator $i\left[{{H}_{\text{LP}}},n\left(t,\mathbf{x}\right)\right]$, whereas the remaining contributions incorporate the dissipative parts. Taking the average of this relation and expressing the expectation values as Keldysh functional integrals yields equation (120).

Finally, we can obtain some intuition on the dissipative correction to the current expectation value in mean field theory, where the correlators in equation (120) factorize into products of single field expectation values $\langle {{\psi}_{+}}\rangle =\langle {{\psi}_{-}}\rangle ={{\psi}_{0}}$. We then recognize in the non-derivative terms ($\psi _{0}^{\ast}$ times) the LHS of equation (75) (note that due to the factor of $1/\sqrt{2}$ in the Keldysh rotation (32) the expectation values are related as $\langle {{\phi}_{c}}\rangle ={{\phi}_{0}}=\sqrt{2}\langle {{\psi}_{\pm }}\rangle$, and that ${{r}_{d}}=\left({{\gamma}_{l}}-{{\gamma}_{p}}\right)/2$), which equals zero in a homogeneous situation within mean field theory. In this way, we see that there is no particle number current on average in a homogeneous driven–dissipative system in stationary state, as expected.

Even though the pumping and loss terms in the Keldysh action in equation (73) break the symmetry with respect to quantum phase rotations U_q, the action is still symmetric under U_c transformations. As a result, even in the absence of particle number conservation there is a possibility of spontaneous symmetry breaking. In particular, a finite average value $\langle {{\phi}_{c}}\rangle \ne 0$ breaks the classical phase rotation symmetry in a non-equilibrium condensation transition. Here we show that this spontaneous symmetry breaking is accompanied by the appearance of a massless mode, i.e. a mode with vanishing frequency and decay rate for zero momentum, which is known as the Goldstone boson [217–219]. We obtain this result by carrying over the corresponding derivation from the equilibrium formalism (see, e.g. [182]) to the Keldysh framework.

The single-particle excitation spectrum is encoded in the poles of retarded and advanced Green's functions, and such a pole at $\omega =q=0$ is dubbed a massless mode. Poles of the Green's functions correspond to roots of the determinant of inverse propagator (see section 2.1). Hence, the presence of a massless mode can be detected by checking whether the determinant of the mass matrix M vanishes. The latter is determined by the inverse propagator

$$\begin{eqnarray}P\left(\omega,\mathbf{q}\right)=\left(\begin{array}{*{35}{c}} 0 & {{P}^{A}}\left(\omega,\mathbf{q}\right) \\ {{P}^{R}}\left(\omega,\mathbf{q}\right) & {{P}^{K}}\left(\omega,\mathbf{q}\right) \end{array}\right)={{G}^{-1}}\left(\omega,\mathbf{q}\right),\end{eqnarray} \tag{ 122 }$$

given exemplarily in equation (76), at $\omega =\mathbf{q}=0$, i.e. M  =  − P(0, 0). Therefore, it is sufficient to consider frequency- and momentum-independent field configurations or, in other words, homogeneous field configurations, so that instead of the full effective action equation (39) we only have to deal with the effective potential $U=-\Gamma{{|}_{\text{hom}.}}/\Omega$, where $\Omega$ is the quantization volume in space-time. Since the effective potential U inherits the symmetries of the effective action, it is a function of precisely these combinations of fields which are invariant under classical or quantum phase rotations. In the basis of classical and quantum fields, these U_c-symmetric combinations are ${{\rho}_{\nu {{\nu}^{\prime}}}}=\phi _{\nu}^{\ast}{{\phi}_{{{\nu}^{\prime}}}}$, where the indices ν and ${{\nu}^{\prime}}$ can take the values c and q. In the following, we show that U_c invariance is sufficient to guarantee the existence of a massless mode.

For convenience we switch to a basis of real fields $\chi =\left({{\chi}_{c,1}},{{\chi}_{c,2}},{{\chi}_{q,1}},{{\chi}_{q,2}}\right)$ which correspond to the real and imaginary parts of the complex classical and quantum fields, i.e. ${{\phi}_{\nu}}=\frac{1}{\sqrt{2}}\left({{\chi}_{\nu,1}}+\text{i}{{\chi}_{\nu,2}}\right)$ for $\nu =c,q$. In this basis, the mass matrix reads

$$\begin{eqnarray}&&{{M}_{ij}}=\frac{{{\partial}^{2}}U}{\partial {{\chi}_{i}}\partial {{\chi}_{j}}}{{|}_{\text{ss}}}=\underset{a}{\sum}\,\,{{\left[\frac{{{\partial}^{2}}{{\rho}_{a}}}{\partial {{\chi}_{i}}\partial {{\chi}_{j}}}\frac{\partial U}{\partial {{\rho}_{a}}}\right]}_{\text{ss}}}+\underset{a,b}{\sum}\,\,{{\left[\frac{\partial {{\rho}_{a}}}{\partial {{\chi}_{i}}}\frac{\partial {{\rho}_{b}}}{\partial {{\chi}_{i}}}\frac{{{\partial}^{2}}U}{\partial {{\rho}_{a}}\partial {{\rho}_{b}}}\right]}_{\text{ss}}},\end{eqnarray} \tag{ 123 }$$

where the subscript $\text{ss}$ indicates that the fields should be set to their average values in the stationary state. The indices i and j label the components of the four-vector χ defined above (i.e. ${{\chi}_{1}}={{\chi}_{c,1}}$, ${{\chi}_{2}}={{\chi}_{c,2}}$ etc.), and a and b are double indices, taking the values cc, cq, qc, qq. Let us consider the first term on the RHS of equation (123): in the ordered phase, the classical field has a finite expectation value in the stationary state. Without loss of generality we assume that this value is real. Then, the field equations

$$\begin{eqnarray}&&\frac{\partial U}{\partial {{\chi}_{i}}}{{|}_{\text{ss}}}=\underset{a}{\sum}\,\,{{\left[\frac{\partial {{\rho}_{a}}}{\partial {{\chi}_{i}}}\frac{\partial U}{\partial {{\rho}_{a}}}\right]}_{\text{ss}}}=0,\end{eqnarray} \tag{ 124 }$$

which actually determine the average value of the fields ${{\chi}_{i}}$ in stationary state, have the solution ${{\chi}_{i}}{{|}_{\text{ss}}}=\sqrt{2{{\rho}_{0}}}{{\delta}_{i,1}}$. Performing the derivatives $\partial {{\rho}_{a}}/\partial {{\chi}_{i}}$ in equation (124) explicitly and inserting ${{\chi}_{i}}{{|}_{\text{ss}}}$, we obtain the following conditions:

$$\begin{eqnarray}&&\frac{\partial U}{\partial {{\rho}_{cc}}}{{|}_{\text{ss}}}=\frac{\partial U}{\partial {{\rho}_{cq}}}{{|}_{\text{ss}}}=\frac{\partial U}{\partial {{\rho}_{qc}}}{{|}_{\text{ss}}}=0.\end{eqnarray} \tag{ 125 }$$

Therefore, only the derivative ${{u}_{qq}}={{\left[\partial U/\partial {{\rho}_{qq}}\right]}_{\text{ss}}}$ contributes to the first term on the RHS of equation (123). Denoting mixed second derivatives of the effective potential as ${{u}_{ab}}={{u}_{ba}}={{\left[{{\partial}^{2}}U/\partial {{\rho}_{a}}\partial {{\rho}_{b}}\right]}_{\text{ss}}}$, we find that the mass matrix can be written as

$$\begin{eqnarray}M=\left(\begin{array}{c} 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & {{u}_{qq}} & 0 \\ 0 & 0 & 0 & {{u}_{qq}} \end{array}\right)+{{\rho}_{0}}\left(\begin{array}{c} 0 & 0 & {{u}_{cc,cq}}+{{u}_{cc,qc}} & \text{i}\left({{u}_{cc,cq}}-{{u}_{cc,qc}}\right) \\ 0 & 0 & 0 & 0 \\ {{u}_{cc,cq}}+{{u}_{cc,qc}} & 0 & \frac{1}{2}\left({{u}_{cq,cq}}+2{{u}_{cq,qc}}+{{u}_{qc,qc}}\right) & \frac{\text{i}}{2}\left({{u}_{cq,cq}}-{{u}_{qc,qc}}\right) \\ \text{i}\left({{u}_{cc,cq}}-{{u}_{cc,qc}}\right) & 0 & \frac{\text{i}}{2}\left({{u}_{cq,cq}}-{{u}_{qc,qc}}\right) & -\frac{1}{2}\left({{u}_{cq,cq}}-2{{u}_{cq,qc}}+{{u}_{qc,qc}}\right) \end{array}\right).\end{eqnarray} \tag{ 126 }$$

An entry u_{cc, cc} is not ruled out by symmetry; however, it must vanish due to conservation of probability, see section 2.2.1. Crucially, the retarded and advanced sectors feature one zero eigenvalue. This proves the existence of a massless mode as a consequence of spontaneous symmetry breaking for a theory with U_c invariance.

While this analysis guarantees the existence of a zero mode (i.e. the complex dispersion relation has the property $\omega (q=0)=0$), it does not provide a statement on the functional dependence $\omega (q)$. The latter could be inferred along the lines of [220–224]. As we found in section 2.2.2 in Bogoliubov approximation, in an open system without (i.e. broken) U_q symmetry and spontaneously broken U_c symmetry, the leading behavior at low momenta is diffusive, $\omega (q)\sim -\text{i}D{{q}^{2}}$, with a real diffusion coefficient D. This has to be contrasted to closed systems, in which microscopically both U_q and U_c symmetry are present. There, the leading behavior in a phase with spontaneously broken U_c symmetry is that of coherent sound waves, $\omega \left(\mathbf{q}\right)\sim cq$, with real speed of sound c. A phenomenological justification of this behavior is given in [1], where this can be understood as a consequence of a coupling to additional slow modes relating to particle number conservation. In these phenomenological models, U_q symmetry is absent. This suggests a scenario of an additional spontaneous breakdown of U_q symmetry upon coarse graining, but this issue seems not to be settled to date.

The single-mode Dicke Hamiltonian (133) has the property that all atomic variables couple to the cavity photon mode via the same coupling constant g. Consequently, the coupling of the photons to the sum of the individual Pauli matrices can be replaced by the coupling of the cavity photon mode to a large spin,

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} H={{\omega}_{\text{p}}}\,{{a}^{\dagger}}a+{{\omega}_{z}}{{S}^{z}}+\frac{2g}{\sqrt{N}}{{S}^{x}}\left({{a}^{\dagger}}+a\right), \end{array}\end{eqnarray} \tag{ 136 }$$

where ${{S}^{z}},{{S}^{x}}$ are spin operators in a spin-N/2 Hilbert space. In order to find a path integral representation for the Hamiltonian (136), these spin operators can be transformed into bosonic operators via the common Schwinger-boson or Holstein–Primakoff transformations. For weak coupling g, the system remains strongly polarized $|\langle {{S}_{z}}\rangle |\gg 1$ and the Holstein–Primakoff transformation around the non-interacting ground state, which has the eigenvalue S_z  =  −N/2, is the most natural choice. It reads

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{S}^{z}}={{b}^{\dagger}}b-\frac{N}{2},\,\,\,{{S}^{x}}=\left(\sqrt{N-{{b}^{\dagger}}b}\right)b+{{b}^{\dagger}}\left(\sqrt{N-{{b}^{\dagger}}b}\right), \end{array}\end{eqnarray} \tag{ 137 }$$

where $b,{{b}^{\dagger}}$ are bosonic operators. In the thermodynamic limit $N\to \infty$, the square root in the S^x operator is expanded in powers of 1/N, and the Hamiltonian is subsequently formulated on the Keldysh contour. To zeroth order in 1/N, the corresponding Keldysh action is

$$\begin{eqnarray}S={{{\int}^{}}_{\omega}}{{\Phi}^{\dagger}}(\omega )\,\left(\begin{array}{c} 0 & {{\left(G_{4\times 4}^{A}\right)}^{-1}} \\ {{\left(G_{4\times 4}^{R}\right)}^{-1}} & \Sigma_{4\times 4}^{K} \end{array}\right)\,\Phi(\omega ),\end{eqnarray} \tag{ 138 }$$

with the combined Nambu–Keldysh spinor

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \Phi(\omega )={{\left({{a}_{c}}(\omega ),a_{c}^{\ast}(-\omega ),{{b}_{c}}(\omega ),b_{c}^{\ast}(-\omega ),{{a}_{q}}(\omega ),a_{q}^{\ast}(-\omega ),{{b}_{q}}(\omega ),b_{q}^{\ast}(-\omega )\right)}^{T}}, \end{array}\end{eqnarray} \tag{ 139 }$$

the inverse retarded Green's function in Nambu space

$$\begin{eqnarray}{{\left(G_{4\times 4}^{R}\right)}^{-1}}=\left(\begin{array}{c} \omega -{{\omega}_{\text{p}}}+\text{i}\kappa & 0 & -g & -g \\ 0 & -\omega -{{\omega}_{\text{p}}}-\text{i}\kappa & -g & -g \\ -g & -g & \omega -{{\omega}_{z}} & 0 \\ -g & -g & 0 & -\omega -{{\omega}_{z}} \end{array}\right)\end{eqnarray} \tag{ 140 }$$

and the Keldysh self-energy

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{\Sigma}^{K}}=2\text{i}\kappa \,\text{diag}(1,1,0,0). \end{array}\end{eqnarray} \tag{ 141 }$$

Due to the photon decay terms $\sim \kappa$ and the atom–photon coupling g, the above theory is regularized even without infinitesimal imaginary contributions in the atomic sector. When integrating out the photons, the latter infinitesimal contributions are overwritten in any case by the finite imaginary part of the photon Green's function, and it is therefore reasonable to leave them out from the start.

The excitation spectrum of the atom–photon system is encoded in the retarded Green's function, and the poles, marking the excitation energies, fulfill the requirement

$$\begin{eqnarray}&&0\overset{{}}{{\overset{!}{{=}}\,}}\,\det {{\left(G_{4\times 4}^{R}\right)}^{-1}}=\left({{\omega}^{2}}-\omega _{z}^{2}\right)\left[{{\left(\omega +\text{i}\kappa \right)}^{2}}-\omega _{\text{p}}^{2}\right]-4{{g}^{2}}{{\omega}_{\text{p}}}{{\omega}_{z}}.\end{eqnarray} \tag{ 142 }$$

In the absence of atom–photon coupling, i.e. for g  =  0, these are the non-interacting poles ${{\omega}_{1,2}}=\pm {{\omega}_{z}}$, corresponding to the atomic transition frequencies and ${{\omega}_{3,4}}=\pm {{\omega}_{\text{p}}}+\text{i}\kappa$, corresponding to the energy of a photon ${{\omega}_{\text{p}}}$ and its decay rate κ. For g  >  0, the modes begin to hybridize and the excitation energies are slightly modified compared to their non-interacting values, see figure 7. For small coupling strength, the elementary excitations can still be seen as weakly dressed atoms and photons, with excitation energies close to the non-interacting values, while they are strongly hybridized, inseparable degrees of freedom for strong coupling strengths.

Zoom In Zoom Out Reset image size

Above a critical coupling strength g_c, the ground state of the system breaks the Z₂ symmetry and the atoms form a macroscopic spin aligned in the x-direction, $\langle {{S}_{x}}\rangle =O(N)$, which is accompanied with a coherent macroscopic population of the intra-cavity mode $\langle a\rangle \ne 0$. This symmetry-broken phase is called the superradiant phase of the cavity system. In the case of a macroscopic expectation of S_x, the orthogonal S_z component can no longer be macroscopically large, which renders any expansion of the Holstein–Primakoff operators (137) in 1/N invalid. This is expressed by an instability of the quadratic theory at the superradiance transition, revealed by the presence of a zero energy mode, i.e. a pole at $\omega =0$ in the excitation spectrum. According to equation (142), this happens at

$$\begin{eqnarray}&&{{g}_{c}}=\sqrt{\frac{\left({{\kappa}^{2}}+\omega _{\text{p}}^{2}\right){{\omega}_{z}}}{4{{\omega}_{\text{p}}}}}.\end{eqnarray} \tag{ 143 }$$

The mode structure in the vicinity of the transition has been analyzed in [112], where it has been found that in the presence of photon decay, the critical mode becomes purely imaginary before the transition happens, and the corresponding critical dynamics corresponds to a classical finite temperature transition, see also figure 7. This is mirrored by the fact that the photons effectively thermalize in the low energy regime and their correlations can be described by a low energy effective temperature ${{T}_{\text{eff}}}$. The effective temperature can be obtained via a fluctuation dissipation relation, as discussed in equation (72) but promoted to Nambu space (see section 5.2 for a discussion of the FDR in Nambu space). Solving this equation yields the photon distribution function in Nambu representation

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{F}_{\text{ph}}}={{\sigma}^{z}}+{{\sigma}^{x}}\frac{2{{g}^{2}}}{{{\omega}_{z}}\omega}\frac{\omega _{z}^{2}}{{{\omega}^{2}}-\omega _{z}^{2}}. \end{array}\end{eqnarray} \tag{ 144 }$$

For large frequencies $\omega \gg \frac{{{g}^{2}}}{{{\omega}_{z}}}$, this corresponds to the zero temperature distribution of the non-interacting system $F={{\sigma}^{z}}$, which is diagonal in the photon modes. The approach to this limit is, however, algebraic $\sim {{\omega}^{-3}}$, in contrast to an exponential approach according to large frequency behavior of the Bose-distribution function in equilibrium. On the other hand, for small frequencies $\omega \ll \frac{{{g}^{2}}}{{{\omega}_{z}}}$, the occupation becomes essentially thermal and purely off-diagonal $F=\frac{2{{T}_{\text{eff}}}}{\omega}{{\sigma}^{x}}$, with low energy effective temperature ${{T}_{\text{eff}}}=\frac{{{g}^{2}}}{{{\omega}_{z}}}$, see equation (84). In this low frequency limit, the elementary excitations are strongly hybridized polariton modes, which are diagonal in the photon quadratures x and p. The absence of a global temperature scale in the present system reflects the fact that due to the Markovian dissipation, this interacting system with a discrete set of degrees of freedom is not able to achieve detailed balance between the particular modes, i.e. the x and p quadratures, for all energy scales.

One should note that ${{T}_{\text{eff}}}$ is not proportional to the decay rate κ and consequently the zero temperature equilibrium limit is not simply obtained by taking the limit $\kappa \to 0$. The reason is that in the present setting, the atom–photon coupling g leads to a competition of unitary and dissipative dynamics, which is reflected in the fact that for g  =  0 the eigenstate of the Hamiltonian is a steady state of the dynamics while this is no longer the case for finite g. As a consequence, the low energy effective temperature depends on g and vanishes in the limit $g\to 0$.

A further consequence of the dissipative nature of the transition is that the critical exponents correspond to the classical finite temperature equilibrium transition. One such critical exponent of the superradiant transition is the so-called photon flux exponent ${{\nu}_{x}}$, which describes the divergence of the cavity photon number at the transition, i.e. $n\sim |g-{{g}_{c}}{{|}^{-{{\nu}_{x}}}}$. It is obtained by frequency integration over the photonic correlation function

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} 2n=\left({{{\int}^{}}_{\omega}}\text{i}G_{\text{ph}}^{K}(\omega )\right)-1\sim |g-{{g}_{c}}{{|}^{-1}}. \end{array}\end{eqnarray} \tag{ 145 }$$

The exponent ${{\nu}_{x}}=1$ found for the present non-equilibrium transition coincides with the classical, finite temperature exponent for the equilibrium Dicke model in line with the discussion above. On the other hand, at zero temperature, the exponent of the corresponding quantum phase transition is found to be ${{\nu}_{x}}=1/2$.

The scaling behavior of the critical mode at the transition is expressed in terms of the dynamical critical exponent ${{\nu}_{t}}$. In the presence of dissipation, the excitations in the system decay exponentially in time, which results in an asymptotic exponential decay of real-time correlation functions

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \langle \left\{a(t),{{a}^{\dagger}}(0)\right\}\rangle =G_{\text{ph}}^{K}(t)\sim {{\text{e}}^{-t/{{\xi}_{t}}}}. \end{array}\end{eqnarray} \tag{ 146 }$$

The characteristic decay time ${{\xi}_{t}}$ is determined by the slowest decaying mode in the system, i.e. by the mode with the smallest imaginary part in the frequency spectrum. In the vicinity of the phase transition, this is the critical mode, which shows scaling

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{\xi}_{t}}\sim \frac{1}{|g-{{g}_{c}}|}, \end{array}\end{eqnarray} \tag{ 147 }$$

i.e. a dynamical critical exponent ${{\nu}_{t}}=1$. This behavior is in stark contrast to a zero temperature quantum phase transition, where correlation functions do not decay over time but oscillate with characteristic frequencies ${{\omega}_{c}}$, the smallest of which indicates the distance to the phase transition and encodes the critical scaling behavior. The fact that in the present setting the critical mode becomes purely imaginary is a generic feature of dissipative phase transitions, which is discussed further in section 4.

All of the above-discussed results are encoded in the quadratic Keldysh action, described by equations (138)–(141). The resulting critical behavior is then that of a non-interacting system, and is associated accordingly to a Gaussian fixed point of the renormalization group flow. This is in contrast to critical behavior in interacting systems, which are described by an interacting Wilson–Fisher fixed point, not smoothly connected to the previous one (see the discussion in section 4.2). This action describes the problem up to 1/N corrections resulting from the expansion of the Holstein–Primakoff bosons in 1/N. It contains correctly the physics in the thermodynamic limit $N\to \infty$ and captures the essential dynamics of the superradiance transition. For finite systems, the higher order terms in 1/N can not be neglected and have to be taken into account properly. Analytical approaches, relying on an expansion of the Holstein–Primakoff operators up to first order in the ratio 1/N, have shown that the first order correction term leads to a quartic contribution in the action, whose self-consistently determined one-loop corrections are already in very good agreement with numerical results [112]. The latter have been obtained from a Monte-Carlo wave function (MCWF) simulation of the master equation and agreed with the analytical results even on the level of N  =  10 atoms.

The Dicke model defined by the master equation (134) obeys the common Ising Z₂ symmetry, which is spontaneously broken by the ground state of the superradiant phase. Another approach to express the Dicke model in a Keldysh path integral approach is therefore to represent the spin operators in terms of real Ising field variables, ${{\sigma}_{x}}(t)\to \phi (t)$. In contrast to the Holstein–Primakoff transformation applied in the previous section, the Ising representation does not require a single, large spin to have a well defined field theory representation, and is therefore also applicable in situations in which the atom–photon coupling is spatially dependent $g\to g(x)$ and a large spin representation becomes impossible. Furthermore, we use this approach to describe the symmetry-broken phase of the problem.

The atomic sector of the Hamiltonian (133) describes Ising spins in a transverse field and the corresponding universal low energy description for frequencies below the level spacing ${{\omega}_{z}}$ is obtained by transforming the quantum spin operators to classical real fields. According to the quantum to classical mapping [255] this is possible in the scaling regime of the corresponding, i.e. for the present case in the vicinity of the 2nd order superradiance and glass transitions. The present model represents a generalization of the zero-dimensional Ising model in a transverse field, for which the quantum to classical mapping is known [255] and consists of the replacements

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \sigma _{i}^{x}(t)\to {{\phi}_{i}}(t) \end{array}\end{eqnarray} \tag{ 148 }$$

and

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \sigma _{i}^{z}(t)\to 1-\frac{2}{\omega _{z}^{2}}{{\left({{\partial}_{t}}{{\phi}_{i}}(t)\right)}^{2}}. \end{array}\end{eqnarray} \tag{ 149 }$$

The dynamic constraint ${{\left(\sigma _{i}^{x}(t)\right)}^{2}}=1$ is implemented via the non-linear constraint

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \delta \left(\phi _{i}^{2}(t)-1\right)={\int}^{}\mathcal{D}{{\lambda}_{l}}(t)~{{\text{e}}^{\text{i}{{\lambda}_{l}}(t)\left(\phi _{i}^{2}(t)-1\right)}}, \end{array}\end{eqnarray} \tag{ 150 }$$

which introduces dynamical Lagrange parameters ${{\lambda}_{i}}(t)$ for each spin variable. In this framework, the purely atomic part of the action on the $(\pm )$-contour is

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{S}_{\text{at}}}=\frac{1}{{{\omega}_{z}}}{\int}_{\omega}\underset{\begin{array}{c} \alpha =\pm \\ i=1,..,N \end{array}}{\sum}\,\alpha \left[{{\lambda}_{\alpha i}}\left(\phi _{\alpha i}^{2}-1\right)-{{\left({{\partial}_{t}}{{\phi}_{\alpha i}}\right)}^{2}}\right]. \end{array}\end{eqnarray} \tag{ 151 }$$

The action for the atom–photon coupling and the bare photonic part is straightforwardly derived according to the previous sections, and the corresponding Keldysh action is

$$\begin{eqnarray}S=\frac{1}{{{\omega}_{z}}}{{{\int}^{}}_{t}}\underset{i}{\sum}\,\left({{\phi}_{ci}},{{\phi}_{qi}}\right)\,\left(\begin{array}{c} {{\lambda}_{qi}}(t) & {{\lambda}_{ci}}(t)+\partial _{t}^{2} \\ {{\lambda}_{ci}}(t)+\partial _{t}^{2} & {{\lambda}_{qi}}(t) \end{array}\right)\,\left(\begin{array}{c} {{\phi}_{ci}} \\ {{\phi}_{qi}} \end{array}\right)\\ &&+\frac{g}{\sqrt{N}}{{{\int}^{}}_{t}}\underset{i}{\sum}\,\left({{\phi}_{ci}}\left(a_{q}^{\ast}+{{a}_{q}}\right)+{{\phi}_{qi}}\left(a_{c}^{\ast}+{{a}_{c}}\right)\right)\\ +{{{\int}^{}}_{t}}\left(a_{c}^{\ast},a_{q}^{\ast}\right)\left(\begin{array}{c} 0 & \text{i}{{\partial}_{t}}-{{\omega}_{\text{p}}}-\text{i}\kappa \\ \text{i}{{\partial}_{t}}-{{\omega}_{\text{p}}}+\text{i}\kappa & 2\text{i}\kappa \end{array}\right)\left(\begin{array}{c} {{a}_{c}} \\ {{a}_{q}} \end{array}\right).\end{eqnarray} \tag{ 152 }$$

It is quadratic in the fields $\phi,a,{{a}^{\ast}}$ as was the case for the Holstein–Primakoff action in the large-N limit. However, the non-linear constraint, i.e. the fact that ${{\lambda}_{q,c}}$ is a dynamic variable, introduces higher order terms in the Ising fields. The approach corresponding to the large-N expansion of the Holstein–Primakoff fields in the present formalism (in the sense that the action becomes quadratic) is to treat the Lagrange parameters as static mean field variables, i.e. ${{\lambda}_{ci}}(t)\to {{\lambda}_{c}}$ and ${{\lambda}_{qi}}(t)\to {{\lambda}_{q}}$. Due to causality, a static mean quantum field must be zero, ${{\lambda}_{q}}=0$. On the other hand, the value of ${{\lambda}_{c}}$ has to be chosen such that the non-linear constraint is preserved on average

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \langle {{\left(\sigma _{i}^{x}(t)\right)}^{2}}\rangle =\langle \phi _{c}^{2}(t)\rangle ={{{\int}^{}}_{\omega}}\text{i}G_{\text{at}}^{K}(\omega )\overset{{}}{{\overset{!}{{=}}\,}}\,1. \end{array}\end{eqnarray} \tag{ 153 }$$

This is equivalent to a vanishing variation of the action with respect to the Lagrange parameters ${{\lambda}_{c,q}}$ and is known as the 'soft-spin' approach to spin models (i.e. the constraint is treated on the mean field level). It becomes exact in the thermodynamic limit $N\to \infty$ [256, 257].

Integrating out the N individual atoms leads to the effective photon action

$$\begin{eqnarray}S=\frac{1}{2}{{{\int}^{}}_{\omega}}{{A}^{\dagger}}(\omega )\left(\begin{array}{c} 0 & {{\left(G_{2\times 2}^{A}\right)}^{-1}} \\ {{\left(G_{2\times 2}^{R}\right)}^{-1}} & 2\text{i}\kappa \,\text{diag}(1,1) \end{array}\right)\,A(\omega ),\end{eqnarray} \tag{ 154 }$$

with the Nambu spinor

$$\begin{eqnarray}&&A(\omega )=\left(\begin{array}{c} {{a}_{c}}(\omega ) \\ a_{c}^{\ast}(-\omega ) \\ {{a}_{q}}(\omega ) \\ a_{q}^{\ast}(-\omega ) \end{array}\right)\end{eqnarray} \tag{ 155 }$$

and the inverse retarded Green's function

$$\begin{eqnarray}{{\left(G_{2\times 2}^{R}\right)}^{-1}}=\left(\begin{array}{c} \omega -{{\omega}_{\text{p}}}+\text{i}\kappa -\frac{2{{g}^{2}}{{\omega}_{z}}}{{{\omega}^{2}}-\lambda} & -\frac{2{{g}^{2}}{{\omega}_{z}}}{{{\omega}^{2}}-\lambda} \\ -\frac{2{{g}^{2}}{{\omega}_{z}}}{{{\omega}^{2}}-\lambda} & \omega -{{\omega}_{\text{p}}}-\text{i}\kappa -\frac{2{{g}^{2}}{{\omega}_{z}}}{{{\omega}^{2}}-\lambda} \end{array}\right).\end{eqnarray} \tag{ 156 }$$

The poles of the Green's function are again determined by the zeros of the determinant of the inverse Green's function, i.e. by the roots of the equation

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} 0\overset{{}}{{\overset{!}{{=}}\,}}\,\left({{\lambda}_{c}}-{{\omega}^{2}}\right)\left({{\left(\omega +\text{i}\kappa \right)}^{2}}-\omega _{\text{p}}^{2}\right)+4{{\omega}_{\text{p}}}{{\omega}_{z}}{{g}^{2}}. \end{array}\end{eqnarray} \tag{ 157 }$$

For the non-interacting theory, i.e. for g  =  0, this determines the poles of the atomic sector to be $\omega =\pm \sqrt{{{\lambda}_{c}}}$ and in turn fixes ${{\lambda}_{c}}=\omega _{z}^{2}$ to be consistent with the microscopic theory. As long as the system is not superradiant, the atom–photon interaction does not modify the integral (153), and for the entire parameter range of the normal phase, we find ${{\lambda}_{c}}=\omega _{z}^{2}$ as for the non-interacting case. This again yields the critical value of the coupling strength (143) we obtained above from the Holstein–Primakoff approach. Consequently, the results from the previous section are recovered in the Ising representation of the Keldysh path integral.

We now turn to the description of the ordered phase. Directly at the superradiant transition, the system becomes unstable, which is indicated by a critical mode with frequency $\omega =0$. In order to stabilize the system, for coupling strengths g  >  g_c the steady state breaks the Z₂ symmetry of the action, which is expressed in terms of a finite, symmetry breaking order parameter ${{\phi}_{c}}(\omega )\to \psi \delta (\omega )+{{\phi}_{c}}(\omega )$, and equivalently in the photon basis ${{a}_{c}}(\omega )\to a\delta (\omega )+{{a}_{c}}(\omega )$. These order parameters correspond to finite expectation values $\langle \sigma _{i}^{x}(t)\rangle =\psi$ and $\langle a(t)\rangle =a$ in the system's ground state. They modify the soft-spin condition (153) according to

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \langle {{\left(\sigma _{i}^{x}(t)\right)}^{2}}\rangle ={{\psi}^{2}}+{{{\int}^{}}_{\omega}}\text{i}G_{\text{at}}^{K}(\omega )\overset{{}}{{\overset{!}{{=}}\,}}\,1, \end{array}\end{eqnarray} \tag{ 158 }$$

where $G_{\text{at}}^{K}$ is the Keldysh Green's function of the fluctuating variable ϕ. As a consequence, in order to fulfill equation (158), the Lagrange multiplier ${{\lambda}_{c}}$ becomes a continuous function of the order parameter, and therefore an implicit function of the coupling g in the superradiant phase. The dependence of ${{\lambda}_{c}}$ can be determined from equation (157). At the superradiant transition at g  =  g_c, this equation has one critical solution, i.e. the propagator has a pole at $\omega =0$, see figure 7. For larger coupling strengths, this pole crosses the origin and obtains a positive imaginary part, thereby rendering the system unstable. In order to compensate for this mechanism, ${{\lambda}_{c}}$ is modified such that the pole does not cross the real axis, i.e. remains at its value $\omega =0$ in the superradiant phase. Consequently ${{\lambda}_{c}}=\frac{4{{g}^{2}}{{\omega}_{\text{p}}}{{\omega}_{z}}}{{{\kappa}^{2}}+\omega _{\text{p}}^{2}}$ for g  >  g_c. Via equation (158), this reveals the scaling of $\psi \sim \sqrt{|g-{{g}_{c}}|}$ when the transition is approached from inside the superradiant phase, which is the same critical behavior of the order parameter as for the equilibrium transition.

The results on the critical properties of the Ising field theory at the superradiance transition conclude the discussion of the single mode Dicke model in the framework of the Keldysh path integral. We have shown that both the Holstein–Primakoff approach in terms of complex bosonic variables and the Ising approach in terms of real Ising variables are equivalent on the level of the quadratic, large-N limit of the theory. The critical point of the transition, as well as the mode structure and the critical scaling behavior have been directly derived from the quadratic Keldysh action, which illustrates the strength of the present field theory approach. In the following section, we discuss the extension of the Dicke model to multiple cavity modes, which contains the possibility of frustrated atom–atom interactions and the formation of a spin glass in the cavity system.

The corresponding Keldysh action is similar to equation (152) and reads

$$\begin{eqnarray}S=\frac{1}{{{\omega}_{z}}}{{{\int}^{}}_{t}}\underset{s}{\sum}\,\left({{\phi}_{cs}},{{\phi}_{qs}}\right)\,\left(\begin{array}{c} 0 & {{\lambda}_{c}}+\partial _{t}^{2} \\ {{\lambda}_{c}}+\partial _{t}^{2} & 0 \end{array}\right)\left(\begin{array}{c} {{\phi}_{cs}} \\ {{\phi}_{qs}} \end{array}\right)\\ &&+{{{\int}^{}}_{t}}\underset{sl}{\sum}\,{{g}_{sl}}\left({{\phi}_{cs}}\left(a_{ql}^{\ast}+{{a}_{ql}}\right)+{{\phi}_{qs}}\left(a_{cl}^{\ast}+{{a}_{cl}}\right)\right)\\ +{{{\int}^{}}_{t}}\underset{l}{\sum}\,\left(a_{cl}^{\ast},a_{ql}^{\ast}\right)\,\left(\begin{array}{c} 0 & \text{i}{{\partial}_{t}}-{{\omega}_{l}}-\text{i}\kappa \\ \text{i}{{\partial}_{t}}-{{\omega}_{l}}+\text{i}\kappa & 2\text{i}\kappa \end{array}\right)\left(\begin{array}{c} {{a}_{cl}} \\ {{a}_{ql}} \end{array}\right).\end{eqnarray} \tag{ 161 }$$

Here, the soft spin approximation has been performed and the scaling of the coupling constants g_sl in the thermodynamic limit is implicit, i.e. ${{g}_{sl}}\sim {{N}^{-\frac{1}{2}}}$. The thermodynamic limit is reached by taking an extensive number of atoms $N\to \infty$ but leaving the number of photon modes $M<\infty$ finite since the consideration of an extensive number of photon modes is physically unrealistic. In order to obtain a large-N effective action, the photon degrees of freedom are integrated out, leading to the atomic action in frequency space

$$\begin{eqnarray}S=\frac{1}{{{\omega}_{z}}}{{{\int}^{}}_{\omega}}\underset{s}{\sum}\,\left({{\phi}_{cs}},{{\phi}_{qs}}\right)\left(\begin{array}{c} 0 & {{\lambda}_{c}}-{{\omega}^{2}} \\ {{\lambda}_{c}}-{{\omega}^{2}} & 0 \end{array}\right)\left(\begin{array}{c} {{\phi}_{cs}} \\ {{\phi}_{qs}} \end{array}\right)\\ +\underset{si}{\sum}\,{{J}_{si}}{{{\int}^{}}_{\omega}}\left({{\phi}_{cs}},{{\phi}_{qs}}\right)\left(\begin{array}{c} 0 & {{\Lambda}^{A}}(\omega ) \\ {{\Lambda}^{R}}(\omega ) & {{\Lambda}^{K}}(\omega ) \end{array}\right)\left(\begin{array}{c} {{\phi}_{ci}} \\ {{\phi}_{qi}} \end{array}\right).\end{eqnarray} \tag{ 162 }$$

The frequency dependent terms are the symmetrized photon Green's functions

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{\Lambda}^{R}}(\omega ) & =\frac{-{{\omega}_{\text{p}}}}{{{\left(\omega +\text{i}\kappa \right)}^{2}}-\omega _{\text{p}}^{2}}, \\ {{\Lambda}^{K}}(\omega ) & =\frac{2\text{i}\kappa \left({{\omega}^{2}}+{{\kappa}^{2}}+\omega _{\text{p}}^{2}\right)}{|{{\left(\omega +\text{i}\kappa \right)}^{2}}-\omega _{\text{p}}^{2}{{|}^{2}}} \end{array}\end{eqnarray} \tag{ 163 }$$

with an average photon frequency of the cavity ${{\omega}_{\text{p}}}$. The effective coupling

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{J}_{si}}=\underset{l=1}{\overset{M}{\sum}}\,\frac{{{g}_{sl}}{{g}_{il}}}{4} \end{array}\end{eqnarray} \tag{ 164 }$$

fluctuates between the values $-\frac{M{{g}_{0}}}{4}\leqslant {{J}_{si}}\leqslant \frac{Mg_{0}^{2}}{4}$ and can be either ferromagnetic J_si  <  0 or antiferromagnetic J_si  >  0 depending on the s, i. Assuming g_sl to be independently and equally distributed for each atom, the couplings J_si are for sufficiently large M distributed according to a Gaussian distribution function. We set their average value $\bar{J}\equiv \langle {{J}_{si}}\rangle =0$ as it lifts the frustration in the system caused by fluctuating couplings but does not modify the universal behavior of the system at the glass transition [113, 131].

The averaged partition function of the system, including the probability distribution for the couplings J_si is

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} Z={\int}^{}\mathcal{D}\left[{{\phi}_{s}}\right]\mathcal{D}\left[{{J}_{si}}\right]{{\text{e}}^{\text{i}S+\text{i}{{S}_{P}}}}, \end{array}\end{eqnarray} \tag{ 165 }$$

where

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \text{i}{{S}_{P}}=-\frac{N}{2}\underset{s,i}{\sum}\,\frac{J_{si}^{2}}{K} \end{array}\end{eqnarray} \tag{ 166 }$$

is the action of the random couplings J_si, with correlations

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \langle {{J}_{si}}{{J}_{lm}}\rangle =\left({{\delta}_{s,l}}{{\delta}_{i,m}}+{{\delta}_{s,m}}{{\delta}_{i,l}}\right)\,\frac{K}{N}. \end{array}\end{eqnarray} \tag{ 167 }$$

In this sense, the coupling of the atomic degrees of freedom to the random variables J_si has the same structure as the coupling of the atoms to an external bath. However, the significant difference to a Markovian bath—which would be δ-correlated in space and time—is that the quenched disorder, while correlated locally in space, has infinite correlation in time. The latter is expressed by the fact that the variables J_si are time-independent.

Averaging over all realizations of the couplings is done by a Gaussian integration (see appendix B), and transforms the disorder part, i.e. the second line of equation (162), to a time non-local quartic contribution

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{S}_{\text{at}-\text{at}}}=\frac{\text{i}K}{N}{\int}_{\omega,\nu}\underset{l,m}{\sum}\,\Phi_{l,m}^{\Lambda}(\omega )\Phi_{m,l}^{\Lambda}(\nu ), \end{array}\end{eqnarray} \tag{ 168 }$$

with

$$\begin{eqnarray}\Phi_{l,m}^{\Lambda}(\omega )=\left({{\phi}_{cl}},{{\phi}_{ql}}\right)\left(\begin{array}{c} 0 & {{\Lambda}^{A}}(\omega ) \\ {{\Lambda}^{R}}(\omega ) & {{\Lambda}^{K}}(\omega ) \end{array}\right)\left(\begin{array}{c} {{\phi}_{cm}} \\ {{\phi}_{qm}} \end{array}\right).\end{eqnarray} \tag{ 169 }$$

The double sum in (168) is decoupled via a Hubbard–Stratonovich transformation, which introduces the macroscopic fields ${{Q}_{\alpha {{\alpha}^{\prime}}}}$, with $\alpha,{{\alpha}^{\prime}}=c,q$ being Keldysh indices [113, 131]. This transformation is in general not unique but can be made unique by requiring the equivalence

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \frac{\delta}{\delta {{Q}_{\alpha {{\alpha}^{\prime}}}}}Z=0\Leftrightarrow {{Q}_{\alpha {{\alpha}^{\prime}}}}(\omega,\nu )=\frac{1}{N}\underset{l}{\sum}\,\langle {{\phi}_{\alpha l}}(\omega ){{\phi}_{{{\alpha}^{\prime}}l}}(\nu )\rangle. \end{array}\end{eqnarray} \tag{ 170 }$$

This identifies the Hubbard–Stratonovich field with the average atomic correlation function. Since one is interested in the stationary, time translational invariant state,

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{Q}_{\alpha {{\alpha}^{\prime}}}}(\omega,\nu )=2\pi \delta (\omega +\nu ){{Q}_{\alpha {{\alpha}^{\prime}}}}(\omega ). \end{array}\end{eqnarray} \tag{ 171 }$$

After the Hubbard–Stratonovich decoupling, the resulting action is quadratic in the atomic fields and they are integrated out, leading to the macroscopic action

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \mathcal{S}=\text{i}N\text{Tr}\left[K\Lambda Q\Lambda Q-\frac{1}{2}\ln \tilde{G}\right](\omega ). \end{array}\end{eqnarray} \tag{ 172 }$$

Here $\Lambda$ and Q are matrices in Keldysh space and $\tilde{G}$ is defined as

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \tilde{G}(\omega )={{\left(G_{0}^{-1}(\omega )-2K\Lambda(\omega )Q(\omega )\Lambda(\omega )\right)}^{-1}}. \end{array}\end{eqnarray} \tag{ 173 }$$

In the thermodynamic limit, the partition function is determined by the saddle point value of the action, i.e. by the condition

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \frac{\delta \mathcal{S}}{\delta {{Q}_{\alpha {{\alpha}^{\prime}}}}}\overset{{}}{{\overset{!}{{=}}\,}}\,0 \end{array}\end{eqnarray} \tag{ 174 }$$

for all $\alpha,{{\alpha}^{\prime}}$. This yields the values of the fields ${{Q}_{\alpha {{\alpha}^{\prime}}}}(\omega )$ at the saddle-point, which can be identified with the atomic response and correlation functions. The saddle-point equations for the atomic response and correlation functions are

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{Q}^{R}}(\omega )={{\left(\frac{2\left(\lambda -{{\omega}^{2}}\right)}{{{\omega}_{z}}}-4K{{\left({{\Lambda}^{R}}(\omega )\right)}^{2}}{{Q}^{R}}(\omega )\right)}^{-1}} \end{array}\end{eqnarray} \tag{ 175 }$$

and

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{Q}^{K}}(\omega )=\frac{4K|{{Q}^{R}}{{|}^{2}}{{\Lambda}^{K}}\left({{Q}^{A}}{{\Lambda}^{A}}+{{Q}^{R}}{{\Lambda}^{R}}\right)}{1-4K|{{Q}^{R}}{{\Lambda}^{R}}{{|}^{2}}}. \end{array}\end{eqnarray} \tag{ 176 }$$

Additionally, due to equation (170), the soft-spin constraint maps to the Q-fields according to

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \text{i}{{{\int}^{}}_{\omega}}{{Q}^{K}}(\omega )=\frac{2}{N}\underset{l}{\sum}\,\langle {{\phi}_{cl}}(-\omega ){{\phi}_{cl}}(\omega )\rangle =2. \end{array}\end{eqnarray} \tag{ 177 }$$

In the glass phase, the spins attain temporally frozen configurations, expressed by an infinite correlation time of the atomic correlator, which is expressed by a non-zero Edwards–Anderson parameter

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{q}_{\text{EA}}}\equiv \underset{t\to \infty}{{\lim}}\,\frac{1}{N}\underset{l}{\sum}\,\langle \sigma _{l}^{x}(t)\sigma _{l}^{x}(0)\rangle. \end{array}\end{eqnarray} \tag{ 178 }$$

Consequently, the correlation function ${{Q}^{K}}(\omega )$ consists of a regular part, describing the short time dynamics and a non-vanishing contribution at $\omega =0$. It can be expressed via the modified fluctuation dissipation relation [261]

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{Q}^{K}}(\omega )=4\text{i}\pi {{q}_{\text{EA}}}\delta (\omega )+Q_{r}^{K}(\omega ) \end{array}\end{eqnarray} \tag{ 179 }$$

in terms of the order parameter and a regular contribution $Q_{r}^{K}(\omega )$. The soft-spin condition (177) together with equation (176) fixes the value of the Lagrange parameter λ and therefore fully determines the spectrum of the system. Similar to the superradiant transition, the Edwards–Anderson parameter becomes non-zero when the poles of the system become critical (approach the real axis) and its emergence is a mechanism to stabilize the critical modes in the ordered phase.

The variance of the effective, long-range atom–atom interaction K is a measure of the frustration in the system and for sufficiently strong frustration, the system enters a glass phase, described by an infinite autocorrelation time of the spins and the emergence of a non-zero Edwards–Anderson parameter ${{q}_{\text{EA}}}>0$. This goes hand in hand with a critical continuum of modes reaching zero, which is the characteristic feature of a critical phase of matter. The phase diagram for the fully coherent model has been analyzed in [131] and in [113] it was shown that the glass phase persists in the presence of dissipation. However, the dissipative model shows new universal features of the glass phase, which correspond neither to a zero nor to a finite temperature equilibrium transition. This is attributed to the fact that the critical modes are described by poles in the complex plane which are neither purely real as in the zero temperature (or quantum) case nor purely imaginary as for the finite temperature transition.

In the presence of dissipation, a photon that is emitted from an atom can either be absorbed by another atom, leading to an effective atom–atom interaction, which is of infinite range in space, or decay from the cavity. The latter happens on a characteristic time scale ${{\tau}_{p}}=1/\kappa$, which is as well the time scale for which atoms can interact with each other by exchanging a photon before the photon decays. As a consequence, the photon decay reduces the effective range of the atom–atom interaction in time and reduces the strength of frustration in the system. Close to the glass transition and in the glass phase, this leads to the emergence of a crossover scale ${{\omega}_{c}}=\tau _{c}^{-1}$, above which the spectral properties of the system correspond exactly to a zero temperature spin glass. Below this frequency, the spectrum reveals the breaking of time reversal symmetry by the dissipation and formally corresponds to the dynamical universality class of dissipative quantum glasses, other examples of which are spin glasses coupled to an ohmic bath [261–263] or a bath of metallic electrons [257, 259]. The crossover frequency is

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{\omega}_{\text{c}}}=\kappa \,{{\left(\frac{\omega _{\text{p}}^{2}+\frac{1}{2}{{\kappa}^{2}}}{\omega _{\text{p}}^{2}+{{\kappa}^{2}}}+\frac{{{\left(\omega _{\text{p}}^{2}+{{\kappa}^{2}}\right)}^{2}}}{2\sqrt{K}\omega _{z}^{2}}\right)}^{-1}}. \end{array}\end{eqnarray} \tag{ 180 }$$

and the modifications of the spectrum below this scale are due to the damping introduced by the Markovian bath. On the normal or disordered side of the phase diagram, the corresponding low frequency propagator is

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} Q_{n}^{R}(\omega )=Z{{\left({{\left(\omega +\text{i}\gamma \right)}^{2}}-{{\alpha}^{2}}\right)}^{-1}}, \end{array}\end{eqnarray} \tag{ 181 }$$

describing the Ising spins as damped harmonic oscillators with decay rate $\gamma >0$ and frequency α, with the physical meaning of an inverse lifetime and energy gap of the atomic excitations. When the glass transition is approached, the gap, the inverse lifetime and the residue Z scale to zero simultaneously and the atomic response in the glass phase

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} Q_{g}^{R}=\bar{Z}{{\left({{\omega}^{2}}+\bar{\gamma}|\omega |\right)}^{-\frac{1}{2}}} \end{array}\end{eqnarray} \tag{ 182 }$$

is non-analytic and can no longer be interpreted as a harmonic oscillator. The broken time reversal invariance manifests itself in both parameters $\gamma,\bar{\gamma}\ne 0$ being non-zero, which modifies the low frequency dynamics in the entire glass phase towards a dissipative quantum glass. This is for instance expressed by a spectral density which features an anomalous square root behavior $\mathcal{A}(\omega )=-2\text{Im}\left({{Q}^{R}}(\omega )\right)\sim \sqrt{|\omega |}$ for small frequencies and therefore has a non-analytic response at zero frequency, see figure 8. The latter leads to a characteristic algebraic decay of the photon correlation function, as discussed below. The universality class of the present dissipative quantum glass transition is determined by the critical exponents at the transition, which describe the scaling of the parameters $\alpha,\gamma,Z$ as a function of $\delta =\,|K-{{K}_{c}}|$. It is summarized in the equations

$$\begin{eqnarray}\begin{array}{c} {{\alpha}_{\delta}} & = & \frac{\sqrt{2}\left(\omega _{\text{p}}^{2}+{{\kappa}^{2}}\right)}{8\sqrt{{{K}^{3}}}\kappa} & |\frac{\delta}{\log (\delta )}{{|}^{\frac{3}{2}}} \\ {{\gamma}_{\delta}} & = & \frac{\omega _{\text{p}}^{2}+{{\kappa}^{2}}}{16{{K}^{2}}\kappa} & |\frac{\delta}{\log (\delta )}{{|}^{2}} \\ {{Z}_{\delta}} & = & \frac{{{\omega}_{\text{p}}}\left(\omega _{\text{p}}^{2}+{{\kappa}^{2}}\right)}{8\sqrt{{{K}^{5}}}{{\kappa}^{2}}} & |\frac{\delta}{\log (\delta )}{{|}^{3}} \end{array},\end{eqnarray} \tag{ 183 }$$

which show the typical logarithmic correction to scaling at a quantum glass transition and identify the critical exponents. Indicated by the scaling behavior, dissipative and coherent dynamics rival each other when approaching the transition, separating this glass transition from equilibrium transitions and the previously discussed Dicke transition, which are either fully coherent (quantum) with $\gamma =0$ or fully dissipative (classical) with $\alpha =0$ sufficiently close to the transition. The present glass transition has therefore no counterpart in static equilibrium physics and is termed dissipative quantum glass transition. Similar behavior is present in system bath settings in equilibrium, where the bath however not only imprints a finite temperature to the system but as well modifies its spectral properties [259, 261, 262]. What these situations share in common is that the effective theory, after elimination of the bath variables, obeys no time reversal symmetry, which ensures the same asymptotic universal behavior as in the case of the Markovian photon loss in the present system.

Zoom In Zoom Out Reset image size

As typical for many open quantum systems (for exceptions, see the discussion in section 4), the statistical properties of the excitations are described by a low energy effective temperature (LET) and a corresponding thermal distribution of the excitations. However, as has been found for the single mode Dicke model in the normal as well as in the superradiant phase, the LET for the photonic and the atomic subsystems did not coincide, i.e. the subsystems have not thermalized but are held at different temperatures corresponding to their individual coupling to a bath [112]. This dramatically changes in the present system with multiple photon modes. As frustration is increased by driving the variance K towards the critical value K_c, atoms and photons begin to thermalize towards the same, shared LET. The distribution function F of photonic or atomic modes is obtained by solving the fluctuation dissipation relation

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{Q}^{K}}(\omega )={{Q}^{R}}(\omega ){{F}_{\text{at}}}(\omega )-{{F}_{\text{at}}}(\omega ){{Q}^{A}}(\omega ) \end{array}\end{eqnarray} \tag{ 184 }$$

for the atomic Green's functions. It leads to an atomic LET

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{T}_{\text{eff}}}=\frac{\omega _{\text{p}}^{2}+{{\kappa}^{2}}}{4{{\omega}_{\text{p}}}} \end{array}\end{eqnarray} \tag{ 185 }$$

in the glass phase, which coincides with the photonic LET [113]. The thermalization of the two subsystems is a consequence of the disorder induced effective long range interactions, which redistribute energy between the different modes and enable equilibration even in the presence of the Markovian dissipation in the photon sector. In the paramagnetic phase, the distribution functions of atoms and photons are identical for frequencies $\omega >\alpha$ larger than the gap, but deviate from each other for lower frequencies. The elementary excitations above this frequency are strongly correlated and can not be seen as weakly dressed photons or atoms. As a consequence, the observables in the critical glass phase are dominated by the universal low energy behavior of the glass propagator (182) and thermal statistics of the excitations.

The strong light–matter interactions not only lead to thermalization of the atomic and photonic subsystems, they also feature a complete imprint of the glass dynamics of the atoms onto the cavity photons. This results in universal spectral properties of the photonic response, as has been shown in [113], and concerns a complete locking of both subsystems' low energy response properties. The photons form a photon glass state, which features a large, incoherent population of photon modes, signaled by a finite Edwards–Anderson parameter in the photon correlation function

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} G_{\text{ph}}^{K}(\omega )=4\pi \delta (\omega ){{{\tilde{q}}}_{\text{EA}}}+G_{\text{reg}}^{K}(\omega ), \end{array}\end{eqnarray} \tag{ 186 }$$

which shows the same scaling behavior as the atomic Edwards–Anderson parameter close to the glass transition ${{\tilde{q}}_{\text{EA}}}\sim {{q}_{\text{EA}}}$. A finite ${{\tilde{q}}_{\text{EA}}}$ implies long temporal memory in photon autocorrelation functions and an extensive number of photons permanently occupying a continuum of high energy modes. The latter is revealed by an algebraic decay of the system's correlation functions, as for instance for the four point (or g⁽²⁾) correlation function, which is shown in figure 9.

Zoom In Zoom Out Reset image size

The photon response and correlation functions are easily accessible via cavity photon output measurements, and therefore represent the natural observables for the detection of the glass dynamics in cavity QED experiments. In fact, in [113], it has been demonstrated that a complete characterization of the glass phase can be performed by different but standard photon output measurements, such as homodyne detection and intensity correlation measurements. Formally, the photon Green's functions can be obtained from the Keldysh action formalism by introducing source fields ${{\mu}_{\alpha {{\alpha}^{\prime}}}}$ in the microscopic action, which couple to the photon variables. After integrating out the photon variables in the action, the photon Green's functions are then obtained via functional derivatives with respect to these source fields. This standard field theory technique relates the photon Green's functions to the solution for the atomic response and correlation functions [113].

The Keldysh path integral formalism represents a theoretical approach which incorporates in a very straightforward way the coupling of the system to a non-thermal bath. Here we compare quenched disorder and Markovian baths in more detail. As shown above in the present section, the coupling of the system variables to quenched disorder, realized by the variables J_lm, is equivalent to the coupling to a bath with infinite correlation time,

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \langle {{J}_{lm}}(t){{J}_{{{l}^{\prime}}{{m}^{\prime}}}}\left({{t}^{\prime}}\right)\rangle =\frac{K}{N}\left({{\delta}_{l{{l}^{\prime}}}}{{\delta}_{m{{m}^{\prime}}}}+{{\delta}_{l{{m}^{\prime}}}}{{\delta}_{m{{l}^{\prime}}}}\right). \end{array}\end{eqnarray} \tag{ 187 }$$

This represents the opposite limit of a Markovian bath, which has a typical correlation time that is much shorter than the system correlation time. As a consequence, at each single system–bath interaction event, the Markovian bath immediately equilibrates and loses all its memory on the interaction. On the other hand, the quenched bath equilibrates infinitely slowly and keeps memory on each single system–bath interaction. As a consequence, it introduces effective temporally long range interactions between the system variables. These tend to slow down the system, pronouncing its $\omega =0$ correlations as expressed by the Edwards–Anderson parameter.

In an equilibrium formalism, disordered systems have to be dealt with by a computationally demanding replica method [256, 257], which is furthermore physically far less transparent than the Keldysh formalism, which thus represents a convenient approach to disordered systems. In open quantum systems, where disorder, i.e. the coupling to a quenched bath, comes together with the dissipation introduced by a Markovian bath, both types of bath couplings compete with each other, leading in the present glassy system to strongly modified spectral properties, which mirror the effect of the Markovian bath by a crossover from non-equilibrium to equilibrium scaling behavior.

The Wetterich equation (132) describes how the scale-dependent effective action ${{\Gamma}_{\Lambda}}$ evolves from the microscopic action S at $\Lambda={{\Lambda}_{0}}$ to the full effective action $\Gamma$ at $\Lambda\to 0$, as fluctuations with momenta larger than the cutoff $\Lambda$ are integrated out, and the latter is gradually lowered. It is an exact non-linear differential equation for the functional ${{\Gamma}_{\Lambda}}$. In the absence of an exact solution, suitable schemes to approximate the functional differential equation have to be found. For the analysis of critical behavior, progress can be made by choosing an ansatz for ${{\Gamma}_{\Lambda}}$ that has a similar structure to the microscopic action S in equation (73). In this way, the effective action ${{\Gamma}_{\Lambda}}$ is parameterized in terms of the coupling constants appearing in the ansatz, and consequently the functional differential equation for ${{\Gamma}_{\Lambda}}$ can be cast in the form of a set of ordinary differential equations for the couplings, as we describe in detail below.

Any such ansatz will contain only a finite number of couplings, and will hence truncate the most general structure of ${{\Gamma}_{\Lambda}}$. In other words, by choosing an ansatz of a particular form, one makes an approximation, and the question arises, which couplings should be included in the ansatz or truncation in order to obtain meaningful results. For the study of critical phenomena at a second order phase transition, the power counting scheme of section 2.3 provides a guideline: following the arguments given there, we choose a truncation which includes all couplings that are not irrelevant, and which therefore takes the form of the semiclassical action in equation (79):

$$\begin{eqnarray}&&{{\Gamma}_{\Lambda}}={{{\int}^{}}_{t,\mathbf{x}}}\left\{\bar{\phi}_{q}^{\ast}\left[\left(\text{i}{{Z}^{\ast}}{{\partial}_{t}}+{{\bar{K}}^{\ast}}{{\nabla}^{2}}\right){{\bar{\phi}}_{c}}-\frac{\partial {{{\bar{U}}}^{\ast}}}{\partial \bar{\phi}_{c}^{\ast}}\right]+\text{c}.\text{c}.+\text{i}\bar{\gamma}\bar{\phi}_{q}^{\ast}{{\bar{\phi}}_{q}}\right\}.\end{eqnarray} \tag{ 204 }$$

(Recall that the variables of the effective action are the field expectation values ${{\rm{\bar \Phi}}_{\nu}},\nu =c,q$ defined in equation (37).) Here, in addition to the complex prefactor $\bar{K}=\bar{A}+\text{i}\bar{D}$ of the Laplacian, we are including a complex wave-function renormalization $Z={{Z}_{R}}+\text{i}{{Z}_{I}}$. In the semiclassical limit, the homogeneous (i.e. not containing derivatives with respect to time or space) part of the action can be written in terms of an effective potential $\bar{U}$, which is a function of ${{\bar{\rho}}_{c}}=\bar{\phi}_{c}^{\ast}{{\bar{\phi}}_{c}}$. Due to the invariance of the microscopic action under classical phase rotations (see section 2.4.4), which is inherited by the effective action, only this combination of fields is allowed in the potential. The latter is given by

$$\begin{eqnarray}&&\bar{U}\left({{\bar{\rho}}_{c}}\right)=\frac{1}{2}{{\bar{u}}_{2}}{{\left({{\bar{\rho}}_{c}}-{{\bar{\rho}}_{0}}\right)}^{2}}+\frac{1}{6}{{\bar{u}}_{3}}{{\left({{\bar{\rho}}_{c}}-{{\bar{\rho}}_{0}}\right)}^{3}}.\end{eqnarray} \tag{ 205 }$$

Here, both the two-body and three-body couplings, ${{\bar{u}}_{2}}=\bar{\lambda}+\text{i}\bar{\kappa}$ and ${{\bar{u}}_{3}}={{\bar{\lambda}}_{3}}+\text{i}{{\bar{\kappa}}_{3}}$ respectively, are complex. The three-body term is marginal according to power counting, and therefore included in the truncation. In the FRG, it is advantageous to approach the transition from the ordered phase. Then, the form of the effective potential in equation (205) corresponds to an expansion around the stationary condensate density ${{\bar{\rho}}_{0}}$. Indeed, this choice implies that the field equations $\delta {{\Gamma}_{\Lambda}}/\delta \bar{\phi}_{c}^{\ast}=0,\delta {{\Gamma}_{\Lambda}}/\delta \bar{\phi}_{q}^{\ast}=0$ (equations (40) in the absence of sources and evaluated with the scale-dependent effective action ${{\Gamma}_{\Lambda}}$) are solved by ${{\bar{\rho}}_{c}}={{\bar{\rho}}_{0}}$ and ${{\bar{\phi}}_{q}}=0$ on all scales $\Lambda$.

In the truncation equation (204), all couplings—including the condensate density ${{\bar{\rho}}_{0}}$—are scale dependent. As indicated above, by means of such an ansatz for the effective action ${{\Gamma}_{\Lambda}}$, the functional differential equation (132) can be rewritten as a set of ordinary differential equations for these running couplings, with initial conditions given by the microscopic action equation (79). This is achieved by applying projection prescriptions. In the following, we summarize this method for the problem at hand. For details we refer the reader to [181].

The main idea of a projection prescription on a specific coupling is to extract this coupling from the effective action ${{\Gamma}_{\Lambda}}$ by taking appropriate derivatives of the latter with respect to the fields and coordinates, and subsequently setting the fields to their stationary values ${{\bar{\phi}}_{c}}={{\bar{\phi}}_{0}}=\sqrt{{{{\bar{\rho}}}_{0}}}$ and ${{\bar{\phi}}_{q}}=0$ (as noted in section 2.4.6, choosing ${{\bar{\phi}}_{c}}$ to be real does not lead to a loss of generality). Then, applying the very same projection description to the Wetterich equation (132) yields the flow equation for the corresponding coupling. For the actual evaluation of the resulting flow equations, it is convenient to introduce rescaled fields, which are related to the bare ones by absorbing the wave-function renormalization Z in the quantum field:

$$\begin{eqnarray}&&{{\phi}_{c}}={{\bar{\phi}}_{c}},\quad {{\phi}_{q}}=Z{{\bar{\phi}}_{q}}.\end{eqnarray} \tag{ 206 }$$

As a consequence of this transformation, and by rescaling all couplings appropriately, it is possible to obtain a reduced set of flow equations, from which Z is eliminated. In a second step, the number of flow equations can be diminished further by introducing dimensionless renormalized variables, as detailed in section 4.2.3 below.

We start by deriving flow equations for the non-linear couplings in the effective potential defined in equation (205). To see how one can project the Wetterich equation onto flow equations for these couplings, consider the effective action, evaluated for homogeneous, i.e. space- and time-independent 'background fields':

$$\begin{eqnarray}&&{{\Gamma}_{\Lambda,cq}}=-\Omega \left({{\bar{U}}^{\prime}}{{\bar{\rho}}_{cq}}+{{\bar{U}}^{\prime \ast}}{{\bar{\rho}}_{qc}}-\text{i}\bar{\gamma}{{\bar{\rho}}_{q}}\right).\end{eqnarray} \tag{ 207 }$$

Here, the subscript cq in ${{\Gamma}_{\Lambda,cq}}$ indicates that both the classical and the quantum fields are set to constant but non-zero values; $\Omega$ denotes the quantization volume, and we introduced the following products of fields, which are invariant under classical phase rotations (see section 2.4.4): ${{\bar{\rho}}_{cq}}=\bar{\phi}_{c}^{\ast}{{\bar{\phi}}_{q}}=\bar{\rho}_{qc}^{\ast}$ and ${{\bar{\rho}}_{q}}=\bar{\phi}_{q}^{\ast}{{\bar{\phi}}_{q}}$. From the representation equation (207) it becomes immediately clear how to project the flow equation (132) for ${{\Gamma}_{\Lambda}}$ onto a flow equation for the derivative of the potential $\bar{U}$ in equation (205) with respect to ${{\bar{\rho}}_{c}}$: one has to (i) evaluate equation (132) for homogeneous fields, (ii) take the derivative with respect to ${{\bar{\rho}}_{cq}}$, and finally (iii) set the quantum background fields to their stationary state value,

$$\begin{eqnarray}&&{{\partial}_{\ell}}{{\bar{U}}^{\prime}}=-\frac{1}{\Omega }{{\left[{{\partial}_{{{{\bar{\rho}}}_{cq}}}}{{\partial}_{\ell}}{{\Gamma}_{\Lambda,cq}}\right]}_{{{{\bar{\phi}}}_{q}}=\bar{\phi}_{q}^{\ast}=0}}.\end{eqnarray} \tag{ 208 }$$

Here and in the following, we specify flow equations in terms of the logarithmic cutoff scale $\ell =\ln \left(\Lambda/{{\Lambda}_{0}}\right)$. From the potential $\bar{U}$ in equation (205), the couplings ${{\bar{u}}_{2,3}}$ can be obtained by taking further derivatives with respect to ${{\bar{\rho}}_{c}}$. However, instead of projecting the flow equation (208) in this way onto equations for ${{\bar{u}}_{2,3}}$, it is more convenient to introduce rescaled quantities as outlined above. Inserting the representation equation (206) of the bare quantum field in the effective action (207) leads to appearance of a factor 1/Z^* in front of the term involving the effective potential. This factor can be absorbed by introducing a rescaled potential via $\bar{U}=ZU$. The flow equations of the bare and rescaled effective potential are related via

$$\begin{eqnarray}&&{{\partial}_{\ell}}{{\bar{U}}^{\prime}}=Z\left(-{{\eta}_{Z}}{{U}^{\prime}}+{{\partial}_{\ell}}{{U}^{\prime}}\right),\end{eqnarray} \tag{ 209 }$$

where ${{\eta}_{Z}}$ denotes the anomalous dimension associated with the wave-function renormalization (${{\partial}_{\ell}}Z$ is specified below in equation (222)),

$$\begin{eqnarray}&&{{\eta}_{Z}}=-{{\partial}_{\ell}}Z/Z.\end{eqnarray} \tag{ 210 }$$

Then, with ${{\partial}_{{{{\bar{\rho}}}_{cq}}}}=Z{{\partial}_{{{\rho}_{cq}}}}$ and inserting equation (208) on the RHS of equation (211), the flow equation for the renormalized potential becomes (here, primes denote derivatives with respect to ${{\rho}_{c}}=\phi _{c}^{\ast}{{\phi}_{c}}$)

$$\begin{eqnarray}&&{{\partial}_{t}}{{U}^{\prime}}={{\eta}_{Z}}{{U}^{\prime}}+{{\zeta}^{\prime}},\quad {{\zeta}^{\prime}}=-\frac{1}{\Omega }{{\left[{{\partial}_{{{\rho}_{cq}}}}{{\partial}_{t}}{{\Gamma}_{k,cq}}\right]}_{{{\phi}_{q}}=\phi _{q}^{\ast}=0}}.\end{eqnarray} \tag{ 211 }$$

In analogy to equation (205), the renormalized effective potential can be written as

$$\begin{eqnarray}&&U\left({{\rho}_{c}}\right)=\frac{1}{2}{{u}_{2}}{{\left({{\rho}_{c}}-{{\rho}_{0}}\right)}^{2}}+\frac{1}{6}{{u}_{3}}{{\left({{\rho}_{c}}-{{\rho}_{0}}\right)}^{3}},\end{eqnarray} \tag{ 212 }$$

with renormalized couplings defined as ${{u}_{2}}={{\bar{u}}_{2}}/Z=\lambda +\text{i}\kappa$ and ${{u}_{3}}=\bar{u}/Z={{\lambda}_{3}}+\text{i}{{\kappa}_{3}}$. To obtain flow equations for u_n, n  =  2, 3, one simply has to take derivatives of the relation ${{u}_{n}}={{U}^{(n)}}\left({{\rho}_{0}}\right)$ with respect to the cutoff scale $\ell$, taking into account that also ${{\rho}_{0}}$ is a running coupling:

$$\begin{eqnarray}&&{{\partial}_{\ell}}{{u}_{n}}=\left({{\partial}_{\ell}}{{U}^{(n)}}\right)\left({{\rho}_{0}}\right)+{{U}^{(n+1)}}\left({{\rho}_{0}}\right){{\partial}_{\ell}}{{\rho}_{0}}.\end{eqnarray} \tag{ 213 }$$

Inserting equation (211) on the RHS of this relation leads us to

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{\partial}_{\ell}}{{u}_{2}} & ={{\beta}_{{{u}_{2}}}}={{\eta}_{Z}}{{u}_{2}}+{{u}_{3}}{{\partial}_{\ell}}{{\rho}_{0}}+{{\partial}_{{{\rho}_{c}}}}{{\zeta}^{\prime}}{{|}_{\text{ss}}}, \end{array}\end{eqnarray} \tag{ 214 }$$

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{\partial}_{\ell}}{{u}_{3}} & ={{\beta}_{{{u}_{3}}}}={{\eta}_{Z}}{{u}_{3}}+\partial _{{{\rho}_{c}}}^{2}{{\zeta}^{\prime}}{{|}_{\text{ss}}}, \end{array}\end{eqnarray} \tag{ 215 }$$

where we evaluate ${{\zeta}^{\prime}}$ with ${{\rho}_{c}}$ set to its stationary value ${{\rho}_{c}}{{|}_{\text{ss}}}={{\rho}_{0}}$. Finally, flow equations for the real and imaginary parts of u₂ and u₃ can be obtained by taking the real and imaginary parts of equation (214) and (215) respectively.

The flow equation for the stationary density ${{\rho}_{0}}$ cannot be specified without formal ambiguity: as we have already seen in section 2.2.2, both real and imaginary parts of the field equation (75) yield conditions on ${{\rho}_{0}}$. However, we have seen as well that the condition stemming from the real part can always be satisfied by choosing the proper rotating frame. Therefore, the physically correct choice is to assume that ${{\rho}_{0}}$ is actually determined by the imaginary part of the field equation, i.e. by the condition $\text{Im}{{U}^{\prime}}\left({{\rho}_{0}}\right)=0$. Taking the derivative of this condition with respect to the cutoff $\ell$, we find

$$\begin{eqnarray}&&{{\partial}_{\ell}}{{\rho}_{0}}=-\left(\text{Im}{{\partial}_{\ell}}{{U}^{\prime}}\right)\left({{\rho}_{0}}\right)/\text{Im}{{U}^{\prime\prime}}\left({{\rho}_{0}}\right)=-\text{Im}{{\zeta}^{\prime}}{{|}_{\text{ss}}}/\kappa,\end{eqnarray} \tag{ 216 }$$

where ${{\zeta}^{\prime}}$ is the same as in equation (211).

Having illustrated the main idea, the flow equation for the rescaled noise strength $\gamma =\bar{\gamma}/|Z{{|}^{2}}$ can easily be obtained along the lines of the derivation of equations (214) and (215), with the result (for details of the derivation see [181]; ${{\eta}_{ZR}}$ is the real part of ${{\eta}_{Z}}$)

$$\begin{eqnarray}&&{{\partial}_{\ell}}\gamma ={{\beta}_{\gamma}}=2{{\eta}_{ZR}}\gamma -\frac{\text{i}}{\Omega }{{\left[{{\partial}_{{{\rho}_{q}}}}{{\partial}_{\ell}}{{\Gamma}_{\Lambda,cq}}\right]}_{\text{ss}}}.\end{eqnarray} \tag{ 217 }$$

Thus far we have specified how to project the Wetterich equation onto flow equations for the couplings that parameterize the homogeneous part of the effective action given in equation (207), where the classical and quantum fields are set to constant values. In the following we review the derivation of flow equations for the frequency- and momentum-dependent couplings, i.e. the wave-function renormalization Z and the coefficient $\bar{K}$ of the Laplacian in equation (204). This requires us to consider non-constant values of the fields. Moreover, we work in a basis of real fields which we introduced already in section 2.4.6, ${{\bar{\phi}}_{\nu}}=\frac{1}{\sqrt{2}}\left({{\bar{\chi}}_{\nu,1}}+\text{i}{{\bar{\chi}}_{\nu,2}}\right)$ for $\nu =c,q$. Hence, the inverse propagator in this basis is given by the second variational derivative of the effective action with respect to the fields ${{\bar{\chi}}_{i}}$ (see equation (41); to ease the notation, we collect these fields in a vector $\bar{\chi}=\left({{\bar{\chi}}_{c,1}},{{\bar{\chi}}_{c,2}},{{\bar{\chi}}_{q,1}},{{\bar{\chi}}_{q,2}}\right)$; the components of this vector are labeled by $i=1,\ldots,4$), and the flow equation of the inverse propagator reads accordingly:

$$\begin{eqnarray}&&{{\partial}_{\ell}}{{\bar{P}}_{ij}}\left(\omega,\mathbf{q}\right)\delta \left(\omega -{{\omega}^{\prime}}\right)\delta \left(\mathbf{q}-{{\mathbf{q}}^{\prime}}\right)={{\left[\frac{{{\delta}^{2}}{{\partial}_{\ell}}{{\Gamma}_{\Lambda}}}{\delta {{{\bar{\chi}}}_{i}}\left(-\omega,-\mathbf{q}\right)\delta {{{\bar{\chi}}}_{j}}\left({{\omega}^{\prime}},{{\mathbf{q}}^{\prime}}\right)}\right]}_{\text{ss}}}.\end{eqnarray} \tag{ 218 }$$

In particular, the inverse retarded propagator is given by

$$\begin{eqnarray}{{\bar{P}}^{R}}\left(\omega,\mathbf{q}\right)=\left(\begin{array}{*{35}{l}} -\text{i}{{Z}_{I}}\omega -\bar{A}{{q}^{2}}-2\bar{\lambda}{{{\bar{\rho}}}_{0}} & \text{i}{{Z}_{R}}\omega -\bar{D}{{q}^{2}} \\ -\text{i}{{Z}_{R}}\omega +\bar{D}{{q}^{2}}+2\bar{\kappa}{{{\bar{\rho}}}_{0}} & -\text{i}{{Z}_{I}}\omega -\bar{A}{{q}^{2}} \end{array}\right).\end{eqnarray} \tag{ 219 }$$

Note that the Goldstone theorem (i.e. the existence of a zero eigenvalue of ${{\bar{P}}^{R}}\left(\omega =0,\mathbf{q}=0\right)$, see section 2.4.6) is preserved during the flow. At the transition, where ${{\bar{\rho}}_{0}}\to 0$, both branches of the excitation spectrum that is encoded in the zeros of $\det {{\bar{P}}^{R}}\left(\omega,\mathbf{q}\right)$ (see equation (77)) become gapless. As pointed out in section 2.5, in the FRG, the resulting infrared divergences are regularized by introducing an additional contribution $\Delta {{S}_{\Lambda}}$, given in equation (127), in the functional integral. In fact, by choosing the following optimized form of the cutoff function [278, 279] (which obviously satisfies the requirements stated in equations (130) and (131)),

$$\begin{eqnarray}&&{{R}_{\Lambda,\bar{K}}}\left({{q}^{2}}\right)=-\bar{K}\left({{\Lambda}^{2}}-{{q}^{2}}\right)\theta \left({{\Lambda}^{2}}-{{q}^{2}}\right),\end{eqnarray} \tag{ 220 }$$

in the regularized inverse propagator $\Gamma_{\Lambda}^{(2)}+{{R}_{\Lambda}}$ appearing in the Wetterich equation (132), the terms $\bar{A}{{q}^{2}}$ in equation (219) are replaced by

$$\begin{eqnarray}\bar{A}\left[{{q}^{2}}+\left({{\Lambda}^{2}}-{{q}^{2}}\right)\theta \left({{\Lambda}^{2}}-{{q}^{2}}\right)\right]=\left\{\begin{array}{*{35}{l}} \bar{A}{{\Lambda}^{2}} & \text{for}\quad {{q}^{2}}<{{\Lambda}^{2}}, \\ \bar{A}{{q}^{2}} & \text{for}\quad {{q}^{2}}\geqslant {{\Lambda}^{2}}, \end{array}\right. \end{eqnarray} \tag{ 221 }$$

(and there is an analogous replacement for the terms $\bar{D}{{q}^{2}}$). Hence, fluctuations with momenta below the cutoff scale $\Lambda$ acquire a mass $\sim {{\Lambda}^{2}}$, and the infrared divergences are lifted.

It remains to specify the flow equations for Z and $\bar{K}$. They can be obtained by choosing specific values of the indices i and j in the flow equation (218) for the inverse propagator, and by taking derivatives with respect to the frequency ω and the squared momentum q², respectively:

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{\partial}_{\ell}}Z & =-\frac{1}{2}{{\partial}_{\omega}}\text{tr}\left[\left({\mathbb 1}+{{\sigma}_{y}}\right){{\partial}_{\ell}}{{{\bar{P}}}^{R}}\left(\omega,\mathbf{q}\right)\right]{{|}_{\omega =0,\mathbf{q}=0}}, \end{array}\end{eqnarray} \tag{ 222 }$$

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{\partial}_{\ell}}\bar{K} & =-{{\partial}_{{{q}^{2}}}}\left({{\partial}_{\ell}}\bar{P}_{22}^{R}\left(\omega,\mathbf{q}\right)+\text{i}{{\partial}_{\ell}}\bar{P}_{12}^{R}\left(\omega,\mathbf{q}\right)\right){{|}_{\omega =0,\mathbf{q}=0}}. \end{array}\end{eqnarray} \tag{ 223 }$$

Comparison with equation (219) shows, that these are indeed correct projection prescriptions. Note, however, that there is some ambiguity in choosing these projection prescriptions: for example, $\bar{A}$ appears both in $\bar{P}_{11}^{R}$ and $\bar{P}_{22}^{R}$. Our choice extracts $\bar{K}$ corresponding to the Goldstone direction, and mixes Goldstone and gapped directions symmetrically in the projection on Z (see [181] for details). Finally, the flow equation for the renormalized coefficient $K=\bar{K}/Z$ is given by

$$\begin{eqnarray}&&{{\partial}_{\ell}}K={{\beta}_{K}}={{\eta}_{Z}}K+{{\partial}_{\ell}}\bar{K}/Z.\end{eqnarray} \tag{ 224 }$$

While equations (214)–(217) and (224) define a closed system of flow equations for the couplings ${{u}_{2}},{{u}_{3}},\gamma,{{\rho}_{0}},$ and K (as indicated above, the wave-function renormalization Z drops out of these equations), the explicit evaluation of the various projections is rather tedious. For details of this calculation we refer the reader to [181].

At the critical point of a continuous phase transition, the correlation length diverges, $\xi \to \infty$. Then, instead of exponential decay according to $\sim {{\text{e}}^{-r/\xi}},$ correlation and response functions depend on distance as power laws. This algebraic scaling behavior is reflected in the RG flow: indeed, the critical point corresponds to a scaling solution to the RG flow equations. In practice, finding a scaling solution is facilitated by absorbing the scaling factors $\sim {{\Lambda}^{\theta}}$ (with some exponent θ for each coupling) in new variables, which thus take constant values at the critical point. Hence, the latter corresponds to a fixed point of the flow equations for the rescaled couplings. Moreover, by means of a suitable choice of rescaled variables it is often possible to further reduce the number of flow equations, ending up with a minimal set of independent equations. For the present case, the flow can be specified in terms of just six real couplings.

At the beginning of this section, we introduced the quantity λ in equation (188) as a quantitative measure of the deviations from thermal equilibrium conditions. In a similar spirit, the strength of coherent relative to dissipative dynamics, which is encoded in the real and imaginary parts of the couplings in the microscopic Keldysh action (73), is measured by the ratios

$$\begin{eqnarray}&&\mathbf{r}=\left({{r}_{K}},{{r}_{{{u}_{2}}}},{{r}_{{{u}_{3}}}}\right)=\left(A/D,\lambda /\kappa,{{\lambda}_{3}}/{{\kappa}_{3}}\right).\end{eqnarray} \tag{ 225 }$$

The flow equations for these ratios can be obtained straightforwardly by taking the RG scale derivatives, e.g. ${{\partial}_{\ell}}{{r}_{K}}={{\partial}_{\ell}}A/D-A{{\partial}_{\ell}}D/{{D}^{2}}$, and expressing ${{\partial}_{\ell}}A$ as the real part of equation (224) etc. In addition to $\mathbf{r}$, we define another three scaling variables as

$$\begin{eqnarray}&&\mathbf{s}=\left(w,\tilde{\kappa},{{\tilde{\kappa}}_{3}}\right)=\left(\frac{2\kappa {{\rho}_{0}}}{{{\Lambda}^{2}}D},\frac{\gamma \kappa}{2\Lambda{{D}^{2}}},\frac{{{\gamma}^{2}}{{\kappa}_{3}}}{4{{D}^{3}}}\right).\end{eqnarray} \tag{ 226 }$$

The flow equations for the six dimensionless running couplings collected in $\mathbf{r}$ and $\mathbf{s}$ form a closed set. Besides the Gaussian fixed point at which the non-linear couplings vanish, these equations have a non-trivial fixed point corresponding to the driven–dissipative condensation transition at

$$\begin{eqnarray}&&{{\mathbf{r}}_{\ast}}=0,\quad {{\mathbf{s}}_{\ast}}\approx (0.475,5.308,51.383).\end{eqnarray} \tag{ 227 }$$

This fixed point is reached in the RG flow, when the parameters in the microscopic action are chosen such that the system is tuned precisely to the transition point. (Note that this point corresponds to the renormalized value of $w\propto {{\rho}_{0}}$ going to zero, and not the bare one.) What are the physical implications of this fixed point? First, the value ${{\mathbf{r}}_{\ast}}$ indicates, that the effective action at the fixed point is purely dissipative. As we have already mentioned at the beginning of section 4, for vanishing coherent dynamics (or, in the terminology of equilibrium dynamical models [1]: in the absence of reversible mode couplings), the driven–dissipative model reduces to the equilibrium model A. Thus, the values ${{\mathbf{s}}_{\ast}}$ are the same as in model A, and therefore the fixed point itself does not allow to distinguish whether the microscopic starting point of the RG flow was in or out of equilibrium. However, the non-equilibrium nature of the driven–dissipative condensate is witnessed in the RG flow towards this effective equilibrium fixed point. In the following we consider the universal regime of the RG flow, which is reached in the deep IR (i.e. for $\Lambda/{{\Lambda}_{0}}\ll 1$). In this regime, when the couplings are close to their values at the fixed point, the RG flow can be obtained from a linearization of the flow equations in $\delta \mathbf{s}=\mathbf{s}-{{\mathbf{s}}_{\ast}},\delta \mathbf{r}=\mathbf{r}$. The stability matrix governing the linearized flow takes block diagonal form,

$$\begin{eqnarray}&&{{\partial}_{\ell}}\left(\begin{array}{*{35}{l}} \delta \mathbf{r} \\ \delta \mathbf{s} \end{array}\right)=\left(\begin{array}{*{35}{l}} N & 0 \\ 0 & S \end{array}\right)\,\left(\begin{array}{*{35}{l}} \delta \mathbf{r} \\ \delta \mathbf{s} \end{array}\right),\end{eqnarray} \tag{ 228 }$$

with $3\times 3$ submatrices N and S. This block-diagonal structure indicates, that the flow of $\mathbf{r}$ and $\mathbf{s}$ decouples in the IR. Therefore, the flow of $\mathbf{s}$ close to the fixed point is the same as if we would have set $\mathbf{r}=0$ from the very beginning. In other words, not only the values of the couplings $\mathbf{s}$ at the fixed point, but also the critical exponents encoded in the flow of $\mathbf{s}$, which are the correlation length exponent ν, the anomalous dimension η, and the dynamical exponent z, see [181], are the same for both model A and driven–dissipative condensates. (Note, however, that in model F [1], which describes condensation with particle number conservation in equilibrium, the dynamical exponent takes a different value than in model A, where particle number conservation is absent.) This confirms the asymptotic thermalization of correlation functions mentioned in section 4.2 (i). The values of the critical exponents we obtain from the truncation (204) are

$$\begin{eqnarray}&&\nu \approx 0.716,\quad \eta \approx 0.039,\quad z\approx 2.121,\end{eqnarray} \tag{ 229 }$$

and agree reasonably well with results from more sophisticated calculations of ν and η in the context of the static equilibrium problem [280].

All the information on the universal properties of the driven–dissipative transition, which are the same as in the equilibrium model A, are encoded in ${{\mathbf{s}}_{\ast}}$ and the block S of the stability matrix in equation (228). The non-equilibrium nature of the microscopic model, on the other hand, is betrayed by the block N, which describes the flow of $\delta \mathbf{r}$. This block has three positive eigenvalues,

$$\begin{eqnarray}&&{{n}_{1}}\approx 0.101,\quad {{n}_{2}}\approx 0.143,\quad {{n}_{3}}\approx 1.728,\end{eqnarray} \tag{ 230 }$$

indicating that the ratios $\mathbf{r}$ are attracted to the fixed point value ${{\mathbf{r}}_{\ast}}=0$. The general solution to the linearized flow equation for $\mathbf{r}$ reads

$$\begin{eqnarray}&&\mathbf{r}=\underset{i=1}{\overset{3}{\sum}}\,{{\mathbf{u}}_{i}}{{c}_{i}},\end{eqnarray} \tag{ 231 }$$

where ${{\mathbf{u}}_{i}}$ are the eigenvectors of N associated with the eigenvalues n_i in equation (230). The coefficients c_i, which are referred to as scaling fields [178], take the scaling form ${{c}_{i}}\sim {{\text{e}}^{{{n}_{i}}\ell}}\sim {{\Lambda}^{{{n}_{i}}}}$. Hence, for $\Lambda\to 0$, the dominant contribution to $\mathbf{r}$ is given by $\mathbf{r}\sim {{\mathbf{u}}_{1}}{{\Lambda}^{{{n}_{1}}}}={{\mathbf{u}}_{1}}{{\Lambda}^{-{{\eta}_{r}}}}$, where we identified the drive exponent

$$\begin{eqnarray}&&{{\eta}_{r}}=-{{n}_{1}}\approx -0.101.\end{eqnarray} \tag{ 232 }$$

As anticipated in section 4.2 (ii), this exponent governs the universal fade-out of coherent dynamics $\propto \mathbf{r}$ in the driven–dissipative system. Note that the existence of three distinct eigenvalues (230) is due to the non-equilibrium character of the microscopic model. Indeed, in model A with reversible mode couplings, the equilibrium symmetry discussed in section 2.4.1 allows of only one ratio $r={{r}_{K}}={{r}_{{{u}_{2}}}}={{r}_{{{u}_{3}}}}$ (see equation (101) and figure 6). Then, the block N of the stability matrix in equation (228) has only one single entry, which is given by the 'middle' eigenvalue n₂ in equation (230). This shows that also in the equilibrium setting the dynamics becomes purely dissipative at the largest scales [268], however, the value of the critical exponent that governs universal decoherence is different. As pointed out above, this is due to the absence of the equilibrium symmetry in the driven–dissipative case. The fact that the different values can be traced back to a difference in symmetry supports the strength of the result, as different symmetries are known to give rise to quantitatively different critical behavior [86].

Decoherence at large scales has clear physical signatures which facilitate probing the drive exponent in experiments; to wit, it implies that low-momentum excitations are diffusing rather than propagating (note that this is also predicted by mean-field theory, see the discussion below equation (77)). A careful analysis in the scaling regime reveals that the effective dispersion relation of single-particle excitations close to criticality takes the form [181]

$$\begin{eqnarray}&&\omega \sim {{A}_{0}}{{q}^{z-{{\eta}_{r}}}}-\text{i}{{D}_{0}}{{q}^{z}}\sim {{A}_{0}}{{q}^{2.223}}-\text{i}{{D}_{0}}{{q}^{2.121}},\end{eqnarray} \tag{ 233 }$$

where the diffusive contribution is supplemented by a subdominant (by the small difference of ${{\eta}_{r}}$ in the exponent) coherent part. By definition, equation (233) is the location of the pole of the retarded Green's function, and hence the coherent and diffusive parts encode respectively the position and width of the peak of the spectral function defined in equation (61). The latter is probed, e.g. in angle-resolved spectroscopy in exciton-polariton systems [281] or radio-frequency spectroscopy in ultracold atoms [282]. However, the small difference in the scaling of the position and width of the peak predicted by equation (233) poses a challenge to its experimental observation.

For sufficiently small momenta $q\ll 1$, the kinetic equation simplifies considerably. In this case, the second line of equation (284) can be discarded completely, since the integral is performed over a very small momentum interval 0  <  p  <  q. On the other hand, for the integral in the third line of equation (284), all terms $(q+p)\approx p$ for the dominant part of the integral. Therefore

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{\partial}_{t}}{{n}_{q,t}}=\,|q|\left(\frac{\gamma}{2\pi K}+{{\mathcal{J}}_{t}}\right), \end{array}\end{eqnarray} \tag{ 285 }$$

where the time-dependent integral ${{\mathcal{J}}_{t}}$ is q-independent and reads

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{\mathcal{J}}_{t}}=2v_{0}^{2}{\int}_{0<p}\frac{{{p}^{2}}{{n}_{p}}\left({{n}_{p}}+1\right)}{\sigma _{p}^{R}}. \end{array}\end{eqnarray} \tag{ 286 }$$

As a consequence, the resulting change of the phonon density for small momenta $|q|$ is linearly proportional to $|q|$, i.e.

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{n}_{q,t}}={{n}_{q,0}}\,+\,|q|\left(\frac{\gamma t}{2\pi K}+{{{\int}^{}}_{0<{{t}^{\prime}}<t}}{{\mathcal{J}}_{{{t}^{\prime}}}}\right). \end{array}\end{eqnarray} \tag{ 287 }$$

The momentum regime $0<|q|<{{q}_{2}}$, defined as the regime where the non-linear contributions in $|q|$ are negligible, depends on the actual value of the phonon density n_q, t, and has to be determined for each specific scenario individually. However, the existence of q₂ is guaranteed by the above general arguments, and the fact that ${{\partial}_{t}}{{n}_{q=0,t}}=0$ for all times t. The latter is a consequence of the U(1) symmetry of the present setup and the global particle number conservation, as discussed in [44].

The fact that the present system is described by a U(1)-invariant, massless field theory is reflected by the absence of a scale in the self-energy equation (283). One important consequence is that $\sigma _{q,t}^{R}\overset{q\to 0}{{\to}}\,0$ generically, i.e. the generation of a mass gap is forbidden by symmetry. A further consequence of equation (283) is that, whenever the term in brackets obeys a scaling law,

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \frac{{{\partial}_{t}}{{n}_{p}}}{\sigma _{p}^{R}}+2{{n}_{p}}+1\sim {{\gamma}_{n}}|\,p{{|}^{{{\eta}_{n}}}}, \end{array}\end{eqnarray} \tag{ 288 }$$

the solution for the self-energy will be a scaling function $\sigma _{q}^{R}={{\gamma}_{R}}|\,p{{|}^{{{\eta}_{R}}}}$ as well. Inserting this scaling ansatz in equation (283) yields

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{\gamma}_{R}}|\,p{{|}^{{{\eta}_{R}}}}\,=\frac{v_{0}^{2}{{\gamma}_{n}}}{{{\gamma}_{R}}}|\,p{{|}^{4+{{\eta}_{n}}-{{\eta}_{R}}}}{{\mathcal{I}}_{{{\eta}_{n}},{{\eta}_{R}}}} \end{array}\end{eqnarray} \tag{ 289 }$$

and identifies

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{\eta}_{R}}=2+\frac{{{\eta}_{n}}}{2}\,\text{and}\,{{\gamma}_{R}}={{v}_{0}}\sqrt{{{\gamma}_{n}}{{\mathcal{I}}_{{{\eta}_{n}},{{\eta}_{R}}}}}. \end{array}\end{eqnarray} \tag{ 290 }$$

The dimensionless integral ${{\mathcal{I}}_{{{\eta}_{n}},{{\eta}_{R}}}}$ is defined as

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{\mathcal{I}}_{{{\eta}_{n}},{{\eta}_{R}}}}={{{\int}^{}}_{0<x}}{{x}^{{{\eta}_{n}}}}\left(\frac{x(1-x)}{{{x}^{{{\eta}_{R}}}}\,+\,|1-x{{|}^{{{\eta}_{R}}}}}+\frac{x(1+x)}{{{x}^{{{\eta}_{R}}}}+{{(1+x)}^{{{\eta}_{R}}}}}\right), \end{array}\end{eqnarray} \tag{ 291 }$$

and depends only on the exponents ${{\eta}_{R}},{{\eta}_{n}}$. This self-energy integral is dominated by contributions around x  =  1 and converges for all physically reasonable phonon densities. In this sense, the scaling behavior is universal, i.e. robust against the influence from high energy modes [43, 44].

The results of a numerical simulation of the phonon densities are shown in figure 13. One can clearly identify the crossover from an interaction dominated thermal regime with ${{n}_{q,t}}=\frac{T(t)}{u|q|}$ at large momenta to a heating dominated regime with ${{n}_{q,t}}\sim |q|$. The crossover momentum between the two regimes represents the inverse thermal length ${{\left({{x}_{\text{th}}}(t)\right)}^{-1}}\sim {{t}^{-4/5}}$.

Zoom In Zoom Out Reset image size

In the large momentum regime, the dominant contribution to the kinetic equation is the collision term, which establishes an approximate detailed balance between phonon absorption and emission in the presence of the heating. This is expressed by the evolution of the phonon density towards a Bose–Einstein distribution, which is the fixed point of the collision term alone. In this regime ${{n}_{q,t}}={{n}_{\text{B}}}\left(u|q|,T(t)\right)\approx \frac{T(t)}{u|q|}$ with very good agreement and the only indicator of the permanent heating is the continuously increasing temperature $T(t)\sim t$. This is contrasted by the evolution of the phonon density in the low momentum regime, which is dominated by strong phonon production. In this regime, the scattering of high-momentum phonons into low-momentum modes enhances the effect of the heating and leads, due to the structure of the vertex, to an increase of the linear production rate according to equation (285). The pinning of the phonon occupation of the q  =  0 mode is an exact result for the underlying U(1)-symmetric, i.e. particle number conserving dynamics. It can be shown that the phonon number fluctuations in the zero momentum mode are proportional to the fluctuations of the global particle number in the system, see [43], and consequently they are integrals of motion for a U(1)-symmetry preserving dynamics. The pinning effect for the low momentum distribution at ${{n}_{q=0,t}}={{n}_{q=0,0}}$ leads to a very slow, sub-ballistic thermalization of the low momentum regime, since the formation of the typical, thermal Rayleigh–Jeans divergence ${{n}_{q,t}}\sim 1/|q|$ is only achieved by the scaling of the inverse thermal length ${{\left({{x}_{\text{th}}}(t)\right)}^{-1}}$ to zero, instead of a direct filling of these modes.

Away from the crossover scale $q\ne x_{\text{th}}^{-1}$, the phonon density exhibits scaling behavior and can be written as

$$\begin{eqnarray}{{n}_{q,t}}=\left\{\begin{array}{c} c(t)|q| & \text{for}~|q|\ll x_{\text{th}}^{-1} \\ \frac{T(t)}{u|q|} & \text{for}~|q|\gg x_{\text{th}}^{-1} \end{array},\right. \end{eqnarray} \tag{ 292 }$$

where the functions c(t), T(t) have to be determined numerically. In these regimes, according to section 5.3.2, one finds scaling behavior of the phonon self-energy as well. The corresponding scaling exponent is determined by

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{f}_{q,t}}\equiv \frac{{{\partial}_{t}}{{n}_{q,t}}}{\sigma _{q,t}^{R}}+2{{n}_{q,t}}+1\sim |q{{|}^{{{\eta}_{n}}}}. \end{array}\end{eqnarray} \tag{ 293 }$$

In the heating dominated regime, ${{n}_{q,t}}=c(t)|q|$ and therefore

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{f}_{q,t}}=\frac{{{c}^{\prime}}(t)}{{{\gamma}_{R}}}|q{{|}^{1-{{\eta}_{R}}}}\,+\,c(t)|q|\,+\,1\overset{{{\eta}_{R}}>1}{{\to}}\,\frac{{{c}^{\prime}}(t)}{{{\gamma}_{R}}}|q{{|}^{1-{{\eta}_{R}}}} \end{array}\end{eqnarray} \tag{ 294 }$$

for small momenta, since ${{\eta}_{R}}>1$ is guaranteed by the subleading nature of the interactions. This directly implies ${{\eta}_{n}}=1-{{\eta}_{R}}$ and a super-diffusive exponent ${{\eta}_{R}}=\frac{5}{3}$ in this regime, which lacks any equilibrium counterpart.

On the other hand, in the large momentum regime f_q, t is obviously dominated by the term proportional to n_q, t and ${{\eta}_{n}}=-1$, which leads to the known thermal equilibrium exponent ${{\eta}_{R}}=\frac{3}{2}$ [320] and a thermal scaling behavior for large momenta not only in the distribution function but also in the self-energy. The crossover between the two scaling regimes can be estimated by equating the two relevant terms, which tend to dominate f_{q, t} in the corresponding regimes. This yields

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \frac{{{\partial}_{t}}{{n}_{{{q}_{\text{th}}}}}}{2{{n}_{{{q}_{\text{th}}}}}+1}=\sigma _{{{q}_{\text{th}}}}^{R}, \end{array}\end{eqnarray} \tag{ 295 }$$