# Nonequilibrium Dyson equations for strongly coupled light and matter:
spin glass formation in multi-mode cavity QED

Hossein Hosseinabadi  Institut für Physik, Johannes Gutenberg-Universität
Mainz, 55099 Mainz, Germany  Darrick E. Chang  ICFO—Institut de Ciències
Fotòniques, The Barcelona Institute of Science and Technology, 08860
Castelldefels, Spain  ICREA—Institució Catalana de Recerca i Estudis Avançats,
08015 Barcelona, Spain  Jamir Marino  Institut für Physik, Johannes Gutenberg-
Universität Mainz, 55099 Mainz, Germany

###### Abstract

Light-matter interfaces have now entered a new stage marked by the ability to
engineer quantum correlated states under driven-dissipative conditions. To
propel this new generation of experiments, we are confronted with the need to
model non-unitary many-body dynamics in strongly coupled regimes, by
transcending traditional approaches in quantum optics. In this work, we
contribute to this program by adapting a functional integral technique,
conventionally employed in high-energy physics, in order to derive
nonequilibrium Dyson equations for interacting light-matter systems. Our
approach is grounded in constructing ’two-particle irreducible’ (2PI)
effective actions, which provide a non-perturbative and conserving framework
for describing quantum evolution at a polynomial cost in time. One of the aims
of the article is to offer a pedagogical introduction designed to bridge
readers from diverse scientific communities, including those in quantum
optics, condensed matter, and high-energy physics. We apply our method to
complement the analysis of spin glass formation in the context of frustrated
multi-mode cavity quantum electrodynamics, initiated in our accompanying work
[H. Hosseinabadi, D. Chang, J. Marino, arXiv:2311.05682]. Finally, we outline
the capability of the technique to describe other near-term platforms in many-
body quantum optics, and its potential to make predictions for this new class
of experiments.

##  I Introduction

The field of quantum simulation has recently undergone a transformation,
evolving into a new realm of research where the worlds of condensed matter
physics and quantum optics merge. Today, an increasing array of platforms
hosting many-body systems can accommodate both unitary and dissipative
dynamics in a controlled fashion [1, 2, 3, 4, 5, 6]. This circumstance paves
the way for the exploration of phases of matter and strongly correlated
behavior that have no counterparts in either thermodynamic equilibrium or
isolated non-equilibrium conditions.

In contrast to conventional solid-state systems, driven-dissipative condensed
matter (a.k.a many-body quantum optics) systems exhibit a host of innovative
characteristics. In the former, dissipation poses the primary challenge to
quantum coherence, while in the latter, dissipation is at times intentionally
harnessed or even engineered to drive the system into entangled states [7, 8,
9].

Traditional solid-state physics focuses on understanding equilibrium phases of
matter that result from the interplay of interactions and thermal fluctuations
[10, 11, 12, 13]. In the realm of driven-dissipative condensed matter,
instead, emergent behavior can occur in dynamics or in the non-equilibrium
steady state where the system settles. These novel responses usually arise
from the intricate interplay of classical and quantum noise within the
strongly coupled limit of a many-particle system. Finally, the nature of
interactions fostering strong correlations takes on a fundamentally different
character in the domain of quantum many-body optics. Here, short-range
interactions originating from atoms or molecules coexist with long-range or
even all-to-all interactions mediated by light, leading to fundamentally
distinct cooperative mechanisms.

These distinctive features find natural occurrence in cavity quantum
electrodynamics (cQED) experiments, which have occupied a central position in
the field of driven-open quantum simulators for over a decade. In cavity QED,
atoms and photons experience couplings as a consequence of their confinement
within high-finesse optical cavities [14]. This confinement enables repeated
light-matter scattering in nearly isolated conditions, achieving effectively
strong interactions. The advent of ultra-cold quantum gases has made it
possible to couple Bose-Einstein condensates to one or more photonic modes
within such cavities, creating regimes of cooperative dynamics that intertwine
matter waves and light [15]. Leveraging the flexibility, tunability, and
engineering capabilities inherent in quantum simulators, these platforms have
become ideal environments for realizing non-equilibrium phases of matter under
driven-open conditions. These include phenomena such as dynamical
ferromagnetism, super-solids, topological pumps, time crystals, spin networks,
and many others [15]. Beyond fundamental research, applications extend to
quantum information, where cavity QED currently holds the world record for
spin squeezing [16, 17], an entangled state which surpasses classical limits
for metrology and sensing.

In the majority of these experiments, the dynamical behavior of the system can
be reduced to a few degrees of freedom: typically, these are photon amplitudes
and numbers in conjunction with the dynamics of a collective spin that
encapsulates the motion of the entire assembly of atoms [18, 19, 20, 21, 22,
23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34]. This usually results from the
all-to-all nature of photon-mediated interactions in the cavity, which allows
for an effective mean-field description. Quantum fluctuations are sub-leading
in the number of atoms when the dynamics of such macroscopic degrees of
freedom are considered: only weak quantum noise influences the correlation
between light and matter, and the system’s behavior can be effectively
characterized through semi-classical descriptions [23].

This situation presents significant advantages when it comes to modeling
cavity QED experiments, and it has certainly served as the foundation for
constant advancements in joint theoretical and experimental efforts over the
past fifteen years. Simultaneously, it naturally prompts the exploration of
conditions where the strong correlations inherent in these systems become
dominant. In this context, we are witnessing the emergence of a novel category
of experiments, where the true many-body nature of the platform emerges and it
defies a description using only a few macroscopic degrees of freedom. These
experiments encompass a variety of setups, such as Rydberg tweezer arrays
integrated into optical cavities [35], atomic ensembles with adjustable
loading capacities [36, 37], and the fermionic variant of traditional cavity
QED experiments [38].

All these experiments necessitate the integration of methods traditionally
employed for addressing strongly correlated problems with the extra
complication that detailed balance is broken and dynamics are non-unitary.
Accessing the long-term evolution of open many-body systems subject to
external (coherent or incoherent) drives as well as strong interactions, is
one of the most challenging computational frontiers. Nevertheless, it has now
become increasingly important in order to steer this new generation of cavity
QED platforms.

Indeed, as numerous established platforms push the boundaries of the NISQ
(Noisy Intermediate-Scale Quantum) era [39], a cavity QED platform endowed
with strong correlations presents an intriguing opportunity. It would hold the
potential for innovative strategies in quantum processing, leveraging the
combined advantages offered by cooperative behavior [40], the interplay of
long and short-range interactions [35, 41, 42], the manipulation of
controllable quantum fluctuations [43], including the potential to manipulate
decoherence channels [44].

As of today, state of art methods to model these experiments would encompass a
number of options. Refined versions of semi-classical techniques, such as
discrete phase space representations of the Hilbert space, have been tailored
for addressing the driven-dissipative dynamics of spin and bosonic systems
[45, 46, 47, 48, 49, 50]. Tensor network ansätze, which retain the most
informaive correlations to describe dynamics, have been used to study phase
transitions in driven open systems [51, 52, 53, 54, 55, 56, 57, 58]. At the
same time, cluster mean field methods [59] or self-consistent Gaussian
approximations [60] remain a straightforward, and, in some circumstances,
competitive way, to qualitatively capture non-unitary dynamics.

Figure 1: Schematics of the model considered in this work: $N$ clusters (tan
circles) of $N_{s}$ two-level atoms (spins) interact with $M$ cavity modes
(green lines) via static disordered couplings $g_{\alpha,i}$. Cavity modes
mediate effective all-to-all interactions modulated among clusters. The
effective interaction is modulated by static disorder and generates
frustration and glassy dynamics. Each cluster can be taken as a large spin $S$
with amplitude $N_{s}/2$. Cavity modes induce ferromagnetic (FM) interactions
between spins within the same cluster.

In this paper, we introduce a novel method to tackle the dynamics of strongly
interacting light-matter systems in the many-body limit [61, 62]. It is an
adaptation of the nonequilibrium Dyson equations (DE) derived from functional
integrals in high energy physics and cosmology [63, 64, 65, 66, 67, 68, 69,
70, 71, 72, 73, 74, 75], and previously used also in the context of AMO and
condensed matter [76, 77, 78, 79, 80, 81, 82, 83, 84]. Our method involves
deriving an effective quantum action for the system as a function of its
$n$-point correlation functions, from which DE descend following the least
action principle [61, 62]. Here, we concentrate on deriving dynamics from what
is referred to as the two-particle irreducible (2PI) effective action [63, 61,
62]. The term “two-particle irreducible” pertains to the Feynman diagrams
contributing to the kernel of the DE, which cannot be separated into two parts
by cutting at most two of their lines. Such action gives exact equations of
motion of two-point functions, and with some ingenuity also the dynamics of
higher point correlation functions [85, 86]. It is via a set of controlled
non-perturbative approximations (large $N$ limit, dilute expansion, loop
expansion) [61, 62, 79, 82] that the equations are actually solved
numerically.

DE derived from 2PI effective actions offer numerous advantages. They descend
from a variational principle, and therefore they give rise to self-consistent
dynamics of two-point functions free from secular effects [62, 87]. Since
approximations are directly performed at the level of the action, the
resulting DE are _conserving_ , in the sense they cannot spoil the conserved
quantities of the original model. They also can seamlessly incorporate both
coherent and dissipative dynamics, as functional integral methods do not
markedly differentiate between the two [88]. They do not suffer from
limitations when degrees of freedom with unbounded Hilbert space, like photons
or phonons, are included in dynamics. The dynamics governed by DE have
polynomial time costs in system size and they can be run on a personal
computer. The price to pay is formulating educated guesses on the physics of
the problem, which are crucial for selecting the proper approximation scheme.
Furthermore, these DE are known to semi-quantitatively reproduce dynamics, and
they are therefore excellent for elucidating the mechanisms at work in a given
problem of interest, but less suited to fit with accuracy experimental curves,
for instance. With these caveats, the method is highly flexible and applicable
virtually to any driven-dissipative many particle system, made of fermions,
bosons or spins, as we also expand in the conclusions of this manuscript.

In this work, we initiate our research program by applying this method to
examine the dynamics of strongly correlated light-matter systems within multi-
mode cavity QED. Our model is inspired by the experimental setup presented in
Refs. [37, 89, 90, 91] which involves several photonic modes connecting nodes
(Fig. 1). Each of these nodes houses atomic ensembles with adjustable loading
capacities. By manipulating the number of atoms in each node (ranging from a
few to thousands in the experiment) by using optical tweezers, one can
introduce tunable quantum fluctuations in the platform. These fluctuations
enable the exploration of system dynamics, from strongly correlated to semi-
classical regimes. Importantly, this flexibility is not unique to this
platform [36, 92] and represents a promising starting point for delving into
many-body cavity QED beyond the domain of collective dynamic responses [19,
20]. During the final stages of the current work, the experiment detailed in
Ref. [93] verified the existence of a spin glass (SG) phase within the quantum
gas microscope platform of Refs. [37, 89, 90, 91] for the first time through
direct measurement of the configuration of spins in the system. Notably, the
observed spin-glass phase exhibits replica symmetry breaking [94], a
distinctive and intriguing characteristic of certain models of SG, whose
experimental verification is well-known to be a challenging task. In this
paper, we complement the findings presented in our accompanying work on
dynamical spin glass formation [95] by investigating the full spectrum of non-
equilibrium phases and crossovers inherent to the experiment in [37, 89, 90,
91]. Furthermore, we offer a pedagogical introduction to the formalism of DE
derived from effective actions. Our objective is to bridge the domains of AMO
and the many-body community working at the interface of condensed matter and
field theory.

###  I.1 Outline of the paper

The goal of this work is to understand the far from equilibrium dynamics of SG
phases in frustrated cavity QED with strong disorder, where fluctuations
cannot be omitted and mean field treatments are not applicable. We remark that
our approach is distinct from those of Refs. [96, 97] which by construction,
are suitable only for the universal SG behavior at steady state, in two key
aspects. First, the formalism developed here is applicable in far from
equilibrium situations such as quench dynamics, while keeping track of the
quantum nature of spins. Second, the platform of Refs. [37, 89, 90, 91]
naturally includes extra FM interactions which are unimportant in the steady
state of the system [43], but as shown in this work, can qualitatively modify
conventional SG dynamics away from equilibrium. We now briefly outline our key
results which expand upon the results of our accompanying work [95].

_Introduction to the model_ – We commence the paper by introducing the model
in Section II. We will briefly review previous works on the behavior of the
model in different regimes of parameters at the steady state, and will
motivate using our approach to treat its far from equilibrium dynamics.

_Introduction to the method_ – We provide a comprehensive introduction to
effective action methods and non-equilibrium field theory in Section III. This
introduction is anchored in a pedagogical example drawn from basic quantum
mechanics – the dynamics of the anharmonic quantum oscillator. This serves as
a foundation for the detailed derivation of our formalism in the context of
frustrated light-matter interactions in cavity QED, covered in Section IV,
where we develop a versatile approach to address real-time dynamics in a model
for disordered cavity QED given in Fig. 1. To enhance accessibility, this
section and its accompanying appendices are crafted to be reproducible from
scratch by the interested reader. Readers with experience on 2PI methods can
gloss over these two sections.

Figure 2: (a) Cartoon of the dependence of dynamics on spin size $N_{s}$ for
initial states with finite magnetization. Small spins relax fast while for
large spins, after an initial quick decay of FM order, the system enters a
prethermal state with slowly decaying magnetization. SG order is built up
faster for smaller spins. (b) Spin glass order versus cavity detuning
$\omega_{c}$, showing a resonance close to atomic detuning $\Delta$.

_Magnetization dynamics_ – In Section V.1, we show that the mean-field (MF)
approximation, given by the leading order contribution in our approach,
predicts a paramagnetic (PM) to ferromagnetic (FM) phase transition of
effective spin degrees of freedom in the universality class of infinite range
Ising model, but it completely omits the effect of frustrated interactions
generated by static disorder in the system. FM interactions are mediated by
virtual photon exchange processes among spins within the same cluster, and in
contrast to inter-cluster couplings, are not frustrated. In Section V.2 we
demonstrate that, upon the non-perturbative incorporation of disorder and
fluctuations, our approach provides a dramatic improvement of MF results. As
the focus of this work is the far from equilibrium dynamics of this system
after interaction quenches, we first look at the dynamics of simple spin
observables, such as global magnetization. We find that if the system is
initiated in a symmetry broken state with a finite total magnetization, the
relaxation of magnetic order substantially depends on the coupling strength
and the size of atomic ensembles $N_{s}$ or equivalently, the amplitude of
large spins $\hat{S_{i}}$ per each cluster (red curves in Fig. 2a). For weak
couplings (Section V.2.1), global magnetization displays paramagnetic
oscillations which are weakly damped due to the dephasing generated by static
disorder. Upon increasing the coupling (Section V.2.2), spin relaxation
changes from underdamped dynamics to overdamped dynamics without oscillations.
In Section V.2.3 we address the effect of ensemble size $N_{s}$ on
magnetization dynamics. We show that in the overdamped regime and for small
$N_{s}$, global magnetization $\langle\hat{S}^{x}\rangle$ relaxes quickly to
zero while for large $N_{s}$, it experiences an initial quick collapse to a
finite value, and then stays in a transient prethermal state with a slow
spiral decay of the magnetization vector along the axis of temporary FM order.

_Dynamics of SG_ – To probe into the nature of the transition in the system as
magnetization dynamics change from underdamped to overdamped, we consider more
complex and richer spin observables in Section V.3. Particularly, and with the
expectation of exploring SG order in the system, we consider the dynamics of
two different order parameters for SG in Sections V.3.2 and V.3.3. First, we
consider the disorder average of the square of local magnetization of each
spin ensemble, a finite quantity in the glass phase, which measures how much
the frozen spin state is protected from fluctuations. We extend the formalism
of Keldysh field theory in order to access the real-time dynamics of this
order parameter which is done for the first time in this work, although there
is a technical similarity between our approach and the one introduced by Ref.
[98] to study measurement-induced phase transitions. We confirm that the
change in spin relaxation profile from underdamped to overdamped dynamics,
indeed, happens when a PM to SG phase transition takes place in the system. In
PM phase, the SG order parameter relaxes to zero after experiencing small
temporary fluctuations. For quenches into SG phase, on the other hand, the
order parameter relaxes to a finite value. Second, we consider the temporal
correlations of spins over long times, conventionally known as the Edwards-
Anderson (EA) [99, 100] order parameter, corresponding to the overlap of two
snapshots of the system taken at long time intervals. We confirm that EA order
parameter changes from zero in the PM phase, when the system loses its memory
quickly, to a finite value in the SG phase, proving that the system is
glassified. We proceed to show that spin fluctuations violate the fluctuation-
dissipation theorem [101, 87], a phenomenon conjectured to be closely related
to replica symmetry breaking in SG systems [102].

_Effect of spin size_ – Section V.3.4 addresses the effect of ensemble size
$N_{s}$ and initial conditions on the formation of SG phase. It is shown that
the rate at which SG order builds up is slower for initial states with larger
global magnetization, and is fastest if the system is initialized in the
ground-state of non-interacting spins. The dependence on the initial state is
found to be less significant for small $N_{s}$ and more pronounced for large
ensemble sizes (green curves in Fig. 2a). The deceleration of SG formation for
large clusters can be attributed to the long-lasting transient FM order
mentioned above, which is protected by non-frustrated intra-cluster
interactions, and acts as a competing order with SG [103, 104, 105, 106, 107]
over intermediate timescales. For initial states with small net magnetization,
the rate of SG formation is fastest and is not affected by cluster size.
Notice that accessing the dynamics of atomic ensembles with tunable loading
capacities is crucial for exploring the full spectrum of dynamical responses,
not only in the experiment discussed in this work, but also in other cavity
QED platforms [36]. This is in general a challenging task for conventional
methods in quantum optics, while in our field theory approach this can be
straightforwardly included since $N_{s}$ is natural parameter which appears
naturally in our diagrammatic expansion.

_Effect of resonant photons_ – In Section V.3.5, we consider the effect of
photon frequency on the glass phase by looking at SG order parameter in a wide
range of photon frequencies from fast photons to the resonance limit, where
atomic and cavity detunings are close, and below. As shown schematically in
Fig. 2b, SG order is peaked close to the resonance, while it saturates in the
adiabatic limit. At frequencies below the resonance, SG order is dramatically
suppressed, resembling the suppression of various types of order by low
frequency lattice distortions (phonons) in solid state physics. We also
address briefly the spectrum of low-lying excitations in the SG phase, showing
a continuum of sub-Ohmic modes at small energies.

The capability to include in the same set of dynamical equations variable
ranges of coupling, tunable values of $N_{s}$ and active photons, is one of
the key merits of our approach. It allows us to solve for the dynamics of the
full platform without the need to invoke large energy scale separations,
effective descriptions suited only to atomic or photonic degrees of freedom,
or to treat distinctly the quantum and semi-classical regimes. In this regard,
the method has a degree of flexibility that appears promising to treat other
strongly correlated driven-dissipative systems, as we discuss in the
concluding Section VI.

##  II The Model

The experiment in [89, 37, 90, 91] can be modeled by a system of $N$ clusters,
each one containing $N_{s}$ two-level atoms encoded by the spin-1/2 operators
$\sigma_{i\lambda}$, with cluster $1\leq i\leq N$ and atom indices
$1\leq\lambda\leq N_{s}$ , as shown in Fig. 1. The couplings between the atoms
and the $M$ photonic modes of the cavity are spatial-dependent and
uncorrelated from each other, which justifiemodellings their modeling via
random spin-boson couplings [96, 97]. Starting from the same initial state for
all spins, each cluster is equivalent to a single spin
$S_{i}=\sum_{\lambda}\sigma_{i\lambda}/2$ with amplitude $S=N_{s}/2$. The
parameter $S$ can be tuned by loading few or several atoms in each cluster,
and it dictates the strength of quantum fluctuations. For instance, at large
$S$ each cluster would be effectively described by a classical angular
momentum, since its quantum noise would scale down as $1/S$ [19, 108, 109]. A
minimal model for the system is given by the random Dicke model whose
evolution is governed by
$\partial_{t}\rho=-i\commutator{H}{\rho}+\sum_{\alpha=1}^{M}\mathcal{D}\quantity[a_{\alpha}]\rho$,
where

$H=\frac{\Delta}{2}\sum_{i,\lambda}\sigma^{z}_{i\lambda}+\sum_{\alpha}\omega_{\alpha}a^{\dagger}_{\alpha}a_{\alpha}\\\
+\frac{1}{\sqrt{(M+N)N_{s}}}\sum_{\begin{subarray}{c}i,\lambda,\alpha\end{subarray}}g_{\alpha
i}\quantity(a_{\alpha}+a_{\alpha}^{\dagger})\sigma^{x}_{i\lambda},$ (1)

and

$\mathcal{D}\quantity[a_{\alpha}]\rho=\kappa_{\alpha}\quantity(2a_{\alpha}\rho
a_{\alpha}^{\dagger}-\anticommutator{a_{\alpha}^{\dagger}a_{\alpha}}{\rho}).$
(2)

We assume that cavity modes are nearly degenerate such that
$\omega_{\alpha}=\omega_{c}$ and $\kappa_{\alpha}=\kappa$. The couplings
$g_{\alpha i}$ are assumed to be random and chosen from a Gaussian
distribution:

$\overline{g_{\alpha i}}=0,\quad\overline{g_{\alpha i}g_{\beta
j}}=\delta_{\alpha\beta}\delta_{ij}g^{2}.$ (3)

Couplings for spins in the same cluster are similar as we assume that the
spatial size of each cluster is smaller than the wavelength of cavity modes.
The scaling of the interaction term guarantees that the total energy is
extensive in system size [110]. For $\omega_{c}\gg g,\Delta$ photons are the
fastest degree of freedom in the problem; they quickly relax to stationary
value and approximately, they mediate instantaneous interactions among spins.
Assuming $\kappa=0$, photons can be adiabatically eliminated [111, 112, 113,
114, 43] and the model in (1) is mapped to

$\displaystyle H_{\mathrm{eff}}$
$\displaystyle=\sum_{i}H_{i}+H_{\mathrm{int}},$ (4) $\displaystyle H_{i}$
$\displaystyle\equiv\Delta
S^{z}_{i}-\frac{4}{(N+M)N_{s}\omega_{c}}\quantity(\sum_{\alpha}g_{\alpha
i}^{2})\quantity(S^{x}_{i})^{2},$ (5) $\displaystyle H_{\mathrm{int}}$
$\displaystyle\equiv-\frac{4}{(N+M)N_{s}\omega_{c}}\sum_{i\neq
j}\sum_{\alpha}g_{\alpha i}g_{\alpha j}S^{x}_{i}S^{x}_{j},$ (6)

where $S^{x}_{i}=\sum_{\lambda}\sigma^{x}_{i\lambda}/2$ is the total spin
operator for each cluster. For $N_{s}>1$, each cluster in Eq. (5) is an
infinite range quantum Ising model also known as the Lipkin-Meshkov-Glick
[115] (LMG) model with Hamiltonian

$H_{\mathrm{LMG}}=\Delta S^{z}-\frac{J}{N_{s}}(S^{x})^{2}.$ (7)

$H_{\mathrm{LMG}}$ admits an exact solution using mean-field theory in the
limit $N_{s}\to\infty$, and features a paramagnet (PM) to FM phase transition
[116, 117, 118, 24, 119, 120] at $\Delta=J$. The effective interactions
between atoms within the same cluster in Eq. (5) are purely ferromagnetic and,
to leading order, we can identify $J=4g^{2}\eta/(1+\eta)\omega_{c}$ after
disorder averaging, with $\eta=M/N$. Each cluster is further coupled to other
clusters via Eq. (6), which is expected to generate frustration in the system.
For $N_{s}=1$, the ferromagnetic interaction vanishes since
$(S^{x})^{2}=(\sigma^{x})^{2}=1$ and this model becomes the quantum Hopfield
model (QHM) [121]. The QHM has a PM ground-state for sufficiently large
$\Delta$ while for small $\Delta$, the ground-state crucially depends on the
ratio $\eta$ [122]. For small values $\eta<\eta_{c}\sim O(10^{-1})$, the
system is in the memory retrieval phase [123, 124, 108, 43], which is a Dicke
model in disguise with multiple superradiant/FM ground-states. When the number
of photon modes ($M$) surpasses a critical limit $\eta>\eta_{c}$ [122],
frustrations dominate and turn the system into a quantum glass [96, 97, 125,
109], arguably in the same universality class of the quantum Sherrington-
Kirkpatrick [103, 126] (SK) model. In this paper, we are interested in SG
dynamics and will only consider the limit $\eta=1$.

The many-body nature of the model (1) when photons participate in the dynamics
prevents us to use exact diagonalization. Moreover, because of the frustrated
couplings, mean-field (MF) methods or dynamics of cumulants expansions (CE)
[23] are inapplicable. Instead, we attack this problem using methods of non-
equilibrium quantum field theory (NEQFT). In the next section, we will
introduce the method using a simple example first, and then will proceed to
apply it to the model in Eq. (1). A treatment of the Dicke model without
disordered couplings is also provided in Appendix E, for comparison with the
random Dicke model and pedagogical purposes.

##  III Non-Equilibrium Field Theory

###  III.1 The Action Principle for Quantum Mechanics

There are various formulations of classical mechanics [127, 128]. Albeit
equivalent in terms of their results, one of them is often more convenient to
use, depending on the specific problem that is to be solved. This is
especially true when an exact solution for the problem is not available and
approximations have to be made. One example is the problem of planetary
motion, where Hamiltonian mechanics is the formalism of choice as it admits
the powerful method of canonical perturbation theory [128]. On the other hand,
Lagrangian formalism, based on an action principle, is preferred when dealing
with classical field theories and systems with locally conserved charges due
to continuous symmetries.

Quantum mechanics has the counterparts to Hamiltonian and Lagrangian
formulations in classical mechanics. The Hamiltonian picture of quantum
mechanics is given by the Heisenberg equation of motion

$-i\hbar\partial_{t}\hat{\mathcal{O}}=\commutator{\hat{H}}{\hat{\mathcal{O}}}.$
(8)

However, the action principle approach to quantum mechanics is lesser known.
It was introduced originally by Heisenberg and Euler [129] and Schwinger [130]
in the context of quantum electrodynamics (QED) to study the modifications to
the Maxwell theory by quantum fluctuations. For a single-particle system
following unitary evolution, the quantum effective action $\Gamma$ is given by
[131]

$\Gamma\equiv\int_{t_{i}}^{t_{f}}\bra{\psi_{-},t}i\hbar\partial_{t}-\hat{H}\ket{\psi_{+},t}\,dt,$
(9)

where $\ket{\psi_{\pm},t}$ are arbitrary vectors satisfying the constraints

$\displaystyle\ket{\psi_{+},t_{i}}=\ket{\psi_{-},t_{i}}=\ket{\psi_{i}},\quad\innerproduct{\psi_{-},t}{\psi_{+},t}=1,$
(10)
$\displaystyle\bra{\psi_{-},t}\hat{x}\ket{\psi_{+},t}=X(t),\quad\bra{\psi_{-},t}\hat{p}\ket{\psi_{+},t}=P(t),$
(11)

where $\hat{x}$ and $\hat{p}$ are the position and momentum operators. It can
be shown that $\Gamma$ is only a functional of $X$ and $P$ using Lagrange
multipliers [131]. Then, the equation of motion for $X$ is found from the
stationary solution of $\Gamma$:

$\frac{\delta\Gamma}{\delta X}=\frac{\delta\Gamma}{\delta P}=0.$ (12)

Eq. (9) can be extended to non-unitary dynamics given by the Lindblad master
equation [132, 133]

$\partial_{t}\hat{\rho}=\hat{\mathcal{L}}\quantity(\hat{\rho}),$ (13)

with the effective action defined by

$\Gamma\equiv
i\hbar\int_{t_{i}}^{t_{f}}\Tr(\partial_{t}\hat{\rho}-\hat{\mathcal{L}}\hat{\rho})\,dt,$
(14)

where $\rho$ is an arbitrary density matrix satisfying

$\displaystyle\hat{\rho}(t_{i})=\hat{\rho}_{0},$ (15)
$\displaystyle\Tr(\hat{\rho}(t)\hat{x})=X(t),\quad\Tr(\hat{\rho}(t)\hat{p})=P(t),$
(16)

where $\hat{\rho}_{0}$ is the initial density matrix. We have to make two
remarks here. First, the quantum action principle was defined _independently_
from the path integral language as everything above was expressed in terms of
state vectors. However, $\Gamma$ can be also defined using path integrals as
will be shown later. Second, in contrast to classical actions, $\Gamma$ can
have an imaginary part [134, 88, 135, 136]. This happens when dissipation is
present in the system, either because of the explicit dissipative dynamics due
to coupling to a Markovian bath in Eq. (13), or when unitary dynamics involve
processes that result in dissipation, such as the real process of electron-
positron pair production in QED [129], or collisions between (quasi-)
particles in a many-body system.

Needless to say that $\Gamma$ contains the full quantum dynamics of the system
and for a many-body system, an exact calculation of $\Gamma$ and solving the
resulting equations of motion is still a cumbersome task. However, similarly
to the classical case, the effective action admits new methods of
approximation. These include, for instance, the time-dependent variational
principle (TDVP) [131, 137] and field theoretical methods based on quantum
effective actions [63, 61, 87] which are the subject of the next section.

###  III.2 The Keldysh Formalism

We shortly discuss the principal elements of Keldysh formalism of quantum
mechanics, which is essential to the study of quantum systems out of
equilibrium. Another advantage of this formalism for our problem is its
applicability to systems with static disorder, without resorting to
traditional approaches such as the replica trick. Interested readers can use
Refs. [87, 101] for a detailed exposition of the subject. For a system with
unitary time evolution generated by a stationary Hamiltonian $H$, the density
matrix at time $t$ is related to the initial density matrix $\rho_{0}$ via

$\rho(t)=e^{-iHt}\rho_{0}e^{+iHt}.$ (17)

For simplicity, we assume that $H$ describes a single particle system with
canonical operators $(x,p)$ satisfying $\commutator{x}{p}=i$. The eigenstates
of $x$ form a basis for the Hilbert space and $\rho(t)$ can be expanded in
terms of them:

$\bra{x}\rho(t)\ket{x^{\prime}}=\bra{x}e^{-iHt}\rho_{0}e^{+iHt}\ket{x^{\prime}}.$
(18)

We then divide $t$ into infinitesimal small pieces

$e^{\pm iHt}=e^{\pm iH\delta t}e^{\pm iH\delta t}\dots e^{\pm iH\delta
t}e^{\pm iH\delta t}.$ (19)

Substituting it in (18) and inserting resolutions of identity between
operators give

$\bra{x}\rho(t)\ket{x^{\prime}}=\int\prod_{i}^{N}dx^{+}_{i}\,dx^{-}_{i}\bra{x}e^{-iH\delta
t}\ket{x^{+}_{N}}\bra{x^{+}_{N}}e^{-iH\delta
t}\ket{x^{+}_{N-1}}\times\dots{\\\ }\times\dots\bra{x^{+}_{1}}e^{-iH\delta
t}\ket{x^{+}_{0}}\bra{x^{+}_{0}}\rho_{0}\ket{x^{-}_{0}}\bra{x^{-}_{0}}e^{+iH\delta
t}\ket{x^{-}_{1}}\times\dots\times\bra{x^{-}_{N-1}}e^{+iH\delta
t}\ket{x^{-}_{N}}\bra{x^{-}_{N}}e^{+iH\delta t}\ket{x^{\prime}}.$ (20)

After expanding the exponentials and following the same procedure as the
conventional Feynman path integral [138], except for the two chains of
operators on RHS and LHS of $\rho_{0}$, we get

$\bra{x}\rho(t)\ket{x^{\prime}}=\int\mathcal{D}\quantity[x_{+},x_{-}]e^{+iS_{+}}\bra{x_{+}(0)}\rho_{0}\ket{x_{-}(0)}e^{-iS_{-}},$
(21)

where the forward ($S^{+}$) and backward ($S^{-}$) branches of Keldysh action
are defined as

$S_{\pm}=\int_{0}^{t}dt^{\prime}\,\quantity(\frac{1}{2}m\dot{x}_{\pm}^{2}-V(x_{\pm})),$
(22)

together with constraints

$x_{+}(t)=x\,\quad x_{-}(t)=x^{\prime}.$ (23)

The main difference between Keldysh path integral and the ground-state path
integral is that in the former, fields carry an extra $s=\pm$ index which
indicates whether a field lies on the forward $+$ or backward $-$ branches of
the time contour. Therefore, the number of fields is doubled. For unitary time
evolution, Keldysh action is diagonal in the basis of time direction, and
fields on the opposite branches are coupled to each other only by the initial
state through the factor $\bra{x_{+}(0)}\rho_{0}\ket{x_{-}(0)}$ in Eq. (21).
For open systems, such as those described by the Lindblad master equation,
forward and backward fields are directly coupled through extra terms in the
Keldysh action [88].

In accordance with the fields, correlation functions carry extra indices, too.
For example

$G_{ss^{\prime}}(t,t^{\prime})=-i\expectationvalue{x_{s}(t)x_{s^{\prime}}(t^{\prime})}_{c},\quad
s,s^{\prime}=\pm.$ (24)

It can be shown that [101, 87] the correlation functions of the fields are
related to the correlation functions of the original operators in the
following way

$\displaystyle G_{+-}(t,t^{\prime})$
$\displaystyle=-i\expectationvalue{x(t^{\prime})x(t)},$ (25) $\displaystyle
G_{-+}(t,t^{\prime})$ $\displaystyle=-i\expectationvalue{x(t)x(t^{\prime})},$
(26) $\displaystyle G_{++}(t,t^{\prime})$
$\displaystyle=\Theta(t,t^{\prime})G_{-+}(t,t^{\prime})+\Theta(t^{\prime},t)G_{+-}(t,t^{\prime}),$
(27) $\displaystyle G_{--}(t,t^{\prime})$
$\displaystyle=\Theta(t,t^{\prime})G_{+-}(t,t^{\prime})+\Theta(t^{\prime},t)G_{-+}(t,t^{\prime}).$
(28)

where $\Theta(t,t^{\prime})\equiv\Theta(t-t^{\prime})$ stands for the
Heaviside theta. A compact notation is achieved if we assign $\pm$ labels to
the two different branches of a closed time contour in the complex plane, such
that

$x_{s}(t)\to x(t_{c}),$ (29)

where $t_{c}$ specifies both $s$ and $t$. The Keldysh action becomes

$S=S_{+}-S_{-}=\oint dt_{c}\quantity(\frac{1}{2}m\dot{x}^{2}-V(x)),$ (30)

and the same notation is used for correlation functions.

Except for simple cases, for example when $V(x)$ is quadratic in $x$, it is
not possible to work with (30) without any approximations. In the next
section, we will show how one can build a systematic approximation for (30)
and even more complicated problems, to learn about their dynamics.

Figure 3: The leading order diagrams for $\Gamma_{2}$ in $\hbar$ expansion.

###  III.3 –

Basics of 2PI In Sec. III.1, we showed how to find the effective action for a
quantum system in terms of the expectation values of its operators. The
effective action can be generalized to also explicitly depend on the
correlation functions of the theory. In that case, the stationary solution of
$\Gamma$ yields also the equations of motion for the correlation functions
[63]. The simplest case of these effective actions, which is used in this
work, depends only on the 1-point and 2-point functions. For technical reasons
that will be explained below, $\Gamma$ is called the 2-particle irreducible
(2PI) effective action (EA) in this case. Keldysh formalism of the previous
section allows us to efficiently calculate $\Gamma$. A brief description of
the approach is given below, while a comprehensive introduction can be found
in Ref. [62].

####  III.3.1 Basic Definitions

For simplicity, we use the toy model of a single particle in Eq. (30) in an
anharmonic potential to pedagogically introduce the derivation of Dyson
equations (DE) from the 2PI formalism. The potential is given by

$V(x)=\frac{1}{2}\omega^{2}x^{2}+\frac{u}{4!}x^{4}.$ (31)

We also re-scale $x$ such that $m\equiv 1$. In the next step, we introduce
source terms by [62]

$S\to S+\oint dt\,J_{s}(t)x_{s}(t)+\frac{1}{2}\oint\oint
dt\,dt^{\prime}\,R(t,t^{\prime})x(t)x(t^{\prime}).$ (32)

We further define the generating functional $\mathcal{W}[J,R]$ in terms of the
trace of $\rho(t)$ in the presence of sources

$e^{i\mathcal{W}}\equiv\int\mathcal{D}[x]e^{iS+i\int J\cdot
xdt+\frac{i}{2}\int\int dtdt^{\prime}x\cdot R\cdot x}.$ (33)

It is easy to show that

$\frac{\delta\mathcal{W}}{\delta J(t)}=\expectationvalue{x(t)}\equiv X(t),$
(34)

and that the connected 2-point correlation function, Eq. (24, can be obtained
from

$\frac{\delta\mathcal{W}}{\delta
R(t,t^{\prime})}=\frac{1}{2}\Big{(}X(t)X(t^{\prime})+iG(t,t^{\prime})\Big{)}.$
(35)

The effective action in terms of $X$ and $G$ is given by the Legendre
transform of $\mathcal{W}$ [62]

$\Gamma[X,G]=\mathcal{W}[J,R]-\oint dtJ(t)X(t)\\\ -\frac{1}{2}\oint\oint
dt\,dt^{\prime}\,R(t,t^{\prime})\Big{(}X(t)X(t^{\prime})+G(t,t^{\prime})\big{)}.$
(36)

In principle, $J$ and $R$ have to be found in terms of $X$ and $G$ using Eqs.
(34) and (35) and then be put into Eq. (36). However, a direct evaluation of
the functional derivatives shows that

$\displaystyle\frac{\delta\Gamma}{\delta X(t)}$ $\displaystyle=-J(t)-\oint
dt\,R(t,t^{\prime})X(t^{\prime}),$ (37)
$\displaystyle\frac{\delta\Gamma}{\delta G(t,t^{\prime})}$
$\displaystyle=-\frac{i}{2}R(t^{\prime},t),$ (38)

and therefore, in the absence of sources, the equations of motion are given by
$\delta\Gamma/\delta X=\delta\Gamma/\delta G=0$ as expected.

Figure 4: Schematic representation of different processes captured by the
first few terms in the expansion of the 2PI effective action for a single
quantum mechanical particle in an anharmonic potential. Green circles
represent the excitations (’phonons’) generated on top of $X(t)$. Left panel:
At the leading order (LO), the coupling of the classical mean field solution
$X$ to fluctuations leads to time-dependent effective mass for $X$ and a time-
dependent gap for fluctuations. The next-to-leading order (NLO) corrections
take into account the generation and annihilation of fluctuations by the
classical mean field solution $X$ (center panel) and their interactions (right
panel), similar to the scattering of phonons from each other. These effects
lead to dissipation and equilibration of the system in both mean field and
fluctuation sectors.

The main complication of finding $\Gamma$ is the lack of explicit relations
for $(J,R)$ in terms of $(X,G)$. Fortunately, a diagrammatic representation
for $\Gamma$ exists whose derivation can be found in the literature [63, 62,
87]. The final result is given by

$\Gamma[X,G]=S[X]-\frac{i}{2}\Tr\ln
G^{-1}+\frac{i}{2}\Tr(G_{0}^{-1}G)+\Gamma_{2}.$ (39)

For the anharmonic oscillator, $S$ is the action in (30) and $G_{0}$ is the
connected Green’s function of the system given by
$G_{0}(t_{c},t^{\prime}_{c})=-i\expectationvalue{\delta x(t_{c})\delta
x(t^{\prime}_{c})}$ after we expand the action around the expectation value,
$X$, of $x=X+\delta x$, and keep the result to quadratic order in $\delta x$.
$G_{0}^{-1}$ is given by the kernel of the action for a particle in the
potential of Eq. (31) after expanding the action around a classical solution
$X$ and keeping terms up to quadratic order in $\delta x$

$G_{0}^{-1}(t,t^{\prime})=\delta(t,t^{\prime})\quantity(-\partial_{t}^{2}-\omega^{2}-\frac{u}{4}X^{2}).$
(40)

The trace terms in Eq. (39) describe the contribution of Gaussian fluctuations
on top the classical solution given by $X$. Finally, the last term in (39)
contains non-Gaussian contributions. It can be shown that [63, 62, 87]
$\Gamma_{2}$ is given by the sum of all vacuum bubbles, e.g. Feynman diagrams
without external legs, which cannot be split into disconnected parts by
cutting at most two of their lines, hence the name ’two-particle irreducible’
(2PI) for the method. Some of the first few diagrams of $\Gamma_{2}$ are shown
in Fig. 3. The diagrams in Fig. 3 are sorted in terms of their scaling with
the Planck constant. The number of loops in each diagram determines the power
$\hbar$ that diagram scales with [139, 140, 62, 141]. For instance, the first
two diagrams in Fig. 3 are $\mathcal{O}(\hbar^{2})$ and the third one is
$\mathcal{O}(\hbar^{3})$. This suggests an expansion of $\Gamma_{2}$ in powers
of $\hbar$. The limit $\hbar\to 0$ corresponds to the classical solution, and
by including higher order terms in $\hbar$ we take quantum effects into
account. We again remark that this procedure is not a perturbative expansion.
It is equivalent to a resummation of an infinite number of terms, as each
Green’s function in Fig. 3 is renormalized by interactions and has to
evaluated self-consistently. $\hbar$ is not always a good expansion parameter.
For example, in many-body systems quantum effects can be amplified due to the
many degrees of freedom strongly interacting with each other. In these
situations, another control parameter for the expansion is required, such as
the inverse of the number of particle species [61, 64, 62, 79] or the inverse
of spin size, which we will use later in Section IV. For later reference, we
give the expressions for diagrams in Fig. 3

$\displaystyle\Gamma_{A}$
$\displaystyle=+\frac{u}{8}\oint\quantity[G(t,t)]^{2}\,dt,$ (41)
$\displaystyle\Gamma_{B}$ $\displaystyle=-\frac{u^{2}}{12}\oint\oint
X(t)\quantity[G(t,t^{\prime})]^{3}X(t^{\prime})\,dt\,dt^{\prime},$ (42)
$\displaystyle\Gamma_{C}$
$\displaystyle=+\frac{iu^{2}}{48}\oint\oint\quantity[G(t,t^{\prime})]^{4}\,dt\,dt^{\prime}.$
(43)

####  III.3.2 Equations of Motion

The equations of motion are given by Eqs. (37) and (38). For an expansion of
$\Gamma_{2}$ to third order in $\hbar$, we have for the dynamics of $X$

$\big{(}\overbrace{\partial_{t}^{2}+\omega^{2}+\frac{u}{6}X^{2}}^{\sim\hbar^{0}}+\overbrace{\frac{iu}{2}G(t,t)}^{\sim\hbar}\big{)}X(t)=\\\
\overbrace{\frac{u^{2}}{6}\oint\quantity[G(t,t^{\prime})]^{3}X(t^{\prime})\,dt^{\prime}}^{\sim\hbar^{2}}.$
(44)

This equation is similar to the classical equation of motion for $X$ with
extra terms that capture quantum corrections to the dynamics. The dynamics of
$X$ can be imagined as the motion of a classical particle interacting with the
fluctuations (also referred as ’phonons’ for pictorial purposes) around the
mean field solution. At first order (left panel in Fig. 4), the effective mass
of $X$ becomes time-dependent due to a shift given by
$\delta\omega^{2}(t)=iuG(t,t)/2$. Furthermore, the fluctuation gap becomes
dependent on $X$. Since $X$ changes with time, the fluctuation gap also
becomes time-dependent, and $X$ acts as a coherent drive for fluctuations.
Equivalently, the dynamics is approximated by a time-dependent Gaussian ansatz
for the variational state in Eq. (9) [131], which is equivalent to the
Bogoliubov approximation [101] for the dynamics. The effect of higher order
corrections can be imagined as the incoherent creation and absorption of
fluctuations by $X$ (middle panel of Fig. 4), given by the RHS of Eq. (44) as
a memory term. The advantage of 2PI is that it takes into account these
processes non-perturbatively by resumming their contributions to infinite
order. The equation of motion for $G$ always has the following general form

$G_{0}^{-1}G(t,t^{\prime})=\delta_{c}(t,t^{\prime})+\oint\Sigma(t,t^{\prime\prime})G(t^{\prime\prime},t^{\prime})\,dt^{\prime\prime},$
(45)

known as the Dyson equation [101]. $\delta_{c}(t,t^{\prime})$ is the delta
function for contour variables. $\Sigma$ is the self-energy functional and is
defined as

$\Sigma(t,t^{\prime})\equiv 2i\frac{\delta\Gamma_{2}}{\delta
G(t^{\prime},t)}.$ (46)

Pictorially, $\Sigma$ is obtained by removing a single line from diagrams in
$\Gamma_{2}$. Evaluating $\Sigma$ and using Eq. (40) gives

$-\Big{(}\partial_{t}^{2}+\omega^{2}+\frac{u}{2}X^{2}+\frac{iu}{2}G(t,t)\Big{)}G(t,t^{\prime})=\delta_{c}(t,t^{\prime})\\\
+\oint\tilde{\Sigma}(t,t^{\prime\prime})G(t^{\prime\prime},t^{\prime})\,dt^{\prime\prime},$
(47)

where we have separated the non-local part of $\Sigma$ as $\tilde{\Sigma}$.
$\tilde{\Sigma}$ is the sum of $\Gamma_{B}$ and $\Gamma_{C}$ which are
respectively obtained using the functional derivatives of $\Gamma_{B}$ and
$\Gamma_{C}$ in Eqs. (42) and (43)

$\tilde{\Sigma}(t,t^{\prime})=-\frac{iu^{2}}{2}X(t)\quantity[G(t,t^{\prime})]^{2}X(t^{\prime})-\frac{u^{2}}{6}\quantity[G(t,t^{\prime})]^{3}.$
(48)

The last term in Eq. (48) describes the ’collisions’ between excitations
(right panel of Fig. 4) and usually leads to redistribution of energy and
thermalization of system at long times. Eqs. (44) and (47) are called Dyson
equations [87, 101]. Usually, DE can be solved only approximately using
analytical approaches. Approximate methods can usually access dynamics at
early times, where a perturbative expansion around the initial state is valid,
and near equilibrium, where time translation symmetry is approximately true.
To capture far from equilibrium dynamics at all stages of the evolution we
usually need to solve DE numerically without further approximations [62]. The
latter approach has shown promising results for dynamics in various physical
systems out of equilibrium [142, 65, 66, 67, 68, 69, 70, 73, 143, 75, 76, 77,
78, 79, 80, 81, 82, 144, 83, 145]).

####  III.3.3 Comparison with other approaches

There are various methods to explore many-body quantum dynamics theoretically
and each one has its own advantages and shortcomings. Exact diagonalization
(ED) gives accurate results but is often limited to very small system sizes,
especially for open quantum systems where the size of the vector space grows
even faster due to the necessity of working with mixed states. Methods based
on matrix product states (MPS) [146] are mostly limited to one spatial
dimension and systems with local interactions and weak entanglement. Among the
most frequently used methods in the AMO community are cumulants expansion (CE)
and truncated Wigner approximation (TWA) [147] together with its extension,
discrete truncated Wigner approximation (DTWA) [148]. The main advantages of
CE are simplicity and cheap computational cost. On the other hand, it is an
uncontrolled approximation [149, 150]. There is no a priori knowledge of its
domain of applicability before solving the equations and check the physical
consistency of the results. Furthermore, to calculate correlation functions at
different times in CE, one needs to resort to the quantum regression theorem
[114, 151], which complicates the calculations. TWA approximates quantum
dynamics with classical statistical mechanics by sampling the initial
probability distribution function from the system’s initial wave function and
subsequently, evolving the system according to the classical equations of
motion. The advantage of TWA is that it is a controlled approximation, since
it can be expressed as the leading order contribution in the expansion of
dynamics in powers of $\hbar$, which turn out to be a weighted average of
classical dynamics [147]. The drawback of TWA is its limitation to systems
with weak quantum fluctuations and when quantum effects are not built up over
time. For instance, TWA is unable to capture tunneling phenomena [147] even at
the level of qualitative accuracy. While the calculation of non-local
symmetric correlation functions is straightforward in TWA, to evaluate
quantities such as response functions one needs to go to higher order terms in
$\hbar$ [147], considerably increasing the required effort to use the method.
In comparison to CE and TWA, 2PI can be used to perform controlled
approximations, provided that a control parameter exists in the system. In
this case, the 2PI action admits an expansion in powers of the control
parameter [62]. For instance, this parameter can be $\hbar$, similar to the
example given before, or inverse of the components of a vector field or in our
case, the inverse of the spin amplitude per cluster and the number of photon
modes to the number of clusters. Moreover, symmetric and anti-symmetric
correlation functions are the quantities in terms of which the formalism is
built and are the direct outcomes of calculations. 2PI also excels in
capturing quantum effects, mainly because it involves a resummation of the
perturbative expansion to infinite order. For instance, it has been used to
study dynamics in strongly correlated systems of electrons [143, 152, 153] and
phonons [145, 154] with non-Fermi liquid (NFL) behavior and critical
fluctuations. 2PI works well also when small quantum effects are accumulated
over time, leading to drastic changes in the system at long times such as in
tunneling phenomena [65, 155]. Despite numerous advantages, 2PI has some
limitations. First, the approximations that are usually made to retain only a
subset of the diagrams in the effective action, such as $1/N$ expansions [64,
61, 62, 79], give qualitatively valuable results about the universal trend of
the dynamics, but are not tailored to have quantitative accuracy, i.e., they
are not suitable for a point-wise comparison with experimental data.
Sometimes, one needs to make educated guesses about which diagrams have to be
kept to capture a certain aspect of the physics which is of interest. The
exception to these is working in the weak coupling or the dilute limit, where
collisions can be incorporated perturbatively [101, 87]. Second, working with
non-Gaussian initial states is difficult in 2PI as these require the inclusion
of extra interaction vertices [62] that complicate the approximation. We
emphasize that this restriction only holds for initial states. 2PI is not
limited to Gaussian dynamics (such as second order cumulants) and in fact,
captures non-Gaussian correlations generated over time after initializing the
systems in a Gaussian state.

##  IV 2PI for Disordered Cavity-QED

Figure 5: The spin-photon interaction vertex after averaging over the
disorder. The dashed line indicates the average
$\overline{g_{\alpha,i}g_{\beta,j}}$.

$\Gamma_{2}\hskip 5.0pt\approx\hskip
5.0pt\underset{\mbox{\scalebox{1.0}{$\Gamma_{2}^{\mathrm{LO}}\sim
N_{s}$}}}{\vbox{\hbox{\includegraphics[height=21.68231pt]{MF.png}}}}\hskip
5.0pt+$
$\underbrace{\Bigg{(}\vbox{\hbox{\includegraphics[height=34.69038pt]{melon.png}}}\hskip
2.0pt+\hskip
2.0pt\vbox{\hbox{\includegraphics[height=34.69038pt]{2.png}}}\hskip
2.0pt+\hskip
2.0pt\vbox{\hbox{\includegraphics[height=34.69038pt]{3.png}}}\hskip
2.0pt+\hskip
2.0pt\vbox{\hbox{\includegraphics[height=34.69038pt]{4.png}}}\hskip
2.0pt+\,\,\dots\Bigg{)}}_{\mbox{\scalebox{1.0}{$\Gamma_{2}^{\mathrm{NLO}}\sim
N_{s}^{0}$}}}$

Figure 6: The leading-order and next-to-leading-order diagrams in the
$1/N_{s}$ expansion of 2PI action.

###  IV.1 Keldysh Action for Spins

In order to treat Eq. (1) using field theory, we need a path integral
representation for spin operators. Different spin representations include spin
coherent state path integral and representations in terms of bosons or
fermions [156]. In this work, we represent each spin-half operator in terms of
three Majorana fermions $(\psi^{x},\psi^{y},\psi^{z})$ as [157, 158, 159, 160,
161, 162, 163, 164, 165, 166, 167, 168, 169] given by

$\sigma^{\alpha}=-i\epsilon_{\alpha\beta\gamma}\psi^{\beta}\psi^{\gamma},\qquad\anticommutator{\psi^{\alpha}}{\psi^{\beta}}=\delta_{\alpha\beta},$
(49)

where we have assumed summation over repeated indices. It is easy to check
that (49) satisfies spin commutation relations
$\commutator{\sigma^{\alpha}}{\sigma^{\beta}}=2i\epsilon_{\alpha\beta\gamma}\sigma^{\gamma}$.
This representation was used by Refs. [168, 169] to study the onset of
superradiance in the steady state of the Dicke model with different types of
external baths. The Keldysh action for “free” spins, corresponding to the
first term in Eq. (1) has two parts

$\displaystyle S_{\sigma}$ $\displaystyle=S_{\mathrm{B}}+S_{\Delta},$ (50)
$\displaystyle S_{\mathrm{B}}$
$\displaystyle=\sum_{i}^{N}\sum_{\lambda}^{N_{s}}\sum_{\alpha}^{x,y,z}\oint\frac{i}{2}\psi^{\alpha}_{i\lambda}\partial_{t_{c}}\psi^{\alpha}_{i\lambda}\,dt_{c},$
(51) $\displaystyle S_{\Delta}$
$\displaystyle=\sum_{i}^{N}\sum_{\lambda}^{N_{s}}\oint
i\Delta\psi^{x}_{i\lambda}\psi^{y}_{i\lambda}\,dt_{c},$ (52)

where $S_{B}$ is the contribution of the Berry phase of spins to the action
[170]. The action in Eq. (50) can be compactly written as

$S_{\sigma}=\frac{1}{2}\sum_{i}^{N}\sum_{\lambda}^{N_{s}}\oint\Psi^{T}_{i\lambda}\hat{G}_{0}^{-1}\Psi_{i\lambda}\,dt_{c},$
(53)

with
$\Psi^{T}_{i\lambda}\equiv(\psi^{x}_{i\lambda},\psi^{y}_{i\lambda},\psi^{z}_{i\lambda})$.
The inverse bare Green’s function for fermions $\hat{G}_{0}^{-1}$ is defined
as

$\hat{G}_{0}^{-1}\equiv\begin{bmatrix}i\partial_{t_{c}}&i\Delta&0\\\
-i\Delta&i\partial_{t_{c}}&0\\\ 0&0&i\partial_{t_{c}}\end{bmatrix}.$ (54)

Finally, we define the fermion Green’s function and its diagrammatic
representation as

$iG_{i\lambda,j\lambda^{\prime}}^{\alpha,\beta}(t,t^{\prime})\equiv\expectationvalue{\psi^{\alpha}_{i\lambda}(t)\psi^{\beta}_{j\lambda^{\prime}}(t^{\prime})}:\includegraphics{psipsi.pdf}.$
(55)

Note that $G$ is the dressed fermion Green’s function.

###  IV.2 Keldysh Action for Photon Sector

The dissipative Keldysh action for photons can be obtained directly from the
Liouvillian, by following the prescription given in Ref. [88]

$S_{\mathrm{ph}}=\int\big{[}\bar{a}_{+}\quantity(i\partial_{t}-\omega_{c}+i\kappa)a_{+}-\bar{a}_{-}\quantity(i\partial_{t}-\omega_{c}-i\kappa)a_{-}\\\
-2i\kappa\bar{a}_{-}a_{+}\big{]}\,dt.$ (56)

Since the spin-photon coupling in Eq. (1) depends on the combination
$(a+a^{\dagger})$ of photon operators, dealing with interactions is simpler
when we make the following transformation to real-valued photon fields
$(\phi,\pi)$ given by:

$a_{\pm}=\sqrt{\frac{\omega_{c}}{2}}\quantity(\phi_{\pm}+i\frac{\pi_{\pm}}{\omega_{c}}).$
(57)

Substitution in Eq. (56) gives

$S_{\mathrm{ph}}=\frac{1}{2}\sum_{\alpha}^{M}\int\Phi_{\alpha}^{T}\hat{D}_{0}^{-1}\Phi_{\alpha}\,dt,$
(58)

where
$\Phi_{\alpha}^{T}\equiv\quantity(\phi_{\alpha+},\pi_{\alpha+},\phi_{\alpha-},\pi_{\alpha-})$
and $\hat{D}_{0}^{-1}$ is a $4\times 4$ matrix defined as

$\displaystyle\hat{D}_{0}^{-1}$
$\displaystyle=\begin{bmatrix}\quantity(\hat{D}_{0}^{-1})^{++}&\quantity(\hat{D}_{0}^{-1})^{+-}\\\
\quantity(\hat{D}_{0}^{-1})^{-+}&\quantity(\hat{D}_{0}^{-1})^{--}\end{bmatrix},$
(59) $\displaystyle\quantity(\hat{D}_{0}^{-1})^{++}$
$\displaystyle=\begin{bmatrix}-\omega_{c}^{2}+i\kappa\omega_{c}&-\partial_{t}\\\
\partial_{t}&-1+i\frac{\kappa}{\omega_{c}}\end{bmatrix},$ (60)
$\displaystyle\quantity(\hat{D}_{0}^{-1})^{+-}$
$\displaystyle=\begin{bmatrix}-i\kappa\omega_{c}&-\kappa\\\
\kappa&-i\frac{\kappa}{\omega_{c}}\end{bmatrix},$ (61)
$\displaystyle\quantity(\hat{D}_{0}^{-1})^{-+}$
$\displaystyle=\begin{bmatrix}-i\kappa\omega_{c}&\kappa\\\
-\kappa&-i\frac{\kappa}{\omega_{c}}\end{bmatrix},$ (62)
$\displaystyle\quantity(\hat{D}_{0}^{-1})^{--}$
$\displaystyle=\begin{bmatrix}\omega_{c}^{2}+i\kappa\omega_{c}&\partial_{t}\\\
-\partial_{t}&1+i\frac{\kappa}{\omega_{c}}\end{bmatrix}.$ (63)

Photon Green’s functions are defined according to

$iD^{\rho,\rho^{\prime}}_{\alpha,\beta}(t,t^{\prime})\equiv\expectationvalue{\Phi^{\rho}_{\alpha}(t)\Phi^{\rho^{\prime}}_{\beta}(t^{\prime})}_{c},$
(64)

where $\rho,\rho^{\prime}=(\phi,\pi)$. We will see that only the
$D^{\phi\phi}$ component appears explicitly in the diagrams for the effective
action and self-energies. Hence, only $D^{\phi\phi}$ requires a diagrammatic
representation which is given by

$iD^{\phi\phi}_{\alpha,\beta}(t,t^{\prime})=\expectationvalue{\phi_{\alpha}(t)\phi_{\beta}(t^{\prime})}_{c}:\includegraphics{phiphi.pdf}.$
(65)

###  IV.3 Spin-Photon Interaction

Using the conventions introduced above, the Keldysh action for the interaction
term reads

$S_{\mathrm{int}}=2i\sqrt{\frac{2\omega_{c}}{(N+M)N_{s}}}\sum_{\alpha,i,\lambda}\oint
dt_{c}\,g_{\alpha i}\phi_{\alpha}\psi^{y}_{i\lambda}\psi^{z}_{i\lambda}.$ (66)

In principle, we can proceed by taking the average of the Keldysh action over
the random couplings $g_{\alpha i}$. This yields an effective interaction
defined by $e^{iS_{\mathrm{eff}}}=\overline{e^{iS_{\mathrm{int}}}}$ where

$S_{\mathrm{eff}}\equiv-\frac{4ig^{2}\omega_{c}}{(N+M)N_{s}}\sum_{\alpha,i}\sum_{\lambda,\lambda^{\prime}}\oint\oint
dt_{c}\,dt^{\prime}_{c}\,\phi_{\alpha}(t_{c})\phi_{\alpha}(t^{\prime}_{c})\\\
\times\psi^{y}_{i\lambda}(t_{c})\psi^{z}_{i\lambda}(t_{c})\psi^{y}_{i\lambda^{\prime}}(t^{\prime}_{c})\psi^{z}_{i\lambda^{\prime}}(t^{\prime}_{c}).$
(67)

We have assumed that the initial state is not correlated with disorder profile
(see comments in Sections IV.4 and IV.5 for more details). The diagrammatic
form of $S_{\mathrm{eff}}$ is given in Fig. 5. To have a systematic and
controlled approximation in $1/N_{s}$ that captures the frustrated nature of
the problem, we have to keep an infinite subset of 2PI diagrams shown in Fig.
6. A closed form for the corresponding summation can be found, as shown for
example for quantum $O(N)$ model in Refs. [64, 61, 62, 84]. An easier approach
is the auxiliary field method [64, 61, 62], based on the Hubbard-Stratonovich
(HS) transformation [101, 11, 79, 82]. HS transformation finds various
applications in the study of collective effects in many body systems such as
plasmons [11], superconductivity [11, 171], superfluidity [170] and quantum
spin liquids [172]. The basic idea is to introduce a new field which we label
as $\chi$, that mediates the original interaction in $S_{\mathrm{int}}$. In
our case, $\chi$ decouples the interaction between spins and cavity modes as
diagrammatically illustrated in Fig. 7a. The action of $\chi$ and its coupling
to other degrees of freedom are given by (see Appendix A for a mathematical
derivation)

$S_{\mathrm{int}}\to S_{\chi}+S_{\chi\psi}+S_{g\chi\phi},$ (68)

where $S_{\chi}$ is the action of the HS defined as

$S_{\chi}\equiv\frac{1}{2}\sum_{\alpha,\beta}^{M}\sum_{i,j}^{N}\sum_{\sigma,\sigma^{\prime}}^{1,2}\oint\oint
dt_{c}\,dt^{\prime}_{c}\,\\\ \times\chi^{\sigma}_{\alpha
i}(t_{c})\cdot\quantity(\hat{W}^{-1}_{0})^{\sigma,\sigma^{\prime}}_{\alpha
i,\beta j}(t_{c},t^{\prime}_{c})\cdot\chi^{\sigma^{\prime}}_{\beta
j}(t^{\prime}_{c}).$ (69)

We name $\chi$ as the Ising field, as it mediates the Ising-type interaction
amongst spins as we will see below. $\chi$ is a two-component real valued
scalar field defined as

$\vec{\chi}_{\alpha i}\equiv\quantity(\chi^{1}_{\alpha i},\chi^{2}_{\alpha
i}),$ (70)

together with its inverse bare Green’s function

$\quantity(\hat{W}^{-1}_{0})^{\sigma,\sigma^{\prime}}_{\alpha i,\beta
j}(t_{c},t^{\prime}_{c})\equiv\sqrt{N+M}\delta_{\alpha\beta}\delta_{ij}\quantity(\sigma^{x})_{\sigma\sigma^{\prime}}\delta(t_{c},t^{\prime}_{c}),$
(71)

and its full Green’s function

$iW_{\alpha i,\beta
j}^{\sigma,\sigma^{\prime}}(t,t^{\prime})\equiv\expectationvalue{\chi^{\sigma}_{\alpha
i}(t)\chi^{\sigma^{\prime}}_{\beta
j}(t^{\prime})}_{c}:\includegraphics{chichi.pdf}.$ (72)

Note that the “free” part of the action for the Ising field $S_{\chi}$ is
local in time and does not contain any time derivatives of $\chi$. This makes
the equations of motion for $\chi$ algebraic rather than differential. The
latter is the generic case where the equations of motion for correlation
functions form a system of coupled differential equations, and adiabatic
elimination is equivalent to approximately ignoring the time derivatives of
some of the dynamical variables, which are assumed to have a quick response
compared to other timescales in the system. For the HS field, adiabatically
eliminating $\chi$ is exact, which is equivalent to taking the Gaussian
integral over $\chi$ in Eq. (69), and the result is given by
$S_{\mathrm{int}}$ in Eq. (66). The next term in Eq. (68) is $S_{\chi\psi}$
which describes the coupling of fermions to the first component of $\chi$:

$S_{\chi\psi}\equiv-\frac{2i}{\sqrt{N_{s}}}\sum_{\alpha,i,\lambda}\oint
dt_{c}\,\chi^{1}_{\alpha i}\psi^{y}_{i\lambda}\psi^{z}_{i\lambda},$ (73)

and $S_{g\chi\phi}$ describes the disordered interaction of photons with the
second component of $\chi$

$S_{g\chi\phi}\equiv\sqrt{2\omega_{c}}\sum_{\alpha,i}\oint
dt_{c}\,\chi^{2}_{\alpha i}g_{\alpha i}\phi_{\alpha}.$ (74)

The diagrammatic representations of the original vertex in Eq. (66) and the
transformed ones (Eqs. (73) and (74) are given in Fig. 7a. We will show later
that the two components of the Ising field correspond to different physical
quantities. As will be shown in Section IV.8.1, $\chi^{1}$ is related to the
effective magnetic field each cluster experiences and $\chi^{2}$ is connected
to magnetization. Similarly, the Green’s functions of Ising fields are not
just mathematical objects and have physical meanings. $W$ can be expressed in
terms of the original Green’s functions as shown in Fig. 7c. $W^{22}$ is
related to the spin-spin correlation function or equivalently, the 4-point
function of Majorana fermions

$\expectationvalue{S^{x}_{i}(t)S^{x}_{i}(t^{\prime})}=-4\sum_{\lambda,\lambda^{\prime}}^{N_{s}}\expectationvalue{\psi^{y}_{i\lambda}(t)\psi^{z}_{i\lambda}(t)\psi^{y}_{i\lambda^{\prime}}(t^{\prime})\psi^{z}_{i\lambda^{\prime}}(t^{\prime})}$
(75)

whose leading order expansion is given by the same set of diagrams as $W^{22}$
in Fig. 7c, up to multiplication by an overall constant, as given by Eq.
(118). Therefore, spin-spin correlation functions are natural byproducts of
our formalism. Therefore, there is no need to solve the Bethe-Salpeter
equations to obtain 4-point functions of fermions, usually a cumbersome task
particularly for out of equilibrium systems [79, 173, 86].

(a)  (b)

(c)

Figure 7: (a) The original interaction vertex and the decoupled interactions
after HS transformation. (b) The effective interaction between photons and the
Ising field after disorder averaging. (c) Diagrammatic representation of the
Ising (HS) field propagators. The appearing fermion and photon lines are
assumed to be renormalized by interactions. Thick curls are the dressed
Green’s functions of $\chi$ and thin curls are the bare propagators of $\chi$
and only have constant off-diagonal elements. $W^{22}$ coincides with spin-
spin correlation function, up to an overall multiplicative factor.

###  IV.4 Disorder Averaging

We can now take the average over the disordered couplings. In the Keldysh
formalism, the average can be taken without resorting to the replica trick
[101]. The only term in the action depending on $g$ is $S_{g\chi\phi}$. The
effective interaction after disorder averaging is given by
$e^{iS_{\chi\phi}}=\overline{e^{iS_{g\chi\phi}}}$, where

$S_{\chi\phi}=ig^{2}\omega_{c}\sum_{\alpha i}\oint\oint
dt_{c}\,dt^{\prime}_{c}\,\chi^{2}_{\alpha
i}(t_{c})\phi_{\alpha}(t_{c})\chi^{2}_{\alpha
i}(t^{\prime}_{c})\phi_{\alpha}(t^{\prime}_{c}),$ (76)

is shown diagrammatically in Fig. 7b. We have to make an important remark
about the process of disorder averaging. In obtaining Eq. (76) we have assumed
that the initial state of the system is not correlated with the disorder. This
is valid for the initial states we consider in this paper. However, to study
phenomena such as associative memory in multi-mode cavity QED [121, 122, 124,
108, 174, 43], where a significant overlap of the initial spin configuration
is required for memory retrieval, one has to assume that the initial state
depends on $g_{\alpha i}$. In that case, disorder averaging will generate more
terms than Eq. (76), which couple the initial state to the interaction vertex
in Eq. (74).

###  IV.5 Symmetry Considerations

The symmetry structure of the model helps us to simplify the study of its
dynamical response. Originally, the Hamiltonian in Eq. (1) is invariant only
under a global $Z_{2}$ transformation that maps all spins and cavity modes
simultaneously according to

$\sigma^{x}_{i\lambda}\to-\sigma^{x}_{i\lambda},\qquad a_{\alpha}\to-
a_{\alpha}.$ (77)

According to the language of Ref. [175, 88, 176, 177], in the absence of
photon loss this is a quantum symmetry of the system with a conserved $Z_{2}$
charge. A quantum symmetry is a symmetry of the fields on each individual
Keldysh contour, while a classical symmetry is the invariance of the Keldysh
action under a simultaneous transformation of the fields on forward and
backward contours [88]. With photon loss, the quantum symmetry is demoted to a
classical symmetry without a conserved charge. However, starting from a
symmetric initial state, a classical symmetry still guarantees that the
symmetry will remain unbroken in the absence of symmetry breaking
perturbations. The symmetry structure of the model is enriched after disorder
averaging and using the fermion representation in Eq. (49). It can be easily
verified that the disorder-averaged Keldysh action has the following sets of
symmetries

  1. 1.

A local $Z_{2}$ gauge symmetry under the transformation

$\vec{\psi}_{i\lambda}\to-\vec{\psi}_{i\lambda},$ (78)

which holds for each spin separately. This symmetry is an artifact of
representing spins in terms of quadratic fermion operators, and is not
physical. The initial state or external forces cannot break this symmetry. The
important consequence of this symmetry is that

$G^{\alpha,\beta}_{i\lambda,j\lambda^{\prime}}\propto\delta_{ij}\delta_{\lambda\lambda^{\prime}}.$
(79)

  2. 2.

A $Z_{2}$ symmetry for each separate cluster $i$ and Ising fields coupled to
it:

$\displaystyle\sigma^{x}_{i\lambda}$
$\displaystyle\to-\sigma^{x}_{i\lambda},\quad(\lambda=1,...,N_{s}),$ (80)
$\displaystyle\chi^{\sigma}_{\alpha i}$
$\displaystyle\to-\chi^{\sigma}_{\alpha
i},\quad(\sigma=1,2),\quad(\alpha=1,...,M).$ (81)

(81) holds because the effective interaction in Eq. (76) is quadratic in
$\chi$. This symmetry is a classical symmetry with no conserved quantities.

  3. 3.

A $Z_{2}$ symmetry of each photon mode given by

$\phi_{\alpha}\to-\phi_{\alpha},\quad\pi_{\alpha}\to-\pi_{\alpha}.$ (82)

This symmetry is also a result of $S_{\chi\phi}$ being quadratic in photon
fields and is a weak symmetry. This symmetry implies that

$D_{\alpha,\beta}^{\rho,\rho^{\prime}}\propto\delta_{\alpha\beta}.$ (83)

We see that the $Z_{2}$ symmetries of spin and photon sectors are decoupled.
This means that, even if the initial state of spins breaks the symmetry, no
photon coherence will be generated
($\overline{\expectationvalue{a_{\alpha}(t)}}=0$). On the other hand, a finite
value for $\sigma^{x}_{i\lambda}$ results in a finite value for
$\chi^{1,2}_{\alpha i}$. We again remark that the above arguments hold true
only if the initial state of the system is not correlated with the disorder
pattern, such that Eq. (76) is the only outcome of disorder averaging.
Otherwise, a $Z_{2}$ broken initial state can in principle break the $Z_{2}$
symmetry of some of the photon modes. This happens for example, if the initial
spin configuration has a strong overlap with a single disorder pattern
corresponding to the photon mode $\alpha$, such that

$\lim_{N\to\infty}\frac{1}{NN_{s}}\absolutevalue{\sum_{i}^{N}\sum_{\lambda}^{N_{s}}\expectationvalue{\sigma^{x}_{i\lambda}}_{0}g_{\alpha
i}}>0,$ (84)

or if a symmetry breaking perturbation that favors a single pattern such as

$\delta H=\epsilon\sum_{i,\lambda}g_{\alpha,i}\sigma^{x}_{i\lambda},$ (85)

is applied to the system. In this case, one expects that for sufficiently
small $M/N$, the pattern $\alpha$ to be activated and retrieved [122, 108].
For the fully polarized initial states of spins considered in this problem and
in the thermodynamic limit, we can safely take
$\overline{\expectationvalue{\phi_{\alpha}(t)}}=0$ throughout the evolution.
Even starting from a state with $\expectationvalue{\phi}\neq 0$, its value
will decay to zero as it cannot align itself with any of the disorder
patterns.

###  IV.6 2PI Action

The 2PI action for the model given above is a functional of fermion, photon
and Ising field correlation functions together with the expectation values of
Ising fields and has the general form similar to Eq. (39) given by

$\Gamma\quantity[\tilde{\chi},G,D,W]=S_{\chi}\quantity[\chi]+\frac{i}{2}\Tr\ln
G-\frac{i}{2}\Tr\quantity(G_{0}^{-1}G)\\\ -\frac{i}{2}\Tr\ln
D+\frac{i}{2}\Tr\quantity(D^{-1}_{0}D)-\frac{i}{2}\Tr\ln W\\\
+\frac{i}{2}\Tr\quantity(W^{-1}_{0}W)+\Gamma_{2}\quantity[\tilde{\chi},G,D,W].$
(86)

The expressions for $S_{\chi}$, $G_{0}$, $D_{0}$ and $W_{0}$ were respectively
given in Eqs. (69), (54), (59) and (71). $\tilde{\chi}$ is the expectation
value of the Ising field

$\tilde{\chi}^{\sigma}_{\alpha
i}(t)\equiv\overline{\expectationvalue{\chi^{\sigma}_{\alpha i}(t)}},$ (87)

shown by a black circle connected to a spring in (Fig. 8a). The last term in
Eq. (86) captures interactions and as we mentioned in Sec. III.3, is given by
the sum of 2PI diagrams. In order to systematically expand $\Gamma_{2}$, we
need to specify how the expectation values and Green’s functions of the Ising
field scale with parameters of the system. According to Eq. (71) to the
leading order in $(N+M)^{-1}$ we have

$W^{12}_{i\alpha,j\beta}\sim\delta_{ij}\delta_{\alpha\beta}\,\mathcal{O}\quantity(\frac{1}{\sqrt{N+M}}),$
(88)

The diagonal elements $W^{11}$ and $W^{22}$ are zero at the bare level in Eq.
(69). However, they become non-zero when the couplings of $\chi$ to $\psi$ and
$\phi$ are taken into account (Fig. 7c). For $W^{22}$ we have

$W^{22}_{i\alpha,j\beta}\sim\delta_{ij}\delta_{\alpha\beta}\,\mathcal{O}\quantity(\frac{1}{N+M}).$
(89)

As will be shown later, $W^{11}$ has a sub-leading term due to interactions
which contributes at leading order when it is summed over photon modes

$W^{11}_{\alpha\beta}\sim\delta_{ij}\delta_{\alpha\beta}\,\mathcal{O}\quantity(\frac{1}{N+M})\\\
+\delta_{ij}(1-\delta_{\alpha\beta})\,\mathcal{O}\quantity(\frac{1}{(N+M)^{2}}).$
(90)

Furthermore, $\tilde{\chi}$ will have the following scalings (Eqs. (113) and
(112))

$\tilde{\chi}^{1}\sim\mathcal{O}\quantity(\frac{\sqrt{N_{s}}}{N+M}),\quad\tilde{\chi}^{2}\sim\mathcal{O}\quantity(\sqrt{\frac{N_{s}}{N+M}}).$
(91)

At last, the fermion-Ising vertex in Eq. (73) has

$S_{\chi\psi}\sim\frac{1}{\sqrt{N_{s}}}.$ (92)

We now have all of the necessary ingredients to perform a systematic expansion
of $\Gamma_{2}$.

###  IV.7 Diagrammatic Evaluation of 2PI Action

(a)

(b)

Figure 8: (a) The leading-order contribution of interactions to 2PI action.
The dashed line in the right diagram is not a Green’s function and cannot be
cut. (b) A diagram which linearly scales with $N_{s}$ but is not extensive and
neglected. Black circles connected to springs represent the expectation value
of Ising field $\chi$.

We are interested in the thermodynamic limit of the system in Eq. (1) where
$N,M\to\infty$ while the ratio $\eta=M/N$ is kept fixed. Moreover, the number
of spins per cluster $N_{s}$ is assumed to be larger than one and will play
the role of the control parameter for the expansion. Also, we assume that
$N_{s}\ll N$, which is a valid assumption in the thermodynamic limit of the
problem. Similar to the quantum anharmonic oscillator considered in Section
III, the interaction part of the 2PI action, given by $\Gamma_{2}$ in Eq.
(86), admits a diagrammatic expansion in terms of the connected vacuum bubbles
of the theory which cannot be split into half by cutting one or two of their
Green’s function lines, also known as two-particle irreducible (2PI) graphs.
Below, we will classify these diagrams for our system as leading-order (LO)
terms

$\Gamma_{2}^{\mathrm{LO}}\sim N_{s},$ (93)

and next-to-leading-order (NLO) terms

$\Gamma_{2}^{\mathrm{NLO}}\sim N_{s}^{0},$ (94)

and higher order terms which are ignored in this work. We will also ignore
terms which are sub-extensive. The following discussion will also elucidate
how the parameter controlling the strength of quantum fluctuations, $N_{s}$,
enters naturally in the field theory description and in the DE derived from
it. This is one of the key merits of the approach, at variance with more
numerical oriented methods which have to deal with growing computational
complexity as $N_{s}$ is decreased.

####  IV.7.1 Leading-order contributions

The LO terms have linear scaling with $N_{s}$ and at the same time, scale
extensively with system size. Two of these diagrams exist and both involve the
expectation values of Ising fields, as shown in Fig. 8a. Their mathematical
expressions are given by

$\Gamma_{2}^{\mathrm{LO}}=\frac{2}{\sqrt{N_{s}}}\sum_{\alpha}^{M}\sum_{i}^{N}\sum_{\lambda}^{N_{s}}\oint
dt_{c}\,\tilde{\chi}^{1}_{\alpha
i}(t_{c})G^{y,z}_{i\lambda,i\lambda}(t_{c},t_{c})\\\
-g^{2}\omega_{c}\sum_{\alpha}^{M}\sum_{i}^{N}\oint\oint
dt_{c}\,dt^{\prime}_{c}\,\tilde{\chi}^{2}_{\alpha
i}(t_{c})D_{\alpha,\alpha}^{\phi\phi}(t_{c},t^{\prime}_{c})\tilde{\chi}^{2}_{\alpha
i}(t^{\prime}_{c}).$ (95)

Note that the disorder (dashed) line in Fig. 8a cannot be cut, as it is a part
of the disorder averaged interaction vertex in Fig. 7b. Accordingly, the right
diagram in Fig. 8a is not 2-particle reducible. It can be shown, by solving
the resulting equations of motion derived from Eq. (95), that the LO terms
describe a dynamical mean field interaction of spin expectation values
$\overline{\expectationvalue{\sigma^{x}}}$ mediated by photons through their
response function (see Section V.1). Therefore, the LO contribution describes
the LMG coupling in Eq. (5) with the inclusion of retardation effects due to
photon dynamics. We note that there are other terms that scale linearly with
$N_{s}$, such as the one given in Fig. 8b, but all of them scale non-
extensively with system size and can be neglected in the thermodynamic limit.

####  IV.7.2 Next-to-leading-order contributions

Figure 9: Next-to-leading-order contributions to 2PI action.

The NLO diagrams do not scale with $N_{s}$ but still scale linearly with the
number of clusters $N$. There are two NLO diagrams shown in Fig. 9 and their
formulae are given by

$\Gamma_{2}^{\mathrm{NLO}}=-\frac{2}{N_{s}}\sum_{\alpha,\beta}^{M}\sum_{i}^{N}\sum_{\lambda}^{N_{s}}\oint\oint
dt_{c}\,dt^{\prime}_{c}\,W^{11}_{\alpha i,\beta i}(t_{c},t^{\prime}_{c})\\\
\times\Big{(}G^{y,z}_{i\lambda}(t_{c},t^{\prime}_{c})G^{z,y}_{i\lambda}(t_{c},t^{\prime}_{c})-G^{y,y}_{i\lambda}(t_{c},t^{\prime}_{c})G^{z,z}_{i\lambda}(t_{c},t^{\prime}_{c})\Big{)}\\\
-ig^{2}\omega_{c}\sum_{\alpha}^{M}\sum_{i}^{N}\oint\oint
dt_{c}\,dt^{\prime}_{c}\,W^{22}_{\alpha i,\alpha
i}(t_{c},t^{\prime}_{c})D^{\phi\phi}_{\alpha\alpha}(t_{c},t^{\prime}_{c}).$
(96)

The rest of the terms in $\Gamma_{2}$ are either next-to-next-to-leading-order
(NNLO) in $1/N_{s}$ or scale sub-extensively with system size (Fig. 10). In
this work we neglect these terms and take

$\Gamma_{2}\approx\Gamma_{2}^{\mathrm{LO}}+\Gamma_{2}^{\mathrm{NLO}}.$ (97)

Figure 10: Some of the neglected diagrams in this work which appear after NLO
terms. The left diagram is sub-extensive in system size and the rest are NNLO
in $1/N_{s}$.

It is worth showing which diagrams we are keeping in terms of the original
action prior to HS transformation (Fig. 6). As mentioned before, the LO part
describes a retarded self-interaction of the spin expectation value expressed
in terms of fermion Green’s function (according to Eq. 49)

$\overline{\expectationvalue{\sigma^{x}}}=2G_{-+}^{z,y}(t,t),$ (98)

given by closed loops in the first term of Fig. 6, within the same cluster and
mediated by photons. In the regime of fast photons and at steady state, this
reduces to the ferromagnetic interaction in Eq. 6. For NLO terms, the two
diagrams in Fig. 9 are equivalent to the sum of an infinite number of diagrams
in the original representation of the theory, shown in the brackets of Fig. 6.
This can be verified by plugging the diagrammatic expression of $W$ given in
Fig. 7c into the diagrams of Fig. 9, and subsequently substitute the lines for
bare $W$ with simple dots, as the latter is just a constant. This infinite
series is a byproduct of disordered couplings. For the Dicke model without
disorder (see Appendix E), only the first NLO diagram would contribute since
in the absence of disorder, dashed lines would disappear. This makes the rest
of NLO diagrams 2-particle-_reducible_ and hence, forbidden in the expansion
of $\Gamma_{2}$.

###  IV.8 Equations of Motion

The equations of motion for Green’s functions and field expectation values are
obtained similarly to Eqs. (37) and (38) in the absence of sources as

$\displaystyle\frac{\delta\Gamma}{\delta\tilde{\chi}^{\sigma}_{\alpha
i}(t)}=0,\qquad\frac{\delta\Gamma}{\delta W^{\sigma\sigma^{\prime}}_{\alpha
i,\beta j}(t_{c},t^{\prime}_{c})}$ $\displaystyle=0,$ (99)
$\displaystyle\frac{\delta\Gamma}{\delta
D^{\rho\rho^{\prime}}_{\alpha,\beta}(t_{c},t^{\prime}_{c})}=0,\qquad\frac{\delta\Gamma}{\delta
G^{\alpha\beta}_{i\lambda,j\lambda^{\prime}}(t_{c},t^{\prime}_{c})}$
$\displaystyle=0.$ (100)

After taking functional derivatives of Eq. (86), we see that the equations of
motion for Green’s functions can always be cast compactly as

$\displaystyle\hat{G}^{-1}$ $\displaystyle=\hat{G}_{0}^{-1}-\hat{\Sigma},$
(101) $\displaystyle\hat{D}^{-1}$ $\displaystyle=\hat{D}_{0}^{-1}-\hat{\Pi},$
(102) $\displaystyle\hat{W}^{-1}$
$\displaystyle=\hat{W}_{0}^{-1}-\hat{\Omega},$ (103)

known as Dyson equations [62, 87, 101]. The matrices $\hat{\Sigma}$,
$\hat{\Pi}$ and $\hat{\Omega}$ are fermion, photon and Ising field self-
energies, respectively. They are given in terms of the functional derivatives
of $\Gamma_{2}$ as

$\displaystyle\Sigma^{\alpha,\beta}_{i\lambda,j\lambda^{\prime}}(t_{c},t^{\prime}_{c})$
$\displaystyle\equiv-2i\frac{\delta\Gamma_{2}}{\delta
G^{\beta,\alpha}_{j\lambda^{\prime},i\lambda}(t^{\prime}_{c},t_{c})},$ (104)
$\displaystyle\Pi^{\rho,\rho^{\prime}}_{\alpha,\beta}(t_{c},t^{\prime}_{c})$
$\displaystyle\equiv+2i\frac{\delta\Gamma_{2}}{\delta
D^{\rho^{\prime},\rho}_{\beta,\alpha}(t^{\prime}_{c},t_{c})},$ (105)
$\displaystyle\Omega^{\sigma,\sigma^{\prime}}_{\alpha i,\beta
j}(t_{c},t^{\prime}_{c})$
$\displaystyle\equiv+2i\frac{\delta\Gamma_{2}}{\delta
W^{\sigma^{\prime},\sigma}_{\beta j,\alpha i}(t^{\prime}_{c},t_{c})}.$ (106)

Due to their large size, the expanded forms of Eqs. (101-103), required for
numerically solving them, are given in Appendix B. We will only mention the
important details here. The fermion and photon Green’s functions will be
diagonal in the spin-site and photon-mode bases, respectively. Due to
permutation symmetry [19, 20, 178] we also have

$\displaystyle
G^{\alpha,\beta}_{i\lambda,j\lambda^{\prime}}(t_{c},t^{\prime}_{c})$
$\displaystyle=\delta_{ij}\delta_{\lambda\lambda^{\prime}}G^{\alpha,\beta}(t_{c},t^{\prime}_{c}),$
(107) $\displaystyle
D^{\rho,\rho^{\prime}}_{\alpha,\beta}(t_{c},t^{\prime}_{c})$
$\displaystyle=\delta_{\alpha\beta}D^{\rho,\rho^{\prime}}(t_{c},t^{\prime}_{c}).$
(108)

The same is true for their self-energies. The Ising field’s Green’s function
and its self-energy will acquire off-diagonal elements only in the photon
indices. Due to the emergent permutation symmetry after disorder averaging,
all diagonal elements of $W$ and $\Omega$ are the same. This holds also for
the off-diagonal elements of $W$ and $\Omega$:

$\displaystyle W^{\sigma,\sigma^{\prime}}_{\alpha i,\alpha
j}(t_{c},t^{\prime}_{c})$
$\displaystyle=\delta_{ij}V^{\sigma,\sigma^{\prime}}(t_{c},t^{\prime}_{c}),$
(109) $\displaystyle W^{\sigma,\sigma^{\prime}}_{\alpha i,\beta
j}(t_{c},t^{\prime}_{c})$
$\displaystyle=\delta_{ij}U^{\sigma,\sigma^{\prime}}(t_{c},t^{\prime}_{c}),\quad\alpha\neq\beta.$
(110)

A similar argument applies to $\tilde{\chi}$, and it has the same value for
all sites and photon modes:

$\tilde{\chi}^{\sigma}_{\alpha i}(t)=\tilde{\chi}^{\sigma}(t)$ (111)

We proceed similarly to the previous section. We find the equations of motion
at LO first and then consider the NLO corrections.

####  IV.8.1 Equations of motion at leading-order (LO)

For $\tilde{\chi}$, one always finds that it has the same value on forward and
backward Keldysh contours, as it should be since these are classical
variables. From Eqs. (95) and (99) we have

$\frac{\delta\Gamma^{\mathrm{LO}}}{\delta\tilde{\chi}^{1}}=0\to\tilde{\chi}^{2}(t)=-2\sqrt{\frac{N_{s}}{N+M}}G_{-+}^{y,z}(t,t),$
(112)
$\frac{\delta\Gamma^{\mathrm{LO}}}{\delta\tilde{\chi}^{2}}=0\to\tilde{\chi}^{1}(t)=\frac{2g^{2}\omega_{c}}{\sqrt{N+M}}\int_{0}^{t}dt^{\prime}\,D_{R}^{\phi\phi}(t,t^{\prime})\tilde{\chi}^{2}(t^{\prime}),$
(113)

where $D_{R}^{\phi\phi}$ is the photon response function defined as (for
details see Appendix B)

$D_{R}^{\phi\phi}(t,t^{\prime})\equiv-i\Theta(t-t^{\prime})\overline{\expectationvalue{\commutator{\hat{\phi}(t)}{\hat{\phi}(t^{\prime})}}},$
(114)

which describes the response of the photons to an external linear
perturbation. Note that the time variables used in Eq. (114) are normal
variables and are not defined on the time contour. Substituting (112) into
(113) gives

$\tilde{\chi}^{1}(t)=-\frac{4g^{2}\omega_{c}\sqrt{N_{s}}}{N+M}\int_{0}^{t}dt^{\prime}\,D_{R}^{\phi\phi}(t,t^{\prime})G_{-+}^{y,z}(t^{\prime},t^{\prime}),$
(115)

in agreement with the scaling relations for $\tilde{\chi}$ given in Eq. (91).
For $\Sigma$ we have at the leading order

$\Sigma_{\mathrm{LO}}^{\alpha,\beta}(t_{c},t^{\prime}_{c})=\parbox{34.14322pt}{\includegraphics[height=433.62pt]{Sigma_LO.pdf}}=-\frac{2iM}{\sqrt{N_{s}}}\epsilon_{x\alpha\beta}\delta(t_{c},t^{\prime}_{c})\tilde{\chi}^{1}(t_{c}).$
(116)

The temporally local elements of fermion self-energy describe the effective
magnetic field spins experience due to interactions. Accordingly, Eq. (116)
describes a time-dependent magnetic field $B_{x}$ in the $\hat{x}$-direction
experienced by each spin. This sheds light on the physical meaning of the
first component of the HS field: it is the magnetic field experienced by each
spin and its Green’s functions give the fluctuations of this field. This could
be inferred also from Eq. (73) which describes the coupling of
$\sigma^{x}=-2i\psi^{z}\psi^{y}$ to $\chi_{1}$. Similarly, Eqs. (112) and (98)
relate $\chi^{2}$ to $\expectationvalue{\sigma^{x}}$ as

$\overline{\expectationvalue{S^{x}_{i}(t)}}=-\sqrt{(N+M)N_{s}}\,\tilde{\chi}^{2}(t)/2.$
(117)

Furthermore, it can be shown that $W^{22}$ is related to the 2-point
correlation function of the large spins
$S_{i}^{x}=\sum_{\lambda}\sigma^{x}_{i\lambda}/2$ through

$\overline{\expectationvalue{S^{x}_{i}(t_{c})S^{x}_{i}(t^{\prime}_{c})}_{c}}=i(N+M)N_{s}\,V^{22}(t_{c},t^{\prime}_{c})/4.$
(118)

Eq. (118) can be rigorously proven by introducing source terms in the Keldysh
action which are coupled to $\sigma^{x}$, and then taking functional
derivatives with respect to the sources twice [87, 101]. By substituting
$\tilde{\chi}^{1}$ from Eq. (115) into Eq. (116), the effective magnetic field
$B_{x}$ is found to be

$B_{x}(t)=-\frac{4g^{2}\omega_{c}M}{N+M}\int_{0}^{t}dt^{\prime}\,D_{R}^{\phi\phi}(t,t^{\prime})G_{K}^{y,z}(t^{\prime},t^{\prime}).$
(119)

Interpreting Eq. (119) is straightforward now; spins perturb photons whose
displacement is given by $\phi$ which acts as an effective field applied to
spins through the Dicke coupling in Eq. (1), creating a self-interaction for
spins. Last, we investigate the photon sector at LO by calculating the photon
self-energy $\Pi$. Since all terms in $\Gamma_{2}$ only depend on the
$\phi\phi$ component of $D$, the only non-zero element of photon self-energy
is $\Pi^{\phi\phi}$ given by

$\Pi_{\mathrm{LO}}^{\phi\phi}(t_{c},t^{\prime}_{c})=\parbox{34.14322pt}{\includegraphics[height=433.62pt]{Pi_LO.pdf}}\\\
=-2iNg^{2}\omega_{c}\tilde{\chi}^{2}(t_{c})\tilde{\chi}^{2}(t^{\prime}_{c}).$
(120)

Written in terms of normal time variables, it is easy to show that
$\Pi^{\mathrm{LO}}$ does not alter the spectrum or equivalently [101, 11], the
response function of photons. Hence, its only effect is to increase the photon
population. At this order, spins pump photons, but without generating any
finite values for $\expectationvalue{\phi}$. As couplings to different
clusters have different signs, their MF contributions cancel out each other in
the thermodynamic limit and photon pumping is realized only at the level of
fluctuations. Furthermore, no changes of $D^{\phi\phi}_{R}$ at this order of
approximation means that the kernel of the effective interaction between spins
is given by its non-interacting form, i.e. it is the response function of a
damped (for $\kappa\neq 0$) harmonic oscillator. We note that the self-energy
of the Ising field vanishes at this order. We summarize the physics of the
problem in LO approximation. There is an effective Ising-like interaction of
spins within the same cluster with a retarded kernel given by the response of
cavity photons in the non-interacting limit. Photons are coherently pumped by
spins, but their energy levels and loss rates remain unaffected. At this order
of approximation, the model behaves very similar to the MF solution of the
Dicke model. The effective retarded spin-spin interaction is generated also
for the Dicke model, if we formally solve the equation of motion for the
photon mode in terms of $\sigma^{x}$ and then, substitute them back in the
equations of motion for spins.

####  IV.8.2 Equations of motion at next-to-leading-order

As is evident from Eq. (96), $\Gamma^{\mathrm{NLO}}_{2}$ does not depend on
$\tilde{\chi}$. Therefore, Eqs. (112) and (113) describe $\tilde{\chi}$ to
NLO. Accordingly, the general picture of an effective MF interaction between
spins remain unaltered. Although the interaction kernel given by
$D^{\phi\phi}_{R}$ will be renormalized by fluctuations at NLO.

The self-energies at NLO are found from the functional derivatives of
$\Gamma_{2}^{\mathrm{NLO}}$ and are given by

$\Sigma_{\mathrm{NLO}}^{\alpha,\beta}(t_{c},t^{\prime}_{c})=\parbox{34.14322pt}{\includegraphics[height=433.62pt]{Sigma_NLO.pdf}}\\\
=\frac{4iM}{N_{s}}\quantity(V^{11}(t_{c},t^{\prime}_{c})+(M-1)U^{11}(t_{c},t^{\prime}_{c}))\\\
\times\sum_{\gamma,\delta}\epsilon_{x\alpha\gamma}\epsilon_{x\beta\delta}G^{\gamma,\delta}(t_{c},t^{\prime}_{c}),$
(121)
$\Pi_{\mathrm{NLO}}^{\phi\phi}(t_{c},t^{\prime}_{c})=\parbox{34.14322pt}{\includegraphics[height=433.62pt]{Pi_NLO.pdf}}\\\
=2Ng^{2}\omega_{c}\,V^{22}(t_{c},t^{\prime}_{c}).$ (122)

Ising self-energies are non-vanishing at this order:

$(\Omega_{\mathrm{NLO}})_{\alpha\beta}^{11}(t_{c},t^{\prime}_{c})=\parbox{34.14322pt}{\includegraphics[height=433.62pt]{Omega11_NLO.pdf}}\hskip
25.0pt=4ig^{2}\omega_{c}\delta_{\alpha\beta}\\\
\times\Big{(}G^{yy}(t_{c},t^{\prime}_{c})G^{zz}(t_{c},t^{\prime}_{c})-G^{yz}(t_{c},t^{\prime}_{c})G^{zy}(t_{c},t^{\prime}_{c})\Big{)},$
(123)
$(\Omega_{\mathrm{NLO}})_{\alpha\beta}^{22}(t_{c},t^{\prime}_{c})=\parbox{34.14322pt}{\includegraphics[height=433.62pt]{Omega22_NLO.pdf}}\\\
=2g^{2}\omega_{c}\delta_{\alpha\beta}D^{\phi\phi}(t_{c},t^{\prime}_{c}).$
(124)

We see that the off-diagonal element of $W$ given by $U$ is multiplied by an
extra factor of $\mathcal{O}(M)$ in Eq. (121), and it has to be kept though it
is sub-leading compared to the diagonal element $V$. The resulting equations
of motion are a system of 36 coupled integro-differential equations for
different components of Green’s functions and $\tilde{\chi}$. The complete
expressions for these equations are given in Appendix B in terms of the
symmetric (Keldysh) and anti-symmetric (retarded and advanced) correlation
functions.

####  IV.8.3 Evaluation of glass order parameter

The formalism developed so far is sufficient to calculate some of the
correlation functions of our system which are usually the quantities of
interest for quantum dynamics. However, for systems with static disorder, we
can define new types of expectation values depending on the order of
calculating operator expectation values and disorder averaging [100]. As will
be explained in Sec. V.3.2, the quantity

$\displaystyle q(t)$
$\displaystyle\equiv\frac{1}{N_{s}^{2}}\quantity(\overline{\expectationvalue{S^{x}_{i}(t)}^{2}}-\overline{\expectationvalue{S^{x}_{i}(t)}}^{2}),$
(125)

is of particular importance in our model, and is named spin-glass order
parameter. $q(t)$ cannot be evaluated directly in terms of the “normal”
correlation functions, simply because of the non-commutativity of taking
expectation values and disorder averaging in Eq. (125). However, it can still
be calculated thanks to the versatility of the Keldysh approach in dealing
with quenched disorder [11, 101]. The quantum expectation value of an operator
$\hat{\mathcal{O}}(t)$ before disorder averaging can be written as

$\expectationvalue{\hat{\mathcal{O}}(t)}=\int\mathcal{D}\quantity[\phi]\mathcal{O}(t)e^{iS\quantity[\phi,g_{\alpha
i}]},$ (126)

and therefore

$\expectationvalue{\hat{\mathcal{O}}(t)}^{2}=\quantity(\int\mathcal{D}\quantity[\phi]\mathcal{O}(t)e^{iS\quantity[\phi,g_{\alpha
i}]})^{2}\\\
=\int\mathcal{D}\quantity[\phi_{1}]\mathcal{D}\quantity[\phi_{2}]\mathcal{O}_{1}(t)\mathcal{O}_{2}(t)e^{iS\quantity[\phi_{1},g_{\alpha
i}]+iS\quantity[\phi_{2},g_{\alpha i}]}.$ (127)

We simply have written $\expectationvalue{\mathcal{O}}^{2}$ in terms of a
Keldysh path integral with duplicated fields and $\phi$ is a shorthand
notation for all of the fields in the action. Note that up to this point
fields belonging to different copies, $\phi_{1}$ and $\phi_{2}$, do not
interact with each other. We can now straightforwardly find

$\overline{\expectationvalue{\hat{\mathcal{O}}(t)}^{2}}=\int\mathcal{D}\quantity[\phi_{1}]\mathcal{D}\quantity[\phi_{2}]\mathcal{O}_{1}(t)\mathcal{O}_{2}(t)\overline{e^{iS\quantity[\phi_{1},g_{\alpha
i}]+iS\quantity[\phi_{2},g_{\alpha i}]}}.$ (128)

Since the actions of both copies depend on the same realization of $g_{\alpha
i}$, averaging over disorder couples fields of different replicas, and
generates an effective interaction between them. Before applying Eq. (128) to
our problem, we make some remarks about our finding. Clearly, there is a
strong resemblance between our result and the replica trick [179, 100, 94, 11]
as they both involve more than one copy of the system. However, there are also
crucial differences between the two. In the replica trick the number of
replicas is taken to zero via analytical continuation while here we are
strictly working with 2 copies of the system. Naturally, for calculating
higher powers such as $\overline{\langle\mathcal{O}\rangle^{n}}$ we would need
more copies of the system. Furthermore, in the replica formalism there is the
possibility of replica symmetry breaking (RSB), where the correlation
functions of fields from different replicas do not vanish, with important
consequences for the physics of the problem [179, 100, 94, 11]. Depending on
the particular system and other parameters such as temperature, RSB may or may
not occur while the system is nevertheless a glass [100, 126]. In our
approach, the correlation functions of the fields from different copies can be
non-zero regardless of the system being a glass or not. However, as we will
explain later (cf. Sec. V.3.2), if 1-point functions are vanishing such that
Eq. (125) is finite, the system has glassy behavior. Despite these
differences, we call Eq. (128) the “replicated model” for convenience. We note
that a similar replica approach has been used by Ref. [98] to study
measurement-induced phase transitions due to continuous-time measurements. We
apply Eq. (128) to the HS transformed interaction part of the action in Eq.
(68). $S_{g\chi\phi}$ is the only term depending on $g_{\alpha i}$ and has to
be averaged in Eq. (128). The result of disorder averaging are three terms,
two of them couple fields from the same copies and correspond to the effective
vertex $S_{\chi\phi}$ in Eq. (76). The third term gives an interaction between
fields of different copies and reads as

$S^{12}_{\chi\phi}=2ig^{2}\omega_{c}\sum_{\alpha i}\iint
dt_{c}\,dt^{\prime}_{c}\,\chi^{2,1}_{\alpha
i}(t_{c})\phi^{1}_{\alpha}(t_{c})\chi^{2,2}_{\alpha
i}(t^{\prime}_{c})\phi^{2}_{\alpha}(t^{\prime}_{c}),$ (129)

where the extra index specifies the copy to which a field belongs. Using this
“replicated Keldysh field theory”, we can obtain $q(t)$ in terms of the
diagonal ($t=t^{\prime}$) elements of inter-replica correlation functions. To
distinguish inter-replica correlators from the normal ones, we represent the
former with a tilde mark below them such as
$\underaccent{\sim}{D}^{\rho,\rho^{\prime}}_{\alpha,\beta}(t_{c},t^{\prime}_{c})$.
For example, for $q$ we have

$q(t)=i\frac{N+M}{N_{s}}\underaccent{\sim}{V}^{22}(t,t).$ (130)

For the replicated theory, there will be 4 more independent equations of
motion that have to be solved in addition to those of the previous section.
Crucially, these extra equations do not alter the dynamics of replica-diagonal
quantities, as expected, since replicas are just abstractions and they are not
“aware” of each other. Non-replica diagonal quantities, on the other hand,
depend both on replica diagonal and non-replica diagonal correlation
functions. Furthermore, all of the non-replica diagonal response functions
turn out to vanish, as perturbing one replica cannot leave any effects on the
other one. We leave the details of the calculations and the extra equations of
motion to Appendix D.

##  V Results

In this section, we will report our findings in the following order. First, we
discuss in Sec. V.1 the results of approximating the effective action only to
LO. In Sec. V.2, we demonstrate that the NLO corrections significantly change
the dynamics for all values of $N_{s}$, motivating the necessity of keeping
NLO effects. In Sec. V.3, we discuss the formation of SG phase in the system
by studying various physical characteristics of SG in our system. We will also
address the effect of spin size $N_{s}$ and photon frequency on the glassy
behavior of the system. We will study the quench dynamics starting from a
polarized spin state specified by the angles $(\theta_{0},\varphi_{0})$ such
that

$\expectationvalue{\vec{\sigma}}_{0}=(\sin{\theta_{0}}\cos{\varphi_{0}},\sin{\theta_{0}}\sin{\varphi_{0}},\cos{\theta_{0}}),$
(131)

and the vacuum state for photons $\ket{0}$ satisfying $\hat{a}\ket{0}=0$. We
turn on the spin-photon coupling $g$ at $t=0$ and let the system evolve.

###  V.1 Results at LO

Figure 11: Dynamics at the leading-order of approximation describing the
motion of classical spins coupled via a retarded interaction, for
$g/g_{c}=1.27$ and $(\theta_{0},\varphi_{0})=(0.7\pi,0)$. Top panel: dynamics
of spins in the adiabatic limit $\omega_{c}/\Delta=5.0$ and for $\kappa=0$
(dark curve) and $\kappa/\Delta=0.5$ (light curve). Bottom panel: the same as
top panel only for slow photons $\omega_{c}=\Delta$, showing tunneling events
when photon loss is weak.

The equations of motion for the diagonal elements of fermion Green’s functions
become decoupled from the non-diagonal elements at LO. This allows us to write
the dynamics of magnetization in a transparent form as (Appendix D.1)

$\displaystyle\derivative{t}\expectationvalue{\sigma^{x}(t)}$
$\displaystyle=-\Delta\expectationvalue{\sigma^{y}(t)},$ (132)
$\displaystyle\derivative{t}\expectationvalue{\sigma^{y}(t)}$
$\displaystyle=\Delta\expectationvalue{\sigma^{x}(t)}$
$\displaystyle-J\omega_{c}^{2}\int_{0}^{t}\expectationvalue{\sigma^{z}(t)}D^{\phi\phi}_{R}(t,t^{\prime})\expectationvalue{\sigma^{x}(t)}\,dt^{\prime},$
(133) $\displaystyle\derivative{t}\expectationvalue{\sigma^{z}(t)}$
$\displaystyle=J\omega_{c}^{2}\int_{0}^{t}\expectationvalue{\sigma^{y}(t)}D^{\phi\phi}_{R}(t,t^{\prime})\expectationvalue{\sigma^{x}(t)}\,dt^{\prime},$
(134)

where $J$ is the coupling of LMG model defined in Eq. (4). The photon response
function $D^{\phi\phi}_{R}$ is given by its bare value at LO which is (the
minus of) the response function of a damped harmonic oscillator

$D^{\phi\phi}_{R}(t-t^{\prime})=-\frac{\Theta(t-t^{\prime})}{\omega_{c}^{2}}e^{-\kappa(t-t^{\prime})}\sin{\omega_{c}(t-t^{\prime})}.$
(135)

Figure 12: Evolution of photon population per each mode normalized by cluster
size $N_{s}$, for $g/g_{c}=1.27$ and $(\theta_{0},\varphi_{0})=(0.7\pi,0)$.
Top panel: the adiabatic limit $\omega_{c}/\Delta=5.0$ and for $\kappa=0$
(dark curve) and $\kappa/\Delta=0.5$ (light curve). Bottom panel: the same as
top panel only for slow photons $\omega_{c}=\Delta$.

Eqs (132)-(134) describe the motion of classical angular momentum variables
with a conserved vector length $\sqrt{\tilde{\sigma}\cdot\tilde{\sigma}}=1$.
Adiabatic elimination of photons [180] amounts to approximate the integral in
Eq. (134) as

$\approx\expectationvalue{\sigma^{y}(t)}\expectationvalue{\sigma^{x}(t)}\int_{0}^{t}D^{\phi\phi}_{R}(t-t^{\prime})\,dt^{\prime}\approx\expectationvalue{\sigma^{y}(t)}\expectationvalue{\sigma^{x}(t)}\\\
\times\int_{0}^{\infty}D^{\phi\phi}_{R}(t^{\prime})\,dt^{\prime}\approx-\frac{1}{\omega_{c}(\omega_{c}^{2}+\kappa^{2})}\expectationvalue{\sigma^{y}(t)}\expectationvalue{\sigma^{x}(t)}.$
(136)

Substituting (136) into Eqs. (133) and (134) results in the MF equations of
motion for LMG model with a coupling modified by photon loss according to
$J\to J\omega_{c}^{2}/(\omega_{c}^{2}+\kappa^{2})$. From a physical point of
view, the above derivation shows that at LO, our approximation maps each
cluster to a classical LMG system _without coupling to other clusters_.
Therefore, the LO approximation describes the FM to PM transition of an
infinite range Ising model with retarded interactions. The critical coupling
of this system is determined according to the condition $\Delta=J$ [116] to be

$g_{c}=\frac{1}{2}\sqrt{\frac{\Delta(\omega_{c}^{2}+\kappa^{2})}{\omega_{c}}\frac{N+M}{M}}.$
(137)

Although only valid in the adiabatic regime of the MF solution, we will use
$g_{c}$ throughout this paper and scale the coupling $g$ according to it when
comparing our results for different values of $\Delta$, $\omega_{c}$ or
$\kappa$. Clearly, dynamics at LO do not have any features unexplored in the
past, and we only report the results of our simulations as a consistency check
of our approach. In Fig. 11 we have shown the dynamics of spins initiated
close to an equilibrium state of Eqs. (132)-(134) inside the FM phase. We see
that for fast photons (top panel of Fig. 11), the system remains close to the
minimum with oscillations which are smoothed by photon loss. For slow photons
when adiabatic elimination does not work (bottom panel of Fig. 11), the
fluctuations induced by photons’ dynamics relieve FM correlations and can
create rare tunneling events. With photon loss, the destructive effect of slow
photons on FM order is reduced due to the effective damping that slows down
the spins towards the bottom of the nearby energy minimum.

Figure 13:  Comparison of spin dynamics at LO (dashed lines) and NLO (solid
lines) approximations for a large spin size ($N_{s}=10^{5}$), following a PM
quench ($g/g_{c}=0.71$) and in the limit of fast photons $\omega_{c}/\Delta=5$
with ($\kappa/\Delta=0.5$, light curves) and without (dark curves) photon
loss. Top, middle and bottom rows show the evolution of
$\expectationvalue{\sigma^{x}}$, $\expectationvalue{\sigma^{z}}$ and the
expectation value of spin size, respectively. Spins dephase to a single point
due to random static couplings in the system.

Having explored spin dynamics at LO, we now discuss the evolution of photon
sector. As we said before, the response of photons given by the retarded
function $D^{\phi\phi}_{R}$ is unaltered by interactions at LO. However, the
symmetric correlation function of photons and photon population are affected
by the coupling to spins. Spins pump cavity modes due to the presence of the
transverse field $\Delta$, and as shown in Fig. 12, create large oscillations
(if $\kappa=0$) or a superradiant burst (if $\kappa\neq 0$) in the population
of photons $n(t)$. We do not discuss dynamics of the PM phase in LO
approximation here due to their trivial nature, and will show them in
comparison with NLO results later. The insufficiency of the LO approximation
is clear by above observations. At this order of approximation, the size of
each cluster $N_{s}$ only appears as an overall scaling factor in photon
number $n(t)\propto N_{s}$, and spin dynamics do not depend on $N_{s}$. The
latter is not physically valid, as we expect fluctuations to affect spin
dynamics noticeably when the spin per each cluster
$\hat{S}_{i}=\sum_{\lambda}\hat{\sigma}_{i\lambda}/2$ becomes smaller as
$N_{s}$ is decreased. The other important shortcoming of LO approximation is
visible in the dynamics of photon population in the adiabatic limit without
loss, where the amplitude of the oscillations in $n(t)$ does not change and
remains constant (Fig. 12). In reality, photons experience dissipation even in
a perfect cavity, as they can be reabsorbed by spins over longer timescales
and the system is expected to equilibrate eventually. It is clear from the
discussion above that LO approximation fails to capture this effect. As we
will see next, NLO corrections take into account fluctuations and the
reabsorption of photons, in addition to non-trivial regimes of dynamics such
as glassy behavior.

###  V.2 Results at NLO for Large Spins

Below we will report on the dynamics at NLO for quenches to PM and FM phases
of the system, in the limit of large spin per cluster ($N_{s}=10^{5}$). We
note that we label these phases according to the behavior of the system at the
MF level. The dynamics may change significantly when going beyond MF theory,
making these labels inaccurate a posteriori.

####  V.2.1 Dynamics in the paramagnetic phase

Figure 14:  Comparison of spin dynamics at LO (dashed lines) and NLO (solid
lines) approximations for a large spin size ($N_{s}=10^{5}$), following a PM
quench ($g/g_{c}=0.71$) and in the limit of resonant photons
$\omega_{c}=\Delta$ with ($\kappa/\Delta=0.5$, light curves) and without (dark
curves) photon loss. Top, middle and bottom rows show the evolution of
$\expectationvalue{\sigma^{x}}$, $\expectationvalue{\sigma^{z}}$ and the
expectation value of spin size, respectively. Spins dephase to a single point
due to random static couplings in the system.  Figure 15:  Comparison of spin
dynamics at LO (dashed lines) and NLO (solid lines) approximations for a large
spin size ($N_{s}=10^{5}$), following a FM quench ($g/g_{c}=1.27$) and in the
limit of fast photons $\omega_{c}/\Delta=5$ with ($\kappa/\Delta=0.5$, light
curves) and without (dark curves) photon loss. Top panel: following a quick
collapse, the decay of FM order parameter strongly slows down and a prethermal
plateau is formed. The middle and bottom panels show the evolution of
$\expectationvalue{\sigma^{z}}$ and the radius of the Bloch sphere,
respectively.

We initialize the system in the ground-state of photons and the spin state
with $(\theta_{0},\varphi_{0})=(0.7\pi,0)$ in Eq. (131). We choose a coupling
below $g_{c}$ given in Eq. (137), and consider both cases of perfect and lossy
cavities. The results for fast photons ($\omega_{c}/\Delta=5$) are shown in
Fig. (13), in comparison with LO results. We see that the evolution of
$\expectationvalue{\sigma^{x}}$ is considerably altered by taking fluctuations
into account. Spin dynamics is now damped even in the absence of photon loss.
This is expected, as dissipation is a natural byproduct of interactions in a
many-body system. We also see that photon loss has a minor impact on spin
dynamics as it is weaker than the fluctuation-induced dissipation. As shown in
the middle panel of Fig. 13, $\sigma^{z}$ shows a surprising behavior at NLO
by not relaxing to its minimum value, in contrast to the LO result which at
least, when $\kappa\neq 0$, approaches $-1$. However, there is a clear
explanation for this phenomenon. The value of $g^{2}$ gives the variance of
the disordered coupling in the system. This means that, while most of the
individual couplings for each realization of the disorder are smaller than
$g$, some of them are still large enough to weaken the PM configuration of the
ground-state without causing a phase transition in the system. The fact that
this behavior is purely due to random interactions is supported by our
observation that the NLO approximation for non-disordered Dicke model results
in the decay of spins into a state with $\sigma^{z}=-1$ in the PM phase. We
also have shown the evolution of the spin vector size defined by

$R(t)\equiv\sqrt{\expectationvalue{\sigma^{x}}^{2}+\expectationvalue{\sigma^{y}}^{2}+\expectationvalue{\sigma^{z}}^{2}},$
(138)

in Fig. 13. While $R$ is a constant of motion of Eqs. (132)-(134), it changes
when fluctuations are considered. The final state of the system is a PM with a
smaller spin size. Note that $R(t\to\infty)$ is independent of the initial
state and is always smaller than one even at the lowest temperatures, only due
to the frustrated nature of the system. In the following discussion, $R(t)$
will be a useful proxy for assessing the impact of correlations in dynamics.
It is the radius of the Bloch sphere, given at $t=0$ by the spin coherent
state in Eq. (131), and a constant of motion for all-to-all interacting spin
systems with homogeneous couplings and collective dissipation [119, 181, 19].
It stays constant over time because in this class of systems, the MF
approximation is exact in the thermodynamic limit and therefore no higher
order cumulants are formed, which would make $R(t)$ shrink, signaling the
onset of a strongly correlated regime. Spin dynamics at NLO approximation for
the case of resonant photons $\omega_{c}=\Delta$ are illustrated in Fig. 14.
As expected, dynamics become more irregular both at LO and NLO and photon loss
has a more dramatic effect on the dynamics.

####  V.2.2 Dynamics in the ferromagnetic phase

Figure 16: Evolution of spin vectors projected to the $x\hat{O}z$ plane for
large cluster size $N_{s}=10^{5}$, in the adiabatic limit
$\omega_{c}/\Delta=5$ (top panel) and resonant $\omega_{c}=\Delta$ (bottom
panel) limits, for $\kappa/\Delta=0$ (dark curves) and $\kappa/\Delta=0.5$
(light curves). In the adiabatic limit, spins feature a spiral relaxation
around the axis of FM ordering (dashed line) and towards the origin (parallel
to the arrow).

We quench the coupling to $g=1.27g_{c}$ where a FM state is realized at LO as
shown previously,. The initial spin vector is again taken to be
$(\theta_{0},\varphi_{0})=(0.7\pi,0)$, such that spins are close to the
ground-state of the classical model in Eqs. (132)-(134) at the chosen coupling
strength. Similar to the PM case, we take a large spin size ($N_{s}=10^{5}$)
where the effect of NLO corrections is supposed to be small. It will become
clear that this is an incorrect assumption and NLO contributions are
significant. The results of the numerics for the adiabatic limit
$\omega_{c}/\Delta=5$ are depicted in Fig. 15 at LO and NLO approximations. At
LO, $\sigma^{x}$ and $\sigma^{z}$ show small oscillations around the
equilibrium of the MF dynamics and the spin vector is confined to the surface
of the Bloch sphere ($R=1$). Fluctuations captured at NLO drastically alter
the dynamics. All components of spin decay and the spin vector shrinks toward
the center of the Bloch sphere. The spin decay features an interesting
profile. For an initial period, the relaxation is quick and spins experience a
collapse to a smaller but finite value. Following this, the relaxation becomes
very slow and the spin vector spirals around the axis of FM order of the LO
solution (Fig. 16). This situation is similar to the phenomenon of
prethermalization in the quench dynamics of many-body quantum systems [67, 24,
79, 182, 183, 184, 185, 186, 187, 188, 189], although here the system can be
open. The prethermal behavior is also seen in the time evolution of photon
number $n(t)$ shown in Fig. 17a. During the prethermal plateau of spins,
$n(t)$ has a nearly stationary value with a slow growth (visible in the
overall slope of the light curve in Fig. 17b) towards the true equilibrium
state. Photon losses only qualitatively affect spin dynamics, by weakly
accelerating the process of relaxation.

(a)

(b)

Figure 17: (a) Photon population for a quench into the ordered phase with
$g/g_{c}=1.27$, $\omega_{c}/\Delta=5$ and $\kappa=0$ at LO (dark curve) and
NLO (light curve). (b) Photon number per mode for $g/g_{c}=1.27$ and for
$N=1000,\,200,\,20,\,5$ (light to dark curves). Relaxation of late-time
oscillations is faster for smaller $N_{s}$.

Spin dynamics in the regime of resonant photons ($\omega_{c}=\Delta$) is shown
in Fig. 18. We see that the relaxation is faster in this regime compared to
the adiabatic regime as expected, since photon radiation is kinematically more
favorable in this limit. Also, the prethermal plateau does not appear any
longer due to the enhancement of relaxation. The spiral evolution of spin
vector changes to an irregular form (right panel in Fig. 16), but still
approaches the origin at long times.
As we emphasized before, labeling the system as a FM for $g>g_{c}$ is only
valid at the MF level. It is clear from the previous section that the system
is not truly a FM as the magnetization decays with time. In the upcoming
sections, we will provide evidence that the system is, in fact, a spin glass
(SG) in this regime. However, it is not possible to uniquely determine a SG by
only looking at single spin observables. The only indirect evidence for SG
that we have so far is the slow relaxation of magnetization. We again remark
that our results are expected to be correct for $\eta\gtrsim 1$ where the
model in Eq. (1) is expected to host a SG [96, 123].

Figure 18: Spin dynamics at LO (dashed lines) and NLO (solid lines)
approximations for $N_{s}=10^{5}$, following a FM quench ($g/g_{c}=1.27$) and
in the limit of dynamically active photons $\omega_{c}=\Delta$ with
($\kappa/\Delta=0.5$, light curves) and without (dark curves) photon loss.

####  V.2.3 Dependence of the spin dynamics on cluster size

The number of spins per cluster $N_{s}$ is not only a control parameter for
the expansion of the effective action, but it also determines the strength of
fluctuations in the system. This assumption stems from the conventional
understanding that larger spins often exhibit a more classical behavior
compared to small spins [24]. Here, the size of the spin for each cluster is
determined by the number of spins inside each cluster $N_{s}$. Our solution of
the problem also supports the presence of stronger fluctuations for smaller
values of spins per cluster. This can be seen, for example, in Eq. (121) where
the fermion self-energy scales as $N_{s}^{-1}$. The dynamics of spins at NLO
are shown in fig. 19 for different cluster sizes. The main effect of lowering
$N_{s}$ or equivalently, increasing fluctuations is a faster relaxation of
spins. As $N_{s}$ is decreased, the transient FM state lasts shorter. For
sufficiently small $N_{s}$, the prethermal state is completely bypassed and
the system experiences a quick relaxation towards equilibrium (Fig. 20).
Similarly, photon number relaxes faster for smaller $N_{s}$, while for larger
cluster sizes it has weakly damped oscillations (Fig. 17b).

###  V.3 Dynamics of spin glass

The results of NLO approximation given in previous sections indicate that the
strong coupling limit of the model in Eq. (1) is not a FM and corrections due
to frustration and fluctuations drastically change the MF phase diagram. From
the $N_{s}=1$ limit of the problem [123, 96, 97], we know that the system
hosts a SG at sufficiently low temperatures and strong couplings. The main
questions are the following: (i) What are the physical signatures of SG in and
out of equilibrium? (ii) Can our formalism capture the far from equilibrium
dynamics of SG? (iii) How does the SG phase behave for different atomic
ensemble sizes? In the upcoming sections, we will discuss two direct measures
of SG phase. First, we consider the disorder average of the square of local
magnetization of each spin ensemble which gives a proxy of the presence of
frozen spin configurations. The second one is related to the ability of the
system to retain its memory of the far past and is connected to the phenomenon
of aging, encoded in 2-point correlation functions of the system for large
time separations. We will demonstrate that 2PI can access both of these
measures and is able to track their evolution in real time. At the end, we
will discuss how the value of $N_{s}$ affects the SG phase.

Figure 19: Spin dynamics for different numbers of cluster sizes $N_{s}$,
following a FM quench ($g/g_{c}=1.27$) and in the limit of fast photons
$\omega_{c}/\Delta=5$ without losses. Relaxation is faster for smaller cluster
sizes due to enhanced fluctuations. The prethermal FM state disappears sooner
as $N_{s}$ is decreased.

####  V.3.1 Characteristics of spin glass

According to Landau’s theory of phase transitions (PT), different phases of
matter can be classified according to their symmetries [190, 170, 191]. A
phase transition corresponds to a change in the symmetry group of a physical
system. For example, the PM to FM transition of an Ising FM corresponds to the
spontaneous breaking of the $Z_{2}$ symmetry of the Hamiltonian. The
occurrence of the PT is signaled by the continuous growth of a parameter which
vanishes on one side of the transition. Such parameter is termed as the order
parameter and the way it grows with the control parameter of the transition
provides us with important information about the nature of the PT [191, 170,
11]. While the PM to SG phase transition fits into the symmetry breaking
paradigm, the proper definition of an order parameter for SG has been a
subject of debate for decades, with various candidates [100, 94, 99]. In the
following, we will consider two of these proposed order parameters to study
the formation of SG in our system.

####  V.3.2 Local magnetization

Figure 20: Comparison of spin vectors projected to the $xz$ plane for small
and large cluster sizes, in adiabatic $\omega_{c}/\Delta=5$ and for
$\kappa/\Delta=0$. The FM spiral of large spin ensembles crosses over to a
quicker smooth relaxation for smaller values of $N_{s}$.

(a)

(b)

Figure 21: Formation of SG order for different values of system-bath coupling
$\kappa$. (a) For $\theta_{0}=0.5\pi$, the energy density after the quench is
large and the final temperature is too high to realize a SG, unless the system
is allowed to cool down by coupling it to the bath. (b) For
$\theta_{0}=0.7\pi$, energy density is smaller and the system enters a weak SG
state without coupling to the bath, but it saturates to a smaller value of
$q$, and hence, glass order is weaker in this case. The other parameters are
$g/g_{c}=1.27$, $\omega_{c}/\Delta=5.0$, $N_{s}=5$ and $\eta=1$.

It is straightforward to define an order parameter for a FM. For instance, the
quantum Ising model given by

$H_{\mathrm{Ising}}=-\sum_{\expectationvalue{ij}}\sigma^{x}_{i}\sigma^{x}_{j}-g\sum_{i}\sigma^{z}_{i}-h\sum_{i}\sigma^{x}_{i},$
(139)

is in FM phase if the thermal average $\expectationvalue{\sigma^{x}_{i}}$ is
(i) finite and (ii) has the same sign for all $i$ if we take the thermodynamic
limit and then, take $h\to 0$. One can also fix $h=0$ and start from an
initial state with broken symmetry such that
$\expectationvalue{\sigma^{x}}_{0}\neq 0$ and let the system evolve towards a
FM equilibrium state. However, the symmetry will remain broken only in the
thermodynamic limit. On the other hand, the PM phase is characterized by
$\expectationvalue{\sigma^{x}_{i}}=0$ for all sites. At face value, a SG
corresponds to a situation where the $Z_{2}$ symmetry of the Hamiltonian is
broken such that $\expectationvalue{\sigma^{x}_{i}}$, _but its value and sign
are different for each site $i$_, due to disordered interactions. In the
thermodynamic limit, the average of
$m_{i}\equiv\expectationvalue{\sigma_{i}^{x}}$ over all sites is zero

$\lim_{N\to\infty}\frac{1}{N}\sum_{i}^{N}m_{i}=0.$ (140)

Since the total magnetization is zero both in SG and PM phases, it cannot
serve as an order parameter. However, we can define another quantity [100]

$q=\lim_{N\to\infty}\frac{1}{N}\sum_{i}^{N}m_{i}^{2},$ (141)

which is by definition finite for SG while still vanishing in the PM state. So
far we assumed that we are working with a single realization of the disorder
in the system. If we also average over the disorder, then Eq. (140) can be
replaced with

$\overline{m_{i}}=0,$ (142)

while instead of (141) we can use

$q=\overline{m_{i}^{2}}-\overline{m_{i}}^{2},$ (143)

where following the notation of Section IV, the overline represents disorder
averaging. We also have subtracted $\overline{m_{i}}^{2}$ which is vanishing
at equilibrium, but can be finite when looking at non-equilibrium dynamics. We
explained how to extract $q$ using the 2PI formalism in Sec. IV.8.3. In the
following we will show the dynamics of $q(t)$ as the light-matter coupling $g$
and the photon loss $\kappa$ are varied.

Figure 22:  Dynamics of SG order parameter (top panel) and global
magnetization (bottom panel) for different values of coupling strength.
$g_{c}$ is the critical coupling in the MF theory. This plot shows that the
observed value of $g_{c}$ is smaller than the one obtained at LO
approximation. The other parameters are $\kappa/\Delta=0.5$,
$\omega_{c}/\Delta=5.0$, $N_{s}=5$ and $\eta=1$.

(a)

(b)

Figure 23: (a) Symmetric correlation function $C$ at different waiting times
for quenches from an initial state with high energy density
$\theta_{0}=\pi/2$. For $\kappa=0$ (dashed lines), the system heats up to PM
state and $C$ quickly decays to zero for $t\to\infty$ (lines are almost on top
of each other). For $\kappa/\Delta=0.5$ (solid lines) the system cools down to
SG where $C$ remains finite for $t\to\infty$. (b) $C$ for an initial state
with low energy density $\theta_{0}=0.7\pi$. Regardless of $\kappa$ the system
ends up in the glass phase, however, SG order is stronger for
$\kappa/\Delta=0.5$ (solid lines) where the system is cooled down, compared to
$\kappa=0$ (dashed lines). For both figures $\omega_{c}/\Delta=5.0$,
$g/g_{c}=1.27$, $\eta=1$ and $N_{s}=5$.

We initialize the system in a state with $\expectationvalue{\sigma^{x}_{i}}=1$
and quench $g$ from zero to a sufficiently large value. We monitor the
dynamics of $q$ for different values of system-bath coupling from $\kappa=0$
(closed system) to $\kappa=0.5\Delta$. As can be seen in Fig. 21a, after some
fluctuations caused by the quench, $q$ approaches zero at long times when the
system is not coupled to the bath, while it grows to a finite value with a
faster speed as $\kappa$ is increased. An explanation for this behavior is
that the initial energy $E_{i}$ of the system with respect to the post-quench
Hamiltonian is large enough to put the system in a high temperature
equilibrium PM state if the system is closed. By coupling the system to the
bath, we allow the system to lose its excessive energy by leaking the photons
generated by the quench out of the cavity. Hence, the system cools down and
enters the SG phase with finite $q$.

(a)  (b)

Figure 24: Integrated response $\chi$ versus $C$ for different waiting times
$\tau$. (a) In PM phase, the slope is constant and is proportional to the
inverse of temperature. (b) In SG state, dynamics have multiple time scales
and the slope changes with $t$.

The explanation given above for the role of photon loss in the dynamics is
supported by the behavior of quenches with smaller $E_{i}$. We note that
$E_{i}$ can be calculated easily for the initial states considered in this
work:

$E_{i}=\frac{1}{2}NN_{s}\Delta\cos\theta_{0}.$ (144)

since we have $g=0$ in the pre-quench hamiltonian. By initializing spins
closer to $\theta_{0}=\pi$, we can reduce the energy density and attain a
lower temperature final state for isolated quenches. For instance, for
$\theta_{0}=0.7\pi$ the system enters the SG phase even without losses (Fig.
21b), although $q$ is smaller in this case due to the absence of cooling. We
consider the effect of the coupling strength $g$ on glassy dynamics. Starting
from $\theta_{0}=0.7\pi$ while connecting the system to an external bath by
switching on the photon loss, we look at $q$ for quenches to different values
of $g$ as shown in the top panel of Fig. 22. For weak couplings, we see that
$q(t)$ approaches zero at long times, while showing small fluctuations at
intermediate stages. We note that negative values for $q$ are not unphysical
when $\expectationvalue{\sigma^{x}}$ is finite according to Eq. (143). When
$g$ becomes large enough, $q$ saturates at a finite value at long times. We
see in the bottom panel of Fig. 22 that spin dynamics changes from underdamped
to overdamped at the same coupling where $q$ becomes finite. In Fig. 22 the
coupling $g$ is scaled with $g_{c}$ which was the critical coupling of the MF
limit in Eq. (137). A finite value of $q$ for $g/g_{c}<1$ indicates that NLO
contributions shift the boundary of PM and SG phases towards SG.

####  V.3.3 Temporal correlations and aging phenomena

As discussed before, the diagnosis of SG phase requires the evidence of a
frozen spin configuration which breaks the symmetry of the Hamiltonian and at
the same time, has no global magnetization. In the previous section we
discussed (the square of) the instantaneous magnetization for each site and
showed that it is finite in the glass phase while magnetization relaxes to
zero. The latter feature distinguished SG from FM. Another way to confirm a
frozen spin state in the system is to look at the temporal correlation between
the magnetization of each site at large time separations. For example, we take
a cluster and its total spin operator $S^{x}$ at two different times $t$ and
$t^{\prime}$. If spins are frozen, they do not fluctuate strongly in time and
the overlap of magnetization for the same site at two different times is
large. In fact, it remains finite even if
$\absolutevalue{t-t^{\prime}}\to\infty$. For classical degrees of freedom one
has

$q_{\mathrm{EA}}\equiv\lim_{\absolutevalue{t-t^{\prime}}\to\infty}\overline{\expectationvalue{S^{x}(t)S^{x}(t^{\prime})}}>0,$
(145)

which is the Edwards-Anderson (EA) order parameter [99, 100]. It is clear that
the above condition can also imply FM order in the system. In order to exclude
FM, the disorder averaged magnetization $\overline{\expectationvalue{S^{x}}}$
should also vanish. For the case of quantum operators in our problem, the
correlator in Eq. (145) can be generalized to

$C(\tau,\tau-t)\equiv\frac{1}{N_{s}^{2}}\sum_{\lambda,\lambda^{\prime}}^{N_{s}}\overline{\expectationvalue{\anticommutator{\sigma_{i\lambda}^{x}(\tau)}{\sigma_{i\lambda^{\prime}}^{x}(\tau-t)}}}_{c}.$
(146)

We note that the off-diagonal elements $\lambda\neq\lambda^{\prime}$ in (146)
can be non-zero. In the following, we consider the behavior of $C$ for the
same quenches discussed in the previous section. It will be shown that $C$
changes its behavior at the same points where $q$ in Eq. (125) becomes finite.
Therefore, both $C$ and $q$ are consistent measures of SG order in the system.
We again look at the effect of photon loss on quenches from hot (with
$\theta_{0}=0.5\pi$) and cold (with $\theta_{0}=0.7\pi$) initial states of the
system, according to Eq. (144). We consider $C(\tau,\tau-t)$ for different
values of the “waiting time” $\tau$ after the quench as a function of $t$. In
Fig. 23a, we have shown $C$ for $\theta_{0}=\pi/2$ and for zero and finite
$\kappa$. We see that $C$ behaves differently depending on whether the system
is isolated or not. When the system is closed (dashed lines in Fig. 23a), it
ends up in a high temperature PM state with correlations that relax quickly to
equilibrium and show weak dependence on the waiting time after the quench
$\tau$. The correlations also decay quickly with $t$, meaning that local
magnetization loses the memory of its past quickly in the PM phase. When the
system is allowed to cool down by emitting photons to the outside of the
cavity, a SG phase emerges (solid lines in Fig. 23a). $C$ becomes strongly
dependent on the waiting time $\tau$ and there is a crossover in the behavior
of $C$ for $t\ll\tau$ and $t\approx\tau$. The latter behavior is termed aging
[192, 193] and is another feature of quenches into SG phases. In aging, the
decay of $C(\tau,\tau-t)$ with $t$ becomes slower for “older” systems with
larger $\tau$. In this picture, there is a plateau in $C$ for $t\ll\tau$ (Fig.
23a) whose height gives the EA order parameter. The plateau is more visible in
figures shown in upcoming sections.

Figure 25: Dependence of $C$ on coupling strength for $g/g_{c}=1.27$ (solid
lines) where aging and SG are realized, and for $g/g_{c}=0.57$ (dashed lines)
with a PM profile. The other parameters are $\kappa/\Delta=0.5$,
$\omega_{c}/\Delta=5.0$, $\theta_{0}=\pi$, $\eta=1$ and $N_{s}=5$.

Aging is closely related to the lack of thermalization in spin glasses [192,
194, 102, 195]. The system is unable to reach equilibrium because of its
highly rugged energy landscape in which adjacent energy configurations,
connected by local spin flips, are separated by large energy barriers. As a
result, ergodicity is broken and system will not thermalize [194]. The
relation between aging and thermalization can be made more explicit by using
the fluctuation dissipation theorem (FDT) and fluctuation dissipation ratio
(FDR) defined below. At thermal equilibrium, the symmetric correlation
function $C(t,t^{\prime})$ and the response function $R(t,t^{\prime})$ of an
observable $\mathcal{O}$ are connected to each other via FDT [101, 87, 11]:

$R(t-t^{\prime})=-i\int\frac{d\omega}{2\pi}e^{-i\omega(t-t^{\prime})}\tanh(\frac{\omega}{2T})\tilde{C}(\omega),$