††thanks<EMAIL_ADDRESS>
# Condensation in hybrid superconducting cavity-microscopic spins systems with
finite-bandwidth drive

R. Au-Yeung  Advanced Technology Institute, University of Surrey, Guildford
GU2 7XH, United Kingdom  M. H. Szymańska  Department of Physics & Astronomy,
University College London, London WC1E 6BT, United Kingdom  E. Ginossar
Advanced Technology Institute, University of Surrey, Guildford GU2 7XH, United
Kingdom

###### Abstract

Using Keldysh field theory, we find conditions for non-equilibrium
condensation in the open Tavis-Cummings model under a direct finite-bandwidth
incoherent cavity drive. Experimentally, we expect the condensation transition
to be easily accessible to hybrid superconducting systems coupled to
microscopic spins, as well as to many other incoherently driven light-matter
systems. In our theoretical analysis, we explicitly incorporate the drive’s
spectral distribution into the saddle-point description. We show that the
injected incoherent photons create a drive-dependent effective coupling
between spin-1/2 particles. The condensation transition arises at a critical
regime of driving which we can now accurately predict. Our results also
provide important guidelines for future quantum simulation experiments of non-
equilibrium phases with hybrid devices.

##  I Introduction

The last few decades have seen enormous advances in experimental realization
and theoretical understanding of quantum condensation in a variety of physical
systems. They range from condensed matter models like microcavity exciton-
polariton systems Kasprzak2006 ; Szymanska2007 , electron-hole plasmas of
highly excited semiconductors Kremp2008 , and quantum magnets Zapf2014 ; to
circuit quantum electrodynamics (CQED) devices that demonstrate photon
condensation Marcos2012 , the Dicke states Wu2017 , and dressed collective
qubit states Fink2009 .

Among these, condensation in light-matter cavity QED systems offers a path
into deeper understanding of non-equilibrium phases of matter. Here we show
that this phenomenon could be studied with an analogue quantum simulation
Cirac2012 ; Xiang2013 ; Georgescu2014 ; Zagoskin2018 , which offers a highly
promising avenue for testing quantitative predictions of non-equilibrium
quantum theory Altman2021 . Analogue simulators can mimic the dynamics of
many-body quantum optical systems by reconstructing their Hamiltonian under
precisely controlled conditions Kurizki2015 ; Fitzpatrick2017 ; Blais2021 . In
particular, we consider the much studied hybrid quantum system consisting of a
superconducting resonator coupled to solid state spin centers Bonizzoni2018 ;
Mi2018 ; Clerk2020  as the ideal quantum simulator. Experimental realizations
can include a superconducting resonator coupled to NV- spin centers or other
impurities Schuster2010 ; Kubo2010 ; Amsuss2014 ; Liu2016 . Programmable
signal generators and sensitive amplifiers provide highly controllable driving
and probing in the microwave regime, even to the extent that quantum state
tomography can reconstruct the cavity states of microwave resonators
Kirchmair2013 ; Weiss2019 . Although CQED setups have well-understood
Hamiltonians, the question of how to achieve condensation with these systems
has not been addressed before. A key element is energy pumping (driving) which
we focus on here. We analyse this system in the presence of an incoherent,
finite-bandwidth, typically microwave frequency drive with a general spectral
distribution.

The non-equilibrium response to incoherent drive is challenging to analyze
because the drive induces interactions between the two-level atoms. This is
considered an open problem which was previously approached with an effective
phenomenological (or Markovian) description in the system’s equations of
motion Sieberer2016  or, where appropriate, with the introduction of fermionic
driving Szymanska2007 . In the latter it was shown that the system can access
both condensation and lasing regimes Kasprzak2006 ; Deng2010 ; Byrnes2014 ;
Deveaud2015  which are connected by a smooth crossover Keeling2005 ;
Keeling2010 ; Szymanska2013 .

In contrast, here we analyze the condensation transition based on a model with
a direct photonic drive. We use Keldysh field theory techniques Kamenev2011 ;
Sieberer2016  to study the consequences of a bosonic drive and obtain the
condensation phase diagram with directly applicable predictions. Technically,
our analysis is based on the saddle-point solution of the Keldysh action which
is standard procedure for studying condensation phenomena Sieberer2016 ;
Altland2010 . Note also that because the drive couples to the cavity photon
modes, population inversion of the two-level atoms cannot occur. Hence it is
physically impossible for the system to behave as a laser.

One of the central aims of our work is to determine the parameter space and
photon drive’s distribution function where the system can access the
condensation regime. We hypothesise that driving the cavity with Markovian
noise (flat spectral density) cannot reach the phase transition because it
strongly decoheres and thermalizes the system. Indeed, we see that the
threshold for condensation is only realistically achievable for a finite
bandwidth drive, a case which can be analyzed in our formalism.

In this work, we use the term “non-Markovian” driving to mean frequency-
dependent driving. We only observe spontaneous condensation when
preferentially driving the low frequency part of the spectrum. In the regime
of condensation, we observe that the condensation threshold is sensitive to
the width of the spectral profile of the drive. Larger widths require higher
drive strengths. Hence, extremely large spectral profiles are experimentally
unpractical. In fact, it was established in earlier work Szymanska2001  that
condensation cannot be achieved within the Markov approximation, which drives
the high and low energy modes with equal strength. The drive would cause
strong decoherence that impedes the formation of the condensate Szymanska2001
. Hence we cannot analyze the problem using a Markovian master equation, which
is a more appropriate method for quantum optical problems where energy
dependence of decay or drive are irrelevant.

##  II System and model

In this work, we use Keldysh field theory Kamenev2011 ; Sieberer2016  to
predict the onset of spontaneous coherence (via non-equilibrium condensation)
in many-body light-matter systems under direct incoherent driving with finite
bandwidth and amplitude. Standard quantum optics techniques cannot easily
treat frequency dependent (non-Markovian) drive. However the field theoretical
approach offers a natural language for including general spectral densities
and exploring how this affects the steady-state of the system.

We focus on the open Tavis-Cummings model (Fig. 1) featuring an ensemble of
spin-1/2 particles (or pseudo-spins representing two-level atoms) which are
dipole-coupled to one quantized cavity mode. The cavity is coupled to an
incoherent photonic drive. Unlike previous Keldysh-based theoretical studies
on microcavity polariton systems Szymanska2007 , we model the photonic drive
as a bosonic reservoir (bath) which feeds incoherent photons into the cavity.
We derive an extension to the self-energy of the spins which arise from
fluctuations in the cavity field caused by the incoherent drive. In addition,
the theory includes a low temperature reservoir to represent the non-radiative
environment.

Figure 1: Sketch of possible experimental setup uses NV- vacancy centers
(spin-1/2 particles) in diamond coupled to superconducting resonator
interacting with different types of environment. Frequency profile of photonic
drive characterized by photon mode number $n_{B}(\omega)$.

Note that this photonic driving generates a non-equilibrium condensation which
differs from the lasing and Dicke-type phase transitions Kirton2018 ;
Kirton2019 . This work cannot be directly compared with other phenomena such
as lasing and super-radiance. It is not an incremental study, but rather a new
regime for condensation in the driven-dissipative Tavis-Cummings model which
is directly relevant to CQED devices.

The non-equilibrium condensation is characterized by a non-thermal
distribution function. The structured spectral density of the drive and its
statistical properties are described by a customized distribution function.
The signatures of condensation are a sudden appearance of collective phase
coherence in the spins and condensation of the cavity photons in a coherent
state (Fig. 1).

We define the Hamiltonian as
$\hat{H}=\hat{H}_{\textrm{sys}}+\hat{H}_{\textrm{sys,bath}}+\hat{H}_{\mathrm{bath}}$
($\hbar=1$) where

$\displaystyle\hat{H}_{\textrm{sys}}=\sum_{\alpha}\left[\epsilon_{0}(b_{\alpha}^{\dagger}b_{\alpha}-a_{\alpha}^{\dagger}a_{\alpha})+g_{\alpha}(\psi
b_{\alpha}^{\dagger}a_{\alpha}+\mathrm{H.c.})\right]$
$\displaystyle+\omega_{0}\psi^{\dagger}\psi.$ (1)

Spin-up and -down states, represented using fermionic annihilation operators
$b_{\alpha}$ and $a_{\alpha}$ respectively, are dipole-coupled to cavity
photon mode $\psi$ with individual spin-photon coupling is $g_{\alpha}$. We
use a fermionic representation of spin-1/2 particles to facilitate a
straightforward use of the field theoretical formulation (although the action
can be formulated in terms of coherent spin states Nagaosa1995 , it is
technically difficult to treat). For example the operation
$b_{\alpha}^{\dagger}a_{\alpha}$ creates a spin-up particle with energy
$\epsilon_{0}$. The single-occupancy constraint

$b_{\alpha}^{\dagger}b_{\alpha}+a_{\alpha}^{\dagger}a_{\alpha}=1$ (2)

allows the fermionic representation to properly describe a spin-1/2 particle.
This eliminates the unphysical states $|\textrm{vacuum}\rangle$ and
$b_{\alpha}^{\dagger}a_{\alpha}^{\dagger}|\textrm{vacuum}\rangle$, leaving two
possible transitions for the ground
$|0\rangle=a_{\alpha}^{\dagger}|\textrm{vacuum}\rangle=a_{\alpha}^{\dagger}b_{\alpha}|1\rangle$
and excited state
$|1\rangle=b_{\alpha}^{\dagger}|\textrm{vacuum}\rangle=b_{\alpha}^{\dagger}a_{\alpha}|0\rangle$.
The constraint is maintained by considering only states which obey this
condition and bath distribution functions as described below. For simplicity,
the spin-photon coupling is on resonance i.e. the cavity and atomic
frequencies are both $\omega_{0}$, and there is no detuning in the system i.e.
$\omega_{0}-2\epsilon_{0}=0$. The system is coupled to two baths (reservoirs),

$\displaystyle\hat{H}_{\textrm{sys,bath}}=\sum_{\alpha,k}[\Gamma_{\alpha,k}(a_{\alpha}^{\dagger}A_{k}+b_{\alpha}^{\dagger}B_{k}+\mathrm{H.c.})$
$\displaystyle+\zeta_{k}(\psi^{\dagger}\Psi_{k}+\mathrm{H.c.})]$ (3)

given by
$\hat{H}_{\mathrm{bath}}=\sum_{k}\omega_{k}^{\Gamma}(A_{k}^{\dagger}A_{k}+B_{k}^{\dagger}B_{k})+\sum_{k}\omega_{k}^{\zeta}\Psi_{k}^{\dagger}\Psi_{k}$,
where $A_{k}$ and $B_{k}$ are fermionic annihilation operators and $\Psi_{k}$
are bosonic annihilation operators. The fermionic bath modes oscillate at
frequencies $\omega_{k}^{\Gamma}$ and couple to the spins with coupling
constants $\Gamma_{\alpha,k}$. Similarly the bosonic modes have frequencies
$\omega_{k}^{\zeta}$ and coupling constants $\zeta_{k}$.

The photonic bath represents a continuum of external electromagnetic modes,
for example on a superconducting transmission line in the microwave regime. We
assume that the bath modes are populated with photons in a prescribed spectral
distribution which can be realized by controllable external sources. (We shall
use the terms ‘photonic’ and ‘bosonic’ interchangeably throughout.) We assume
that the spectral distribution is broader than the natural linewidth of the
cavity mode and possesses random phases for different modes. Under these
conditions, this bosonic bath cannot impose any definite phase coherence on
the cavity but rather only drive it incoherently. Hence, the condensation
transition would be a result of spontaneous $U(1)$ symmetry breaking rather
than any externally imposed phase coherence.

It is important to note that under broad (Markovian) spectral distributions we
do not expect cavity photons to condense due to the bath’s excessive dephasing
powerSzymanska2001 . Given these conditions, we expect to find condensation
only for cases with a finite bandwidth. Hence we emphasize that the non-
Markovian, frequency dependent nature of the bath is crucial.

The fermionic baths, which are simpler to treat in field theory than the
bosonic bath, represent the microscopic environment of the spins. These baths
are assumed to be cold compared to the spin energy, hence in the absence of
any other processes would induce relaxation of the spins to the ground state.
The choice of symmetric spectral distributions ensures that on average there
are no pair-breaking events Szymanska2007  i.e. $\langle
b_{\alpha}^{\dagger}b_{\alpha}\rangle+\langle
a_{\alpha}^{\dagger}a_{\alpha}\rangle=1$. Although mathematically we define
two fermionic species comprising the fermionic bath, in reality this
corresponds to any kind of mechanism that causes spin relaxation and
dephasing.

##  III Methodology

Standard procedure of Keldysh field theory in its path-integral formulation
Szymanska2007 ; Sieberer2016  enables the Keldysh action to be derived.
Because the bath fields appear in the action at quadratic level, the
functional integral over them can be performed analytically so we simplify the
action by removing the bath degrees of freedom using Gaussian integration
Szymanska2007 ; Dunnett2016 . This yields an effective description in terms of
the photon field. In the non-equilibrium steady state the two-time Green’s
functions are time-translational invariant, $\tau=t-t^{\prime}$. This means we
can Fourier transform with respect to $\tau$ into the frequency representation
and make some standard assumptions about the bath properties.

Each bath couples equally to the cavity mode so that spontaneous emission is
momentum independent. The baths contain many modes, i.e. are much larger than
system. They therefore thermalize rapidly compared to any system interactions,
and their properties (e.g. distribution function) are unaffected by the system
behavior. Coupling strengths $\Gamma_{k}$ and $\zeta_{k}$ are smooth functions
of momentum $k$. The baths have a dense energy spectrum, so $\sum_{k}$
summations can be replaced by $\int\mathrm{d}\omega^{\Gamma}$ and
$\int\mathrm{d}\omega^{\zeta}$ integrals. We take the Markovian approximation
(frequency independence) for the bath density of states
($N^{\zeta}(\omega^{\zeta})=N^{\zeta}$ and
$N^{\Gamma}(\omega^{\Gamma})=N^{\Gamma}$) and spin-bath coupling
($\zeta(\omega^{\zeta})=\zeta$ and $\Gamma(\omega^{\Gamma})=\Gamma$).

In summary, this step introduces the cavity-photonic drive coupling strength
$\kappa=\pi\zeta^{2}N^{\zeta}$ and cavity-dephasing bath coupling strength
$\gamma=\pi\Gamma^{2}N^{\Gamma}$. As there are two different species of
fermions that couple to the upper and lower states of the spin-1/2 particle,
the latter is considered a relaxation rate that models the microscopic
environment. We also introduce three distribution functions: one for the
photonic drive $F_{\Psi}(\omega)$ and two for the dephasing baths
Szymanska2007  $F_{b/a}(\omega)$. Here $F_{\Psi}(\omega)=1+2n_{B}(\omega)$
where the photon occupation number $n_{B}(\omega)$ can have any form. We
discuss this later in the calculation.

The dephasing baths’ distribution functions are

$F_{b/a}(\omega)=1-2n^{b/a}_{F}(\omega).$ (4)

The occupation functions $n^{b/a}_{F}(\omega)$ are externally imposed and can
be chosen to have any form relevant to a particular physical situation. We
choose a thermal reservoir with $n^{b}_{F}$ and $n^{a}_{F}$ set at equal
internal temperature $T_{F}$ but different chemical potentials
$\mu^{b/a}_{B}$. The system must adhere to the single occupancy constraint
(eq. 2), so the chemical potential of baths $A$ and $B$ are related by
$\mu^{a}_{B}=-\mu^{b}_{B}=-\mu_{B}$, and the occupations should satisfy
Szymanska2006 ; Keeling2010  $n^{b}_{F}+n^{a}_{F}=1$. Hence we choose

$F_{b/a}(\omega)=\tanh\dfrac{\beta}{2}(\omega\mp\mu_{B})$ (5)

with inverse effective temperature $\beta=1/T_{F}$ and chemical potential
$\mu_{B}$.

In the regime where $\mu_{B}<0$ ($\mu_{B}>0$), the lower (upper) state is more
likely to be occupied in thermal equilibrium. The $\mu_{B}\to\mp\infty$ limit
corresponds to certainty that the lower/upper state is occupied. For $F_{b/a}$
to behave as a dephasing bath, we choose $\mu_{B}<0$.

For convenience, we arrange the fermionic fields in a Nambu vector
$\phi=[b,a]^{T}$, with $(1,2)$ indices to represent the fermionic fields under
Keldysh rotation Kamenev2011 . The resulting action is

$S=\iint\mathrm{d}t\mathrm{d}t^{\prime}\left(\sum_{\alpha,\alpha^{\prime}}\bar{\phi}_{\alpha}(t)\mathbf{G}_{\mathbf{0}}^{-1}(t,t^{\prime})\phi_{\alpha^{\prime}}(t^{\prime})+\bar{\psi}(t)\begin{bmatrix}0&i\partial_{t^{\prime}}-\omega_{0}-i\kappa\\\
i\partial_{t^{\prime}}-\omega_{0}+i\kappa&2i\kappa
F_{\Psi}(t-t^{\prime})\end{bmatrix}\psi(t^{\prime})\right)$ (6)

where $\mathbf{G_{0}}$ is the Tavis-Cummings Green’s function. We use
abbreviations $\lambda_{cl,q}=\frac{g_{\alpha}}{\sqrt{2}}\psi_{cl,q}$ so that
we may write

$\displaystyle\mathbf{G}_{\mathbf{0}}^{-1}(t,t^{\prime})$
$\displaystyle=\begin{bmatrix}(i\partial_{t^{\prime}}\sigma_{0}-\epsilon_{0}\sigma_{3}-\lambda_{cl}(t)\sigma_{+}-\bar{\lambda}_{cl}(t)\sigma_{-})\delta_{\alpha,\alpha^{\prime}}+i\gamma\sigma_{0}&-(\lambda_{q}(t)\sigma_{+}+\bar{\lambda}_{q}(t)\sigma_{-})\delta_{\alpha,\alpha^{\prime}}+2i\gamma(F_{b}\sigma_{\uparrow}+F_{a}\sigma_{\downarrow})\\\
-(\lambda_{q}(t)\sigma_{+}+\bar{\lambda}_{q}(t)\sigma_{-})\delta_{\alpha,\alpha^{\prime}}&(i\partial_{t^{\prime}}\sigma_{0}-\epsilon_{0}\sigma_{3}-\lambda_{cl}(t)\sigma_{+}-\bar{\lambda}_{cl}(t)\sigma_{-})\delta_{\alpha,\alpha^{\prime}}-i\gamma\sigma_{0}\end{bmatrix}.$
(7)

The calculation up to this point has followed the same steps as in previous
work Szymanska2007 . The rest of the calculation is original work.

To introduce the effects of incoherent driving, we separate the photon field
into its coherent and incoherent parts by substituting
$\psi=\langle\psi\rangle+\delta\psi$ into the action. Gaussian integration
removes the fluctuation degrees of freedom $\delta\psi$, and produces a new
action containing quartic $\sim\phi^{4}$ terms sm .

Hubbard-Stratonovich decoupling Altland2010  is required to simplify the
quartic terms, at the expense of introducing an auxiliary field
$Q^{(c,d)}(t,t^{\prime})_{\alpha,\alpha^{\prime}}\sim\langle\phi^{c}_{\alpha}(t)\bar{\phi}^{d}_{\alpha^{\prime}}(t^{\prime})\rangle$
which may be expressed in matrix form,

$\mathbf{Q}=\begin{bmatrix}Q^{(1,1)}&Q^{(1,2)}\\\
Q^{(2,1)}&Q^{(2,2)}\end{bmatrix}.$ (8)

The decoupling makes the action quadratic in $\phi$ and is analogous to the
BCS theory of superconductivity. The resulting action indicates that the
incoherent drive causes effective spin-spin interactions mediated by photon
exchange. This behavior appears in terms of the form sm

$\sim\iint\mathrm{d}t\mathrm{d}t^{\prime}\sum_{\alpha,\alpha^{\prime}}\bar{\phi}_{\alpha}(t)G_{\psi}^{R,A,K}(\pm\tau)Q^{(c,d)}_{\alpha,\alpha^{\prime}}(\tau)\phi_{\alpha^{\prime}}(t^{\prime})\\\
+G_{\psi}^{R,A,K}(\tau)Q^{(c,d)}_{\alpha,\alpha^{\prime}}(t,t^{\prime})Q^{(d,c)}_{\alpha^{\prime},\alpha}(t^{\prime},t).$

$G_{\psi}$ is the cavity photon Green’s function, and $R,A,K$ denote the
retarded, advanced and Keldysh components respectively Kamenev2011 ;
Sieberer2016 . By definition, the retarded and advanced components contain
information about the spectral properties whereas the Keldysh component
contains statistical properties of photons.

We can show explicitly that the linear terms in $Q$ influence the cavity’s
Tavis-Cummings dynamics sm . The Dyson equation becomes

$\mathbf{G}^{-1}=\mathbf{G}_{\mathbf{0}}^{-1}-\boldsymbol{\Sigma},$ (9)

where the Tavis-Cummings part sm2  is given by

$\mathbf{G}_{\mathbf{0}}^{-1}=\begin{bmatrix}(G_{0}^{-1})^{R}&(G_{0}^{-1})^{K}\\\
0&(G_{0}^{-1})^{A}\end{bmatrix}.$ (10)

The self-energy

$\boldsymbol{\Sigma}=\begin{bmatrix}\Sigma^{(1,1)}&\Sigma^{(1,2)}\\\
\Sigma^{(2,1)}&\Sigma^{(2,2)}\end{bmatrix}$ (11)

comes directly from the incoherent drive and contains linear terms in $Q$.
Note it is written in Keldysh $(1,2)$ and Nambu (particle-hole, or $(b,a)$)
space with non-trivial matrix elements $(i,j)=1,2$, where

$\displaystyle\Sigma^{(i,i)}_{bb/aa}(\tau)$
$\displaystyle=-i\dfrac{g^{2}}{4}\left[G_{\psi}^{K}(\pm\tau)Q^{(i,i)}_{aa/bb}(\tau)_{\alpha,\alpha^{\prime}}+G_{\psi}^{R/A}(\pm\tau)Q^{(i,j)}_{aa/bb}(\tau)_{\alpha,\alpha^{\prime}}+G_{\psi}^{A/R}(\pm\tau)Q^{(j,i)}_{aa/bb}(\tau)_{\alpha,\alpha^{\prime}}\right],i\neq
j$ $\displaystyle\Sigma^{(i,j)}_{bb/aa}(\tau)$
$\displaystyle=-i\dfrac{g^{2}}{4}\left[G_{\psi}^{K}(\pm\tau)Q^{(i,j)}_{aa/bb}(\tau)_{\alpha,\alpha^{\prime}}+G_{\psi}^{R/A}(\pm\tau)Q^{(i,i)}_{aa/bb}(\tau)_{\alpha,\alpha^{\prime}}+G_{\psi}^{A/R}(\pm\tau)Q^{(j,j)}_{aa/bb}(\tau)_{\alpha,\alpha^{\prime}}\right],i\neq
j.$

Following standard procedure, we take the narrow bandwidth limit
$g_{\alpha}\to g$ using the same motivations behind similar workMarchetti2006
, i.e. the individual spin-photon coupling is independent of the spin index
$\alpha$. Then we define the collective coupling as $g\to g\sqrt{N}$. For both
matrices $\mathbf{G}_{\mathbf{0}}^{-1}$ (eq. 10) and $\mathbf{\Sigma}$ (eq.
11), $b$ and $a$ define the particle-hole space as

$\boldsymbol{\mathcal{D}}=\begin{bmatrix}D_{bb}&D_{ba}\\\
D_{ab}&D_{aa}\end{bmatrix}.$

The Keldysh Green’s function $G^{K}_{\psi}$ contains the statistical
properties of the cavity photons which, after Fourier transforming into
frequency, is

$G^{K}_{\psi}(\omega)=\dfrac{-2i\kappa
F_{\Psi}(\omega)}{(\omega-\omega_{0})^{2}+\kappa^{2}}$ (12)

with system-photonic bath coupling strength $\kappa$. The bath distribution
function

$F_{\Psi}(\omega)=1+2n_{B}(\omega)$ (13)

is the energy resolved occupation of (quasi-)particle modes. Its presence in
$\mathbf{\Sigma}$ means that the photon environment, with $F_{\Psi}$ of modes
outside the cavity, now affects both the free photon evolution and the
development of spontaneous coherence. We will now show that the saddle-point
equation includes this distribution function, a significant departure from
previous works.

By adding a spin-spin interaction that we derive separately, we can construct
an augmented action,

$S_{\textrm{aug}}=\iint\mathrm{d}t\mathrm{d}t^{\prime}\sum_{\alpha,\alpha^{\prime}}\bar{\phi}_{\alpha}(t)\mathbf{G}^{-1}_{\alpha,\alpha^{\prime}}(t-t^{\prime})\phi_{\alpha^{\prime}}(t^{\prime})+S_{Q}+S_{\psi}$
(14)

where $\mathbf{G}$ is the Dyson equation defined in eq. 9 and $S_{Q}$ contains
terms quadratic in $Q$. Then integrating out the fermionic fields $\phi$
produces an effective action written in terms of the photon and auxiliary
fields,

$S_{\textrm{aug}}^{\textrm{eff}}=-i\sum_{\alpha}\mathrm{tr}\ln\mathbf{G}^{-1}_{\alpha,\alpha}+S_{Q}+S_{\psi}$
(15)

where tr traces over all parameters (Keldysh space, particle-hole space, time,
$\alpha$ site indices).

It is difficult to invert the $\mathbf{G}^{-1}$ matrix when the photon field
$\psi(t)$ has an arbitrary time dependence. Hence we take two further
simplifications: the single-frequency ansatz

$\psi(t)=\psi e^{-i\mu_{S}t}$ (16)

and saddle-point approximation. Since we are interested in non-equilibrium
steady states, we take the only time dependence of the photon field to be
oscillation at a single frequency $\mu_{S}$. This leads to explicit time
dependence within $\mathbf{G}^{-1}$. We want to search for self-consistent
solutions with a steady state, uniform, photon field of the form $\psi(t)$.
The $\mu_{S}$ parameter can also be thought of as the system’s chemical
potential. The mean-field theory of the non-equilibrium system describes a
self-consistent steady state, and inverting the $\mathbf{G}^{-1}$ matrix gives
a self-consistent Dyson equation

$\mathbf{G}(\omega)=\mathbf{G_{0}}(\omega)+\mathbf{G_{0}}(\omega)\cdot\boldsymbol{\Sigma}(\omega)\cdot\mathbf{G}(\omega).$
(17)

In summary, the effects of such a time dependence appear in two places: in the
time derivative terms, which lead to the energy shifts
$\omega_{0}\to\omega_{0}-\mu_{S}$ and $\epsilon_{0}\to\epsilon_{0}-\mu_{S}/2$,
and in a gauge transformation of the bath functions $F_{\Psi}(\omega)\to
F_{\Psi}(\omega+\mu_{S})$ and $F_{b/a}(\omega)\to
F_{b/a}(\omega\pm\mu_{S}/2)$.

It is customary to assume that the dominant contribution to the quantum
partition function arises from the configurations $\psi$ and $Q$ which
minimise the total action. This requires taking the functional derivative of
the action with respect to $\bar{\psi}_{q}$, then setting $\psi_{q}=0$:

$\dfrac{\delta
S_{\textrm{aug}}^{\textrm{eff}}}{\delta\bar{\psi}_{q}}\bigg{|}_{\psi_{q}=0}=0.$
(18)

We get the saddle point equation

$(\omega_{0}-\mu_{S}-i\kappa)\langle\psi\rangle=i\dfrac{g}{2}\sum_{\alpha}\mathrm{tr}\;G^{K}_{ba}(\tau)_{\alpha},$
(19)

where tr traces over all parameters. Note that there always exists a solution
to the saddle point equations where the quantum part is zero, which
corresponds to the purely classical limit. Hence when $\bar{\psi}_{q}=0$, then
$\delta S_{\textrm{aug}}^{\textrm{eff}}/\delta\bar{\psi}_{cl}=0$ trivially.

Similarly, varying the action with respect to $Q$,

$\dfrac{\delta S_{\textrm{aug}}^{\textrm{eff}}}{\delta Q^{(c,d)}}=0$ (20)

gives

$\mathbf{Q}_{\alpha,\alpha^{\prime}}(\tau)=-2\mathbf{G}_{\alpha,\alpha^{\prime}}(\tau).$
(21)

Hence we can redefine the auxiliary field with respect to the Dyson equation.
Substituting $Q$ into the Dyson equation (eq. 9) and following causality
arguments Kamenev2011  inherent to the formalism allows us to simplify the
self-energy. Hence we substitute $Q^{(1,1)}=Q^{R}$, $Q^{(1,2)}=Q^{K}$,
$Q^{(2,1)}=0$, and $Q^{(2,2)}=Q^{A}$. We rewrite the self-energy in matrix
form,

$\mathbf{\Sigma}=\begin{bmatrix}\Sigma^{R}&\Sigma^{K}\\\
0&\Sigma^{A}\end{bmatrix}$

where $\Sigma^{(1,1)}_{bb}=\Sigma^{R}_{bb}$,
$\Sigma^{(1,2)}_{bb}=\Sigma^{K}_{bb}$,
$\Sigma^{(2,1)}_{bb}=\Sigma^{(2,1)}_{aa}=0$,
$\Sigma^{(2,2)}_{bb}=\Sigma^{A}_{bb}$, $\Sigma^{(1,1)}_{aa}=\Sigma^{R}_{aa}$,
$\Sigma^{(1,2)}_{aa}=\Sigma^{K}_{aa}$, $\Sigma^{(2,2)}_{aa}=\Sigma^{A}_{aa}$.

The mean-field properties of a driven cavity are described by a complex
analogue of the Gross-Pitaevskii equation in the BEC regime or equivalently
the gap equation in the BCS regime Keeling2005 . The real part of the self-
consistent saddle-point equation (eq. 19) relates the coherent field to the
system’s nonlinear susceptibility, as in the case of equilibrium condensation,
while the imaginary part reflects how the gain and decay are balanced, as in a
laser. The equation contains two unknown parameters that characterize the
system: the coherent photon field $\langle\psi\rangle$ and its common
oscillation frequency $\mu_{S}$. The non-equilibrium Green’s function is
defined as $G^{K}_{ba}(t-t^{\prime})_{\alpha}=-i\langle
a_{cl}^{\dagger}(t)_{\alpha}b_{cl}(t^{\prime})_{\alpha}\rangle$, where
$a_{cl}=(a_{f}+a_{b})/\sqrt{2}$ and $f$ and $b$ are the forward and backward
branches of the Keldysh time contour Szymanska2007 ; Kamenev2011 .

As the name implies, the mean-field description (eq. 19) should only contain
details of the coherent photons, i.e. the saddle-point approximation filters
out the average field and discards any fluctuations. In contrast, our strategy
is to integrate over fluctuations $\delta\psi$ so that the effects of
incoherent photons appear in the action via corrections to the self-energy.
Then after taking the saddle-point approximation, the photonic drive
distribution function contributes at the mean field level.

The saddle-point equation itself is complex and hence contains two nonlinear
simultaneous equations (real and imaginary parts). It is necessary to use
numerical methods to solve for two unknowns $\mu_{S}$ and
$\langle\psi\rangle$. We use numerical solvers based on trust region
techniques sm ; Nocedal2006 , a simple yet robust method in optimization
problems. For a multi-valued function $Y$, the goal is to search for some
vector $x$ that satisfies the $Y(x)=0$ constraint or $Y(x)\approx 0$ if a
defined boundary condition is satisfied. It is necessary to locally
approximate the least-squares objective and constraint functions with some
simple functions that are easy to optimize. The trust region algorithm uses a
quadratic approximation of the objective function, then optimizes the
functions within some localized trust region around a point $x^{\prime}$ to
obtain a candidate step $y$. If $y$ is an improvement on the previous
iteration of $y$, the algorithms moves closer to the centre of the trust
region. If not, the approximations and trust region is adjusted.

##  IV Non-equilibrium condensation

By considering the external photon mode distribution function (eq. 13) with
specific frequency profiles, we now show that varying the drive parameters
does indeed permit the development of spontaneous coherence. For this work, we
set the photon occupation number $n_{B}$ to be a Lorentzian

$n_{B}(\omega)=h\dfrac{2\Omega}{[\omega-(\xi-\mu_{S})]^{2}+\Omega^{2}},$ (22)

centered at $\xi$ (resonance), with a width $\Omega$ (drive profile width) and
amplitude $h$. For simplicity, all parameters are scaled to $g=1$. As
discussed above, we assume a cold fermionic (dephasing) bath and hence we fix
its internal effective temperature to be $T_{F}=0.1g$. We also let
$\epsilon_{0}=0$ with zero detuning $\omega_{0}-2\epsilon_{0}=0$.

The following phase diagrams show the system properties such as the common
oscillation frequency of cavity photons $\mu_{S}$, photon coherence
$\langle\psi\rangle$, absolute value of polarization $|\langle
a^{\dagger}b\rangle|=\sqrt{(\omega_{0}-\mu_{S})^{2}+\kappa^{2}}|\psi|$
(fraction of condensed spins), and density of spin excitations
$\rho=\frac{1}{2}\sum_{\alpha}(1+b_{\alpha}^{\dagger}b_{\alpha}-a_{\alpha}^{\dagger}a_{\alpha})$.
The minimum $\rho=0$ corresponds to all spins in the ground state, whereas
$\rho=1$ corresponds to all spins in the excited state, i.e. maximum
inversion. In between, there is lasing when $\rho>0.5$, non-equilibrium
condensation when $\rho<0.5$, and a smooth crossover connecting the two
regimes.

Figure 2: System observables of driven system. Parameters
$[\mu_{B},\kappa,\gamma]=[-0.5,0.25,0.75]g$ and Lorentzian distribution
$n_{B}$ with $[\xi,\Omega]$. Solid lines indicate drive profile width
$\Omega=10\kappa$, dashed lines $\Omega=5\kappa$. Larger $\Omega$ requires
stronger driving (larger $h$) for system to condense. Driving off-resonance
$\xi\neq 0$ also requires larger $h$ to reach condensation transition. Markers
’x’ on $\langle\psi\rangle$ and polarization graphs indicate where the
condensation transition occurs.

Fig. 2, 3 and 4 show the system’s steady state observables upon increasing the
photonic drive strength $h$. The condensation transition occurs when
polarization and $\langle\psi\rangle$ become finite at some critical drive
strength. In particular, setting chemical potential $\mu_{B}/g<0$ ensures that
the detuning bath causes spontaneous decay processes instead of populating the
higher energy spin state. As expected, a larger driving amplitude $h$ is
required to reach the condensation transition when the drive is off resonance.
The transition occurs at the same critical drive strength when $\xi=\pm 1$,
indicating a symmetry in the system. However the different values for
$\mu_{S}/g$ and $\langle\psi\rangle$ suggest that off-resonance driving
affects the condensed state differently when $\xi=\pm 1$ (Fig. 2). Decreasing
the drive profile width $\Omega$ increases the number of photons per mode
around the $\xi$ frequency. This effectively allows the condensation
transition to occur at a smaller drive amplitude. We also see that larger
$\Omega$ increases the condensation threshold. Hence we do not expect the
system to condense in the Markovian limit (i.e. large $h$ and
$\Omega\to\infty$).

Figure 3: Driven system with Lorentzian distribution $n_{B}$ of external
photons, parameters
$[\mu_{B},\kappa,\gamma,\Omega]=[-0.5,0.25,0.75,10\kappa]g$. Pink lines
indicate the phase boundary where condensation transition occurs.

At some upper limit of the drive amplitude, the numerical method fails. This
possibly demonstrates that there is a region of inaccessible parameters. It
may be caused by a numerical issue, or the system could be entering some
unknown state that displays exotic, perhaps chaotic dynamics. Since the drive
amplitude controls the level of noise, it is also possible that the dephasing
overpowers the spin-photon dipole coupling and interrupts the condensation
processes.

Figure 4: Driven system with Lorentzian distribution $n_{B}$ of external
photons, parameters $[\mu_{B},\kappa,\gamma,\xi]=[-2,0.7,0.3,0]g$. Pink lines
indicate phase boundary of condensation transition.

This inaccessible region is more obvious in a two-dimensional phase diagram.
Fig. 3 uses the same parameters as in Fig. 2. Here we sweep the drive
amplitude $\xi$ and Lorentzian detuning. The phase diagram has an almost
symmetric phase boundary along the $\xi=0$ axis. It is clear that as the drive
detuning increases and becomes more off resonance, the drive amplitude must
increase in order for the system to condense. When the system is in a
condensed state, a large detuning and strong driving ($h=2$) regime results in
non-trivial differences in the system observables when $\xi=\pm 2$. Similarly,
Fig. 4 shows that as drive profile width $\Omega$ increases (which decreases
the number of photons per mode around $\xi=0$), the system remains in a
condensed state for larger ranges of the drive amplitude before entering the
region of inaccessible parameters. Note also that the spin density is always
$\rho<0.5$. In this regime the system is either in the condensed or normal
state. The density never reaches $\rho\geq 0.5$ which shows that the system
cannot behave as a laser. Instead as drive amplitude $h$ increases, the rate
of increase in photon coherence $\langle\psi\rangle$, polarization, and
density start levelling off before the system enters the region of
inaccessible parameters. In addition to taking a Lorentzian frequency
distribution, we examined a frequency independent or Markovian drive. We used
parameters that would be physically sensible for a system consisting of
superconducting resonator coupled to solid state spin-1/2 particles. We
explored different regimes involving combinations of $\gamma\leq\kappa$ or
$\gamma>\kappa$, and $\mu_{B}\leq 0$ or $\mu_{B}>0$, and gradually increased
the amount of Markovian noise (drive amplitude). We found no reliable
indications of condensation. This is an expected result, as the above phase
diagrams show how increasing the bandwidth $\Omega$ requires increasing the
driving strength for the system to reach the condensation transition. In
Szymańska’s work using Maxwell-Bloch equations Szymanska2001 , when all decay
baths are Markovian and finite, the system transitions from condensate-like
(no inversion) at zero decay to one requiring inversion. With frequency
independent driving and decay, there must be inversion for phase transition to
lasing. One major difference is that Szymańska’s work involved driving the
spins directly by coupling to the bath modes, so strong driving could easily
produce inversion. However if instead we drive the incoherent photons with
large energies up to infinite energies, it is highly unlikely to induce
population inversion or condensation.

##  V Summary and outlook

We used Keldysh path-integral field theory to study the open Tavis-Cummings
model. This allowed us to demonstrate the possibility of condensation and
spontaneous coherence under conditions of finite-bandwidth incoherent driving
with a Lorentzian frequency profile. It is crucial to use a frequency
dependent drive biased towards the lower frequencies. By exclusively driving
the low energy modes, we create ideal conditions for developing spontaneous
coherence. The resulting condensation is caused by the photonic drive which
induces an effective spin-spin coupling mediated by photons. This phenomenon
appears explicitly in the field theory and saddle-point equation. The
interactions arising from this incoherent drive may cause non-linear effects
that are not easily seen in other systems. Because the phase diagrams contain
areas where there are no solutions, the obvious next step would be to
investigate the possibility of different phases in those regimes, and to
explore the influence of fluctuations beyond the mean-field physics reported
heresm3 .

This theory can be used directly in state-of-the-art CQED experiments. Given
the microscopic system parameters and drive profile, quantitative predictions
can be made now for finding the non-equilibrium condensation transition, and
exploring novel non-equilibrium phases with finite-bandwidth drives. Since it
is challenging to find a parameter space where we can observe condensation,
further exploration would be extremely useful, as would understanding the
behavior seen in the phase diagrams. By applying our findings to analogue
quantum simulators, this could pave the way for fundamental research into
many-body physics and powerful new technological applications.

###### Acknowledgements.

R. A. acknowledges funding from University of Surrey’s University Research
Scholarship and Overseas Research Scholarship, and Universities New Zealand’s
Edward & Isabel Kidson Scholarship. E. G. acknowledges support from the
European Commission under project HiTIMe and EPSRC grants (Grant No.
EP/I026231/1 and EP/L02263X/1). We also acknowledge useful discussions with J.
M. J. Keeling. The data underlying this work are available without
restriction. Details of the data and how to request access are available from
the University of Surrey publications repository.

##  Appendix A Keldysh field theory calculation

This Appendix presents technical details related to the analytical procedure.
The bulk of the analysis involves Keldysh field theory in its functional
integral formulation. We use the Hubbard-Stratonovich transformation to
simplify the action, at the expense of introducing an auxiliary field.

The simplification ultimately allows the photon distribution to enter the
mean-field description. In the Keldysh formalism, this is evident in the
extension to the self-energy in the Dyson equation. Further simplification
involves assuming that the cavity photons oscillate at a single frequency
(taking the single-frequency ansatz). The final analytical step is mean field
analysis. This requires taking the functional derivative with respect to the
photon field to produce a complex Gross-Pitaevskii equation, and the auxiliary
field which makes it possible to write it in terms of the Dyson equation.

###  A.1 Building a Tavis-Cummings system

After constructing the Keldysh action using standard procedures Szymanska2007
, we arrange the fermionic fields into a Nambu vector $\phi=[b,a]^{T}$ and
$\phi=[\bar{b},\bar{a}]$. Then the quantum partition function is

$Z=\mathcal{N}\int\sum_{\alpha,k}\mathbf{D}[\bar{\phi}_{\alpha},\phi_{\alpha},\bar{\psi},\psi,\bar{A}_{k},A_{k},\bar{B}_{k},B_{k},\bar{\Psi}_{k},\Psi_{k}]\mathrm{e}^{iS}$
(23)

with normalisation constant $\mathcal{N}$, Keldysh action

$S=\int_{C}\mathrm{d}t\,\langle\Lambda(t)|i\partial_{t}-\hat{H}|\Lambda(t)\rangle,$
(24)

and coherent states specified by the fields in

$|\Lambda(t)\rangle=|\bar{\phi}_{\alpha}(t),\phi_{\alpha}(t),\bar{\psi}(t),\psi(t),\bar{A}_{k}(t),A_{k}(t),\bar{B}_{k}(t),B_{k}(t),\bar{\Psi}_{k}(t),\Psi_{k}(t)\rangle$
(25)

on Schwinger-Keldysh contour $C$. In the Keldysh formalism, the system evolves
from the distant past ($t=-\infty$) to the distant future ($t=\infty$) on the
forwards branch of a closed time contour and returns along the backwards
branch. For bosonic fields, the time evolution on the two branches is written
in terms of separate fields on each branch $\psi_{f,b}$ and rotated into the
classical/quantum basis via the Keldysh rotation Kamenev2011 ; Sieberer2016 ,

$\displaystyle\psi_{cl/q}=\dfrac{1}{\sqrt{2}}(\psi_{f}\pm\psi_{b}),$ (26)
$\displaystyle\bar{\psi}_{cl/q}=\dfrac{1}{\sqrt{2}}(\bar{\psi}_{f}\pm\bar{\psi}_{b}).$
(27)

Similarly the fermionic fields $\phi$ are rotated into the Larkin-Ovchinnikov
basis Kamenev2011

$\displaystyle\phi^{(1,2)}=\dfrac{1}{\sqrt{2}}(\phi_{f}\pm\phi_{b}),$ (28)
$\displaystyle\bar{\phi}^{(1,2)}=\dfrac{1}{\sqrt{2}}(\bar{\phi}_{f}\mp\bar{\phi}_{b}).$
(29)

Integrating out the baths generates a Keldysh action in terms of the photon
and fermion fields, $S=S_{\phi}+S_{\psi}$ with

$S_{\phi}=\iint\mathrm{d}t\mathrm{d}t^{\prime}\sum_{\alpha,\alpha^{\prime}}\bar{\phi}_{\alpha}(t)\mathbf{G}_{\mathbf{0}}^{-1}(t,t^{\prime})\phi_{\alpha^{\prime}}(t^{\prime}),$
(30)

and Tavis-Cummings Green’s function $\mathbf{G}_{\mathbf{0}}^{-1}$.

The photon action is

$S_{\psi}=\iint\mathrm{d}t\mathrm{d}t^{\prime}\,\bar{\psi}(t)\mathbf{G}_{\boldsymbol{\psi}}^{-1}(t,t^{\prime})\psi(t^{\prime})$
(31)

with photon Green’s function

$\displaystyle\mathbf{G}_{\boldsymbol{\psi}}^{-1}(t,t^{\prime})$
$\displaystyle=\begin{bmatrix}0&(G_{\psi}^{-1})^{A}\\\
(G_{\psi}^{-1})^{R}&(G_{\psi}^{-1})^{K}\end{bmatrix}$ (32)
$\displaystyle=\begin{bmatrix}0&i\partial_{t^{\prime}}-\omega_{0}-i\kappa\\\
i\partial_{t^{\prime}}-\omega_{0}+i\kappa&2i\kappa
F_{\Psi}(t-t^{\prime})\end{bmatrix}.$ (33)

As an aside: The Fourier transform Mahan1990  is defined for some arbitrary
complex function $x(\tau)$ in time $\tau=t-t^{\prime}$,

$\displaystyle
X(\omega)=\int_{-\infty}^{\infty}\mathrm{d}\tau\;\mathrm{e}^{i\omega\tau}x(\tau),$
(34) $\displaystyle
x(\tau)=\int_{-\infty}^{\infty}\dfrac{\mathrm{d}\omega}{2\pi}\mathrm{e}^{-i\omega\tau}X(\omega).$
(35)

The single frequency ansatz makes it possible to invert
$\mathbf{G}_{\mathbf{0}}^{-1}$. We let $\psi(t)=\psi\mathrm{e}^{-i\mu_{S}t}$,
hence the photon field $\psi$ time dependenceSzymanska2007  occurs as
oscillation at a single frequency $\mu_{S}$. Then using the Fourier transform
gives the Tavis-Cummings terms in frequency domain,

$\mathbf{G_{0}}=\begin{bmatrix}(G^{R}_{bb})_{0}&(G^{R}_{ba})_{0}&(G^{K}_{bb})_{0}&(G^{K}_{ba})_{0}\\\
(G^{R}_{ab})_{0}&(G^{R}_{aa})_{0}&(G^{K}_{ab})_{0}&(G^{K}_{aa})_{0}\\\
0&0&(G^{A}_{bb})_{0}&(G^{A}_{ba})_{0}\\\
0&0&(G^{A}_{ab})_{0}&(G^{A}_{aa})_{0}\end{bmatrix}$ (36)

where

$\displaystyle G^{R/A}_{0}(\omega)=\dfrac{(\omega\pm
i\gamma)\sigma_{0}+\epsilon\sigma_{3}+g_{\alpha}(\psi\sigma_{+}+\bar{\psi}\sigma_{-})}{\omega^{2}-E_{\alpha}^{2}\pm
2i\gamma\omega-\gamma^{2}}$ (37) $\displaystyle
G^{K}_{bb/aa}(\omega)_{0}=-2i\gamma\dfrac{((\omega\pm\epsilon)^{2}+\gamma^{2})F_{b/a}(\omega)+g_{\alpha}^{2}|\psi_{f}|^{2}F_{a/b}(\omega)}{((\omega-
E_{\alpha})^{2}+\gamma^{2})((\omega+E_{\alpha})^{2}+\gamma^{2})}$ (38)
$\displaystyle G^{K}_{ba}(\omega)_{0}=-G^{K}_{ab}(\omega)_{0}^{*}=-2i\gamma
g_{\alpha}\psi_{f}\dfrac{(\omega+\epsilon+i\gamma)F_{b}(\omega)+(\omega-\epsilon-i\gamma)F_{a}(\omega)}{((\omega-
E_{\alpha})^{2}+\gamma^{2})((\omega+E_{\alpha})^{2}+\gamma^{2})},$ (39)

where $\epsilon=\epsilon_{0}-\mu_{S}/2$ and
$E_{\alpha}=\sqrt{\epsilon^{2}+g_{\alpha}^{2}|\psi_{f}|^{2}}$. (When taking
the saddle-point approximation of the action $S=S_{\phi}+S_{\psi}$, this
requires setting $\psi_{q}=0$ so that
$\psi_{f}=(\psi_{cl}+\psi_{q})/\sqrt{2}=\psi_{cl}/\sqrt{2}$, and we recover
the $RAK$ Keldysh matrix structure.)

Similarly the photon Green’s functions are

$\displaystyle\mathbf{G}_{\boldsymbol{\psi}}=\begin{bmatrix}G_{\psi}^{K}&G_{\psi}^{R}\\\
G_{\psi}^{A}&0\end{bmatrix}$ (40) $\displaystyle
G^{R/A}_{\psi}(\omega)=\dfrac{1}{\omega-\omega_{0}\pm i\kappa},$ (41)
$\displaystyle G^{K}_{\psi}(\omega)=\dfrac{-2i\kappa
F_{\Psi}(\omega)}{(\omega-\omega_{0})^{2}+\kappa^{2}}.$ (42)

###  A.2 Addition of incoherent photon field

After integrating out the bath degrees of freedom, we separate the cavity
photon field into its coherent and incoherent parts
$\psi=\langle\psi\rangle+\delta\psi$, where the latter component is the result
of an externally imposed optical drive.

Again using Gaussian integration to remove the fluctuation fields produces
several terms. The analysis shows that the optical driving produces effective
spin-spin coupling mediated by exchange of photons,

$\displaystyle S_{\phi\phi}$
$\displaystyle=-\dfrac{g^{2}}{2}\iint\mathrm{d}t\mathrm{d}t^{\prime}\sum_{\alpha,\alpha^{\prime}}\biggr{(}\left(\bar{b}_{\alpha}^{(1)}a_{\alpha}^{(1)}+\bar{b}_{\alpha}^{(2)}a_{\alpha}^{(2)}\right)_{t}G_{\psi}^{K}(t,t^{\prime})\left(\bar{a}_{\alpha^{\prime}}^{(1)}b_{\alpha^{\prime}}^{(1)}+\bar{a}_{\alpha^{\prime}}^{(2)}b_{\alpha^{\prime}}^{(2)}\right)_{t^{\prime}}$
$\displaystyle\hskip
80.00012pt+\left(\bar{b}_{\alpha}^{(1)}a_{\alpha}^{(1)}+\bar{b}_{\alpha}^{(2)}a_{\alpha}^{(2)}\right)_{t}G_{\psi}^{R}(t,t^{\prime})\left(\bar{a}_{\alpha^{\prime}}^{(2)}b_{\alpha^{\prime}}^{(1)}+\bar{a}_{\alpha^{\prime}}^{(1)}b_{\alpha^{\prime}}^{(2)}\right)_{t^{\prime}}$
$\displaystyle\hskip
80.00012pt+\left(\bar{b}_{\alpha}^{(1)}a_{\alpha}^{(2)}+\bar{b}_{\alpha}^{(2)}a_{\alpha}^{(1)}\right)_{t}G_{\psi}^{A}(t,t^{\prime})\left(\bar{a}_{\alpha^{\prime}}^{(1)}b_{\alpha^{\prime}}^{(1)}+\bar{a}_{\alpha^{\prime}}^{(2)}b_{\alpha^{\prime}}^{(2)}\right)_{t^{\prime}}\biggr{)}.$
(43)

After reorganising and symmetrising the action, we use the Hubbard-
Stratonovich transformation Altland2010  to decouple the quartic
$\sim\phi^{4}$ terms. The process introduces an additional auxiliary field
$Q^{cd}(t,t^{\prime})_{\alpha,\alpha^{\prime}}\sim\langle\phi^{c}_{\alpha}(t)\bar{\phi}^{d}_{\alpha^{\prime}}(t^{\prime})\rangle$:

$\displaystyle
S_{\phi\phi}=\iint\mathrm{d}t\mathrm{d}t^{\prime}\sum_{\alpha,\alpha^{\prime}}\dfrac{g^{2}}{8}G_{\psi}^{K}(t,t^{\prime})Q^{(1,1)}_{aa}(t,t^{\prime})_{\alpha,\alpha^{\prime}}Q^{(1,1)}_{bb}(t^{\prime},t)_{\alpha^{\prime},\alpha}$
$\displaystyle\hskip
60.00009pt-i\dfrac{g^{2}}{4}\left(\bar{b}_{\alpha}^{(1)}(t)G_{\psi}^{K}(t,t^{\prime})Q^{(1,1)}_{aa}(t,t^{\prime})_{\alpha,\alpha^{\prime}}b_{\alpha^{\prime}}^{(1)}(t^{\prime})+\bar{a}_{\alpha}^{(1)}(t)G_{\psi}^{K}(t^{\prime},t)Q^{(1,1)}_{bb}(t,t^{\prime})_{\alpha,\alpha^{\prime}}a_{\alpha^{\prime}}^{(1)}(t^{\prime})\right)$
$\displaystyle\hskip
40.00006pt+\dfrac{g^{2}}{8}G_{\psi}^{K}(t,t^{\prime})Q^{(1,2)}_{aa}(t,t^{\prime})_{\alpha,\alpha^{\prime}}Q^{(2,1)}_{bb}(t^{\prime},t)_{\alpha^{\prime},\alpha}$
$\displaystyle\hskip
60.00009pt-i\dfrac{g^{2}}{4}\left(\bar{b}_{\alpha}^{(1)}(t)G_{\psi}^{K}(t,t^{\prime})Q^{(1,2)}_{aa}(t,t^{\prime})_{\alpha,\alpha^{\prime}}b_{\alpha^{\prime}}^{(2)}(t^{\prime})+\bar{a}_{\alpha}^{(2)}(t)G_{\psi}^{K}(t^{\prime},t)Q^{(2,1)}_{bb}(t,t^{\prime})_{\alpha,\alpha^{\prime}}a_{\alpha^{\prime}}^{(1)}(t^{\prime})\right)$
$\displaystyle\hskip
40.00006pt+\dfrac{g^{2}}{8}G_{\psi}^{K}(t,t^{\prime})Q^{(2,1)}_{aa}(t,t^{\prime})_{\alpha,\alpha^{\prime}}Q^{(1,2)}_{bb}(t^{\prime},t)_{\alpha^{\prime},\alpha}$
$\displaystyle\hskip
60.00009pt-i\dfrac{g^{2}}{4}\left(\bar{b}_{\alpha}^{(2)}(t)G_{\psi}^{K}(t,t^{\prime})Q^{(2,1)}_{aa}(t,t^{\prime})_{\alpha,\alpha^{\prime}}b_{\alpha^{\prime}}^{(1)}(t^{\prime})+\bar{a}_{\alpha}^{(1)}(t)G_{\psi}^{K}(t^{\prime},t)Q^{(1,2)}_{bb}(t,t^{\prime})_{\alpha,\alpha^{\prime}}a_{\alpha^{\prime}}^{(2)}(t^{\prime})\right)$
$\displaystyle\hskip
40.00006pt+\dfrac{g^{2}}{8}G_{\psi}^{K}(t,t^{\prime})Q^{(2,2)}_{aa}(t,t^{\prime})_{\alpha,\alpha^{\prime}}Q^{(2,2)}_{bb}(t^{\prime},t)_{\alpha^{\prime},\alpha}$
$\displaystyle\hskip
60.00009pt-i\dfrac{g^{2}}{4}\left(\bar{b}_{\alpha}^{(2)}(t)G_{\psi}^{K}(t,t^{\prime})Q^{(2,2)}_{aa}(t,t^{\prime})_{\alpha,\alpha^{\prime}}b_{\alpha^{\prime}}^{(2)}(t^{\prime})+\bar{a}_{\alpha}^{(2)}(t)G_{\psi}^{K}(t^{\prime},t)Q^{(2,2)}_{bb}(t,t^{\prime})_{\alpha,\alpha^{\prime}}a_{\alpha^{\prime}}^{(2)}(t^{\prime})\right)$
$\displaystyle\hskip
40.00006pt+\dfrac{g^{2}}{8}G_{\psi}^{R}(t,t^{\prime})Q^{(1,2)}_{aa}(t,t^{\prime})_{\alpha,\alpha^{\prime}}Q^{(1,1)}_{bb}(t^{\prime},t)_{\alpha^{\prime},\alpha}$
$\displaystyle\hskip
60.00009pt-i\dfrac{g^{2}}{4}\left(\bar{b}_{\alpha}^{(1)}(t)G_{\psi}^{R}(t,t^{\prime})Q^{(1,2)}_{aa}(t,t^{\prime})_{\alpha,\alpha^{\prime}}b_{\alpha^{\prime}}^{(1)}(t^{\prime})+\bar{a}_{\alpha}^{(2)}(t)G_{\psi}^{R}(t^{\prime},t)Q^{(1,1)}_{bb}(t,t^{\prime})_{\alpha,\alpha^{\prime}}a_{\alpha^{\prime}}^{(1)}(t^{\prime})\right)$
$\displaystyle\hskip
40.00006pt+\dfrac{g^{2}}{8}G_{\psi}^{R}(t,t^{\prime})Q^{(1,1)}_{aa}(t,t^{\prime})_{\alpha,\alpha^{\prime}}Q^{(2,1)}_{bb}(t^{\prime},t)_{\alpha^{\prime},\alpha}$
$\displaystyle\hskip
60.00009pt-i\dfrac{g^{2}}{4}\left(\bar{b}_{\alpha}^{(1)}(t)G_{\psi}^{R}(t,t^{\prime})Q^{(1,1)}_{aa}(t,t^{\prime})_{\alpha,\alpha^{\prime}}b_{\alpha^{\prime}}^{(2)}(t^{\prime})+\bar{a}_{\alpha}^{(1)}(t)G_{\psi}^{R}(t^{\prime},t)Q^{(2,1)}_{bb}(t,t^{\prime})_{\alpha,\alpha^{\prime}}a_{\alpha^{\prime}}^{(1)}(t^{\prime})\right)$
$\displaystyle\hskip
40.00006pt+\dfrac{g^{2}}{8}G_{\psi}^{R}(t,t^{\prime})Q^{(2,2)}_{aa}(t,t^{\prime})_{\alpha,\alpha^{\prime}}Q^{(1,2)}_{bb}(t^{\prime},t)_{\alpha^{\prime},\alpha}$
$\displaystyle\hskip
60.00009pt-i\dfrac{g^{2}}{4}\left(\bar{b}_{\alpha}^{(2)}(t)G_{\psi}^{R}(t,t^{\prime})Q^{(2,2)}_{aa}(t,t^{\prime})_{\alpha,\alpha^{\prime}}b_{\alpha^{\prime}}^{(1)}(t^{\prime})+\bar{a}_{\alpha}^{(2)}(t)G_{\psi}^{R}(t^{\prime},t)Q^{(1,2)}_{bb}(t,t^{\prime})_{\alpha,\alpha^{\prime}}a_{\alpha^{\prime}}^{(2)}(t^{\prime})\right)$
$\displaystyle\hskip
40.00006pt+\dfrac{g^{2}}{8}G_{\psi}^{R}(t,t^{\prime})Q^{(2,1)}_{aa}(t,t^{\prime})_{\alpha,\alpha^{\prime}}Q^{(2,2)}_{bb}(t^{\prime},t)_{\alpha^{\prime},\alpha}$
$\displaystyle\hskip
60.00009pt-i\dfrac{g^{2}}{4}\left(\bar{b}_{\alpha}^{(2)}(t)G_{\psi}^{R}(t,t^{\prime})Q^{(2,1)}_{aa}(t,t^{\prime})_{\alpha,\alpha^{\prime}}b_{\alpha^{\prime}}^{(2)}(t^{\prime})+\bar{a}_{\alpha}^{(1)}(t)G_{\psi}^{R}(t^{\prime},t)Q^{(2,2)}_{bb}(t,t^{\prime})_{\alpha,\alpha^{\prime}}a_{\alpha^{\prime}}^{(2)}(t^{\prime})\right)$
$\displaystyle\hskip
40.00006pt+\dfrac{g^{2}}{8}G_{\psi}^{A}(t,t^{\prime})Q^{(2,1)}_{aa}(t,t^{\prime})_{\alpha,\alpha^{\prime}}Q^{(1,1)}_{bb}(t^{\prime},t)_{\alpha^{\prime},\alpha}$
$\displaystyle\hskip
60.00009pt-i\dfrac{g^{2}}{4}\left(\bar{b}_{\alpha}^{(1)}(t)G_{\psi}^{A}(t,t^{\prime})Q^{(2,1)}_{aa}(t,t^{\prime})_{\alpha,\alpha^{\prime}}b_{\alpha^{\prime}}^{(1)}(t^{\prime})+\bar{a}_{\alpha}^{(1)}(t)G_{\psi}^{A}(t^{\prime},t)Q^{(1,1)}_{bb}(t,t^{\prime})_{\alpha,\alpha^{\prime}}a_{\alpha^{\prime}}^{(2)}(t^{\prime})\right)$
$\displaystyle\hskip
40.00006pt+\dfrac{g^{2}}{8}G_{\psi}^{A}(t,t^{\prime})Q^{(2,2)}_{aa}(t,t^{\prime})_{\alpha,\alpha^{\prime}}Q^{(2,1)}_{bb}(t^{\prime},t)_{\alpha^{\prime},\alpha}$
$\displaystyle\hskip
60.00009pt-i\dfrac{g^{2}}{4}\left(\bar{b}_{\alpha}^{(1)}(t)G_{\psi}^{A}(t,t^{\prime})Q^{(2,2)}_{aa}(t,t^{\prime})_{\alpha,\alpha^{\prime}}b_{\alpha^{\prime}}^{(2)}(t^{\prime})+\bar{a}_{\alpha}^{(2)}(t)G_{\psi}^{A}(t^{\prime},t)Q^{(2,1)}_{bb}(t,t^{\prime})_{\alpha,\alpha^{\prime}}a_{\alpha^{\prime}}^{(2)}(t^{\prime})\right)$
$\displaystyle\hskip
40.00006pt+\dfrac{g^{2}}{8}G_{\psi}^{A}(t,t^{\prime})Q^{(1,1)}_{aa}(t,t^{\prime})_{\alpha,\alpha^{\prime}}Q^{(1,2)}_{bb}(t^{\prime},t)_{\alpha^{\prime},\alpha}$
$\displaystyle\hskip
60.00009pt-i\dfrac{g^{2}}{4}\left(\bar{b}_{\alpha}^{(2)}(t)G_{\psi}^{A}(t,t^{\prime})Q^{(1,1)}_{aa}(t,t^{\prime})_{\alpha,\alpha^{\prime}}b_{\alpha^{\prime}}^{(1)}(t^{\prime})+\bar{a}_{\alpha}^{(1)}(t)G_{\psi}^{A}(t^{\prime},t)Q^{(1,2)}_{bb}(t,t^{\prime})_{\alpha,\alpha^{\prime}}a_{\alpha^{\prime}}^{(1)}(t^{\prime})\right)$
$\displaystyle\hskip
40.00006pt+\dfrac{g^{2}}{8}G_{\psi}^{A}(t,t^{\prime})Q^{(1,2)}_{aa}(t,t^{\prime})_{\alpha,\alpha^{\prime}}Q^{(2,2)}_{bb}(t^{\prime},t)_{\alpha^{\prime},\alpha}$
$\displaystyle\hskip
60.00009pt-i\dfrac{g^{2}}{4}\left(\bar{b}_{\alpha}^{(2)}(t)G_{\psi}^{A}(t,t^{\prime})Q^{(1,2)}_{aa}(t,t^{\prime})_{\alpha,\alpha^{\prime}}b_{\alpha^{\prime}}^{(2)}(t^{\prime})+\bar{a}_{\alpha}^{(2)}(t)G_{\psi}^{A}(t^{\prime},t)Q^{(2,2)}_{bb}(t,t^{\prime})_{\alpha,\alpha^{\prime}}a_{\alpha^{\prime}}^{(1)}(t^{\prime})\right),$
(44)

Taking the terms linear in $Q$, we can rewrite this part of $S_{\phi\phi}$ in
the style of (30), such that

$S_{\phi\phi}\bigg{|}_{\textrm{linear in
$Q$}}=\iint\mathrm{d}t\mathrm{d}t^{\prime}\sum_{\alpha,\alpha^{\prime}}\bar{\phi}_{\alpha}(t)\Sigma(t,t^{\prime})\phi_{\alpha^{\prime}}(t^{\prime}).$
(45)

For non-equilibrium systems in the steady state, the two-time Green’s
functions are time invariant ($\tau=t-t^{\prime}$) so self energy can be
written in matrix form,

$\boldsymbol{\Sigma}=\begin{bmatrix}\Sigma^{(1,1)}_{bb}&0&\Sigma^{(1,2)}_{bb}&0\\\
0&\Sigma^{(1,1)}_{aa}&0&\Sigma^{(1,2)}_{aa}\\\
\Sigma^{(2,1)}_{bb}&0&\Sigma^{(2,2)}_{bb}&0\\\
0&\Sigma^{(2,1)}_{aa}&0&\Sigma^{(2,2)}_{aa}\end{bmatrix}$ (46)

with elements

$\displaystyle\Sigma^{(1,1)}_{bb}(\tau)$
$\displaystyle=-i\dfrac{g^{2}}{4}\left(G_{\psi}^{K}(\tau)Q^{(1,1)}_{aa}(\tau)_{\alpha,\alpha^{\prime}}+G_{\psi}^{R}(\tau)Q^{(1,2)}_{aa}(\tau)_{\alpha,\alpha^{\prime}}+G_{\psi}^{A}(\tau)Q^{(2,1)}_{aa}(\tau)_{\alpha,\alpha^{\prime}}\right)$
(47) $\displaystyle\Sigma^{(1,2)}_{bb}(\tau)$
$\displaystyle=-i\dfrac{g^{2}}{4}\left(G_{\psi}^{K}(\tau)Q^{(1,2)}_{aa}(\tau)_{\alpha,\alpha^{\prime}}+G_{\psi}^{R}(\tau)Q^{(1,1)}_{aa}(\tau)_{\alpha,\alpha^{\prime}}+G_{\psi}^{A}(\tau)Q^{(2,2)}_{aa}(\tau)_{\alpha,\alpha^{\prime}}\right)$
(48) $\displaystyle\Sigma^{(2,1)}_{bb}(\tau)$
$\displaystyle=-i\dfrac{g^{2}}{4}\left(G_{\psi}^{K}(\tau)Q^{(2,1)}_{aa}(\tau)_{\alpha,\alpha^{\prime}}+G_{\psi}^{R}(\tau)Q^{(2,2)}_{aa}(\tau)_{\alpha,\alpha^{\prime}}+G_{\psi}^{A}(\tau)Q^{(1,1)}_{aa}(\tau)_{\alpha,\alpha^{\prime}}\right)$
(49) $\displaystyle\Sigma^{(2,2)}_{bb}(\tau)$
$\displaystyle=-i\dfrac{g^{2}}{4}\left(G_{\psi}^{K}(\tau)Q^{(2,2)}_{aa}(\tau)_{\alpha,\alpha^{\prime}}+G_{\psi}^{R}(\tau)Q^{(2,1)}_{aa}(\tau)_{\alpha,\alpha^{\prime}}+G_{\psi}^{A}(\tau)Q^{(1,2)}_{aa}(\tau)_{\alpha,\alpha^{\prime}}\right)$
(50) $\displaystyle\Sigma^{(1,1)}_{aa}(\tau)$
$\displaystyle=-i\dfrac{g^{2}}{4}\left(G_{\psi}^{K}(-\tau)Q^{(1,1)}_{bb}(\tau)_{\alpha,\alpha^{\prime}}+G_{\psi}^{R}(-\tau)Q^{(2,1)}_{bb}(\tau)_{\alpha,\alpha^{\prime}}+G_{\psi}^{A}(-\tau)Q^{(1,2)}_{bb}(\tau)_{\alpha,\alpha^{\prime}}\right)$
(51) $\displaystyle\Sigma^{(1,2)}_{aa}(\tau)$
$\displaystyle=-i\dfrac{g^{2}}{4}\left(G_{\psi}^{K}(-\tau)Q^{(1,2)}_{bb}(\tau)_{\alpha,\alpha^{\prime}}+G_{\psi}^{R}(-\tau)Q^{(2,2)}_{bb}(\tau)_{\alpha,\alpha^{\prime}}+G_{\psi}^{A}(-\tau)Q^{(1,1)}_{bb}(\tau)_{\alpha,\alpha^{\prime}}\right)$
(52) $\displaystyle\Sigma^{(2,1)}_{aa}(\tau)$
$\displaystyle=-i\dfrac{g^{2}}{4}\left(G_{\psi}^{K}(-\tau)Q^{(2,1)}_{bb}(\tau)_{\alpha,\alpha^{\prime}}+G_{\psi}^{R}(-\tau)Q^{(1,1)}_{bb}(\tau)_{\alpha,\alpha^{\prime}}+G_{\psi}^{A}(-\tau)Q^{(2,2)}_{bb}(\tau)_{\alpha,\alpha^{\prime}}\right)$
(53) $\displaystyle\Sigma^{(2,2)}_{aa}(\tau)$
$\displaystyle=-i\dfrac{g^{2}}{4}\left(G_{\psi}^{K}(-\tau)Q^{(2,2)}_{bb}(\tau)_{\alpha,\alpha^{\prime}}+G_{\psi}^{R}(-\tau)Q^{(1,2)}_{bb}(\tau)_{\alpha,\alpha^{\prime}}+G_{\psi}^{A}(-\tau)Q^{(2,1)}_{bb}(\tau)_{\alpha,\alpha^{\prime}}\right),$
(54)

This process shows that we make corrections to $\mathbf{G}^{-1}$ by producing
an additional self-energy. By adding an interaction that we derive separately,
we now construct an augmented action,

$S_{\textrm{aug}}=\iint\mathrm{d}t\mathrm{d}t^{\prime}\sum_{\alpha,\alpha^{\prime}}\bar{\phi}_{\alpha}(t)\mathbf{G}^{-1}_{\alpha,\alpha^{\prime}}(t-t^{\prime})\phi_{\alpha^{\prime}}(t^{\prime})+S_{Q}+S_{\psi}$
(55)

Then integrating out the fermionic fields $\phi$ produces an effective action,

$S_{\textrm{aug}}^{\textrm{eff}}=-i\sum_{\alpha}\mathrm{tr}\ln\mathbf{G}^{-1}_{\alpha,\alpha}+S_{Q}+S_{\psi}$
(56)

where tr traces over all parameters.

###  A.3 Saddle-point analysis

We calculate the mean field equation from the saddle points of the effective
action Szymanska2007 ; Kamenev2011 ; Altland2010  taken relative to the
classical and quantum fields,

$\displaystyle\dfrac{\delta
S_{\textrm{aug}}^{\textrm{eff}}}{\delta\bar{\psi}_{q}}\bigg{|}_{\psi_{q}=0}=0,$
(57) $\displaystyle\dfrac{\delta
S_{\textrm{aug}}^{\textrm{eff}}}{\delta\bar{\psi}_{cl}}\bigg{|}_{\psi_{q}=0}=0\quad\textrm{trivially},$
(58) $\displaystyle\dfrac{\delta S_{\textrm{aug}}^{\textrm{eff}}}{\delta
Q}=0.$ (59)

The $\delta S_{\textrm{aug}}^{\textrm{eff}}/\delta\bar{\psi}_{q}=0$ saddle-
point produces

$(\omega_{0}-\mu_{S}-i\kappa)\psi_{f}=i\dfrac{g}{2}\sum_{\alpha}\mathrm{tr}\;G^{K}_{ba}(\tau)_{\alpha},$
(60)

where $\psi_{f}=\psi_{cl}/\sqrt{2}$. The functional derivative $\delta
S_{\textrm{aug}}^{\textrm{eff}}/\delta Q=0$ produces

$\mathbf{Q}_{\alpha,\alpha^{\prime}}(\tau)=-2\mathbf{G}_{\alpha,\alpha^{\prime}}(\tau).$
(61)

Using causality arguments of retarded and advanced Green’s functions
Kamenev2011  leads to

$\displaystyle Q^{(1,1)}=Q^{R},$ (62) $\displaystyle Q^{(1,2)}=Q^{K},$ (63)
$\displaystyle Q^{(2,1)}=0,$ (64) $\displaystyle Q^{(2,2)}=Q^{A},$ (65)
$\displaystyle\Sigma^{(1,1)}_{bb}=\Sigma^{R}_{bb},$ (66)
$\displaystyle\Sigma^{(1,2)}_{bb}=\Sigma^{K}_{bb},$ (67)
$\displaystyle\Sigma^{(2,1)}_{bb}=\Sigma^{(2,1)}_{aa}=0,$ (68)
$\displaystyle\Sigma^{(2,2)}_{bb}=\Sigma^{A}_{bb},$ (69)
$\displaystyle\Sigma^{(1,1)}_{aa}=\Sigma^{R}_{aa},$ (70)
$\displaystyle\Sigma^{(1,2)}_{aa}=\Sigma^{K}_{aa},$ (71)
$\displaystyle\Sigma^{(2,2)}_{aa}=\Sigma^{A}_{aa}.$ (72)

We take the single frequency ansatz and get a self-consistent Dyson equation

$\mathbf{G}(\omega)=\mathbf{G_{0}}(\omega)+\mathbf{G_{0}}(\omega)\cdot\boldsymbol{\Sigma}(\omega)\cdot\mathbf{G}(\omega).$
(73)

The self-energy

$\displaystyle\boldsymbol{\Sigma}=\begin{bmatrix}\Sigma^{R}_{bb}&0&\Sigma^{K}_{bb}&0\\\
0&\Sigma^{R}_{aa}&0&\Sigma^{K}_{aa}\\\ 0&0&\Sigma^{A}_{bb}&0\\\
0&0&0&\Sigma^{A}_{aa}\end{bmatrix}$ (74)

in frequency is

$\displaystyle\Sigma^{R}_{bb}(\omega)$
$\displaystyle=i\dfrac{g^{2}}{4\pi}\bigg{(}G^{R}_{aa}(\omega)\circ
G_{\psi}^{K}(\omega)$ $\displaystyle\hskip 40.00006pt+G^{K}_{aa}(\omega)\circ
G_{\psi}^{R}(\omega)\bigg{)}$ (75) $\displaystyle\Sigma^{R}_{aa}(\omega)$
$\displaystyle=i\dfrac{g^{2}}{4\pi}\bigg{(}G^{R}_{bb}(\omega)\circ
G_{\psi}^{K}(-\omega)$ $\displaystyle\hskip 40.00006pt+G^{K}_{bb}(\omega)\circ
G_{\psi}^{A}(-\omega)\bigg{)}$ (76) $\displaystyle\Sigma^{A}_{bb}(\omega)$
$\displaystyle=i\dfrac{g^{2}}{4\pi}\bigg{(}G^{A}_{aa}(\omega)\circ
G_{\psi}^{K}(\omega)$ $\displaystyle\hskip 40.00006pt+G^{K}_{aa}(\omega)\circ
G_{\psi}^{A}(\omega)\bigg{)}$ (77) $\displaystyle\Sigma^{A}_{aa}(\omega)$
$\displaystyle=i\dfrac{g^{2}}{4\pi}\bigg{(}G^{A}_{bb}(\omega)\circ
G_{\psi}^{K}(-\omega)$ $\displaystyle\hskip 40.00006pt+G^{K}_{bb}(\omega)\circ
G_{\psi}^{R}(-\omega)\bigg{)}$ (78) $\displaystyle\Sigma^{K}_{bb}(\omega)$
$\displaystyle=i\dfrac{g^{2}}{4\pi}\bigg{(}G^{K}_{aa}(\omega)\circ
G_{\psi}^{K}(\omega)$ $\displaystyle\hskip 40.00006pt+G^{R}_{aa}(\omega)\circ
G_{\psi}^{R}(\omega)$ $\displaystyle\hskip 40.00006pt+G^{A}_{aa}(\omega)\circ
G_{\psi}^{A}(\omega)\bigg{)}$ (79) $\displaystyle\Sigma^{K}_{aa}(\omega)$
$\displaystyle=i\dfrac{g^{2}}{4\pi}\bigg{(}G^{K}_{bb}(\omega)\circ
G_{\psi}^{K}(-\omega)$ $\displaystyle\hskip 40.00006pt+G^{A}_{bb}(\omega)\circ
G_{\psi}^{R}(-\omega)$ $\displaystyle\hskip 40.00006pt+G^{R}_{bb}(\omega)\circ
G_{\psi}^{A}(-\omega)\bigg{)}$ (80)

where convolutions are defined as

$(f\circ g)(\omega)=\int\frac{\mathrm{d}\nu}{2\pi}f(\nu)g(\nu-\omega).$ (81)

Figure B1: Examples of self-energy $\Sigma$ as a function of frequency
$\omega$.

As a first approximation, we let $\mathbf{G}=\mathbf{G_{0}}$. Then we
numerically solve the saddle-point equation using Matlab’s fsolve function.
Through the trust region algorithm (either trust-region or trust-region-
dogleg), we solve for the following nonlinear simultaneous equations

$\displaystyle(\omega_{0}-\mu_{S})\psi_{f}-\Re\left(i\dfrac{g}{2}\int\dfrac{\textrm{d}\omega}{2\pi}G^{K}_{ba}(\omega)_{\alpha}\right)$
$\displaystyle=0$ (82)
$\displaystyle\kappa\psi_{f}-\Im\left(i\dfrac{g}{2}\int\dfrac{\textrm{d}\omega}{2\pi}G^{K}_{ba}(\omega)_{\alpha}\right)$
$\displaystyle=0$ (83)

which gives the numerical solution to $\mu_{S}$ and $\psi_{f}$. For
simplicity, we set $g=1$, $\omega=2\epsilon_{0}=0$, $T_{F}=0.1g$.

##  Appendix B Normal state solutions

We now consider the stability of uniform condensed solutions. The
consideration of stability is important because $\langle\psi\rangle=0$ is
always a solution of the gap equation, so it is necessary to determine which
of the normal and condensed solutions is stable. Secondly, because we
considered only spatially homogeneous fields with a single oscillation
frequency, it is possible that neither $\langle\psi\rangle=0$ nor the single-
frequency ansatz is stable, suggesting more complex behaviour. There is an
important difference in the interpretation of the saddle point equation
between the closed-time-path path-integral formalism used here and the
imaginary-time path integral in thermal equilibrium. In the imaginary time
formalism, extremising the action corresponds to finding configurations which
extremise the free energy. Thus, stable solutions correspond to a minimum of
free energy, and unstable to local maxima.

In contrast, for a classical saddle point (i.e., $\psi_{q}=0$), the action is
always $S=0$ and the saddle point condition corresponds to configurations for
which nearby paths add in phase. Thus, in order to study stability, one must
directly investigate fluctuations about our ansatz, and determine whether such
fluctuations grow or decay. In considering the question of stability, we
discuss stability of the normal state. This shows how the question of whether
fluctuations about the non-equilibrium steady-state grow or decay is directly
related to the instability expected in thermal equilibrium systems when the
chemical potential goes above a bosonic mode. We consider an expansion in
terms of $\psi=\langle\psi\rangle+\delta\psi$ to second order in $\delta\psi$.
To find the spectrum of fluctuations, we consider the effective action
governing fluctuations about either $\psi=\langle\psi\rangle$ or $\psi=0$.
Considering the effective action

$\displaystyle
S_{\textrm{aug}}^{\textrm{eff}}=-i\sum_{\alpha}\mathrm{tr}(\ln\mathbf{G}^{-1})+S_{Q}+S_{\psi},$
(84)
$\displaystyle\mathbf{G}^{-1}=\mathbf{G}_{\mathbf{0}}^{-1}-\boldsymbol{\Sigma}$
(85)

and expanding to second order in $\delta\psi$, one finds a contribution from
the effective photon action, and a contribution from expanding the trace over
fermions. Note $S_{Q}$ contains quadratic terms in $Q$, and
$\boldsymbol{\Sigma}$ contains linear terms in $Q$. Using bosonic basis:

$\mathbf{G}^{-1}=\mathbf{G}_{\mathbf{0}}^{-1}-\boldsymbol{\Sigma}=\begin{bmatrix}0&(G^{-1}_{0})^{A}\\\
(G^{-1}_{0})^{R}&(G^{-1}_{0})^{K}\end{bmatrix}-\begin{bmatrix}0&\Sigma^{A}\\\
\Sigma^{R}&\Sigma^{K}\end{bmatrix}.$ (86)

Expanding to second order in $\delta\psi$ gives contributions from $S_{\psi}$
and $-i\sum_{\alpha}\mathrm{tr}(\ln\mathbf{G}^{-1})$. The latter term gives

$\mathbf{G}_{\alpha,\alpha}^{-1}=(\mathbf{G}^{\mathbf{sp}}_{\alpha,\alpha})^{-1}+\boldsymbol{\delta}\mathbf{G}_{\alpha}^{-1},$
(87)

where $\mathbf{G^{\textrm{sp}}}$ is the saddle-point Green’s function, which
depends on the value of the saddle-point field $\langle\psi\rangle$. The
contribution from fluctuations is

$\displaystyle\boldsymbol{\delta}\mathbf{G}_{\alpha}^{-1}$
$\displaystyle=\boldsymbol{\delta}\mathbf{G}^{-1}_{\textbf{0};\alpha}-\boldsymbol{\delta\Sigma}_{\alpha}$
(88)
$\displaystyle=\dfrac{g}{\sqrt{2}}\left((\delta\bar{\psi}_{q}\cdot\sigma_{-}+\delta\psi_{q}\cdot\sigma_{+})\sigma_{1}^{K}+(\delta\bar{\psi}_{cl}\cdot\sigma_{-}+\delta\psi_{cl}\cdot\sigma_{+})\sigma_{0}^{K}\right)-\boldsymbol{\delta\Sigma}_{\alpha}$
(89)
$\displaystyle=\dfrac{g}{\sqrt{2}}\left((\delta\bar{\psi}_{q}\cdot\sigma_{-}+\delta\psi_{q}\cdot\sigma_{+})\sigma_{1}^{K}+(\delta\bar{\psi}_{cl}\cdot\sigma_{-}+\delta\psi_{cl}\cdot\sigma_{+})\sigma_{0}^{K}\right).$
(90)

The full Dyson equation becomes

$\mathbf{G}^{\mathbf{sp}}_{\alpha,\alpha}=\mathbf{G}^{\mathbf{sp}}_{\mathbf{0};\alpha,\alpha}+\mathbf{G}^{\mathbf{sp}}_{\mathbf{0};\alpha,\alpha}\cdot\boldsymbol{\Sigma}_{\alpha,\alpha}\cdot\mathbf{G}^{\mathbf{sp}}_{\alpha,\alpha}$
(91)

and $\boldsymbol{\Sigma}_{\alpha,\alpha}$ contains convolutions of
$\mathbf{G}^{\mathbf{sp}}_{\alpha,\alpha}$ and
$\mathbf{G}_{\boldsymbol{\psi}}$. Hence the action can be rewritten by
expanding

$\displaystyle-i\sum_{\alpha}\mathrm{tr}\ln(\mathbf{G}_{\alpha,\alpha}^{-1})$
$\displaystyle=-i\sum_{\alpha}\mathrm{tr}\ln((\mathbf{G}^{\mathbf{sp}}_{\alpha,\alpha})^{-1})-i\sum_{\alpha}\mathrm{tr}(\mathbf{G}^{\mathbf{sp}}_{\alpha,\alpha}\cdot\boldsymbol{\delta}\mathbf{G}_{\alpha}^{-1})+\dfrac{i}{2}\sum_{\alpha}\mathrm{tr}(\mathbf{G}^{\mathbf{sp}}_{\alpha,\alpha}\cdot\boldsymbol{\delta}\mathbf{G}_{\alpha}^{-1}\cdot\mathbf{G}^{\mathbf{sp}}_{\alpha,\alpha}\cdot\boldsymbol{\delta}\mathbf{G}_{\alpha}^{-1})$
(92)
$\displaystyle=-i\sum_{\alpha}\mathrm{tr}\ln((\mathbf{G}^{\mathbf{sp}}_{\alpha,\alpha})^{-1})-i\sum_{\alpha}\mathrm{tr}(\mathbf{G}^{\mathbf{sp}}_{\alpha,\alpha}\cdot\boldsymbol{\delta}\mathbf{G}_{0;\alpha}^{-1})+\dfrac{i}{2}\sum_{\alpha}\mathrm{tr}(\mathbf{G}^{\mathbf{sp}}_{\alpha,\alpha}\cdot\boldsymbol{\delta}\mathbf{G}_{0;\alpha}^{-1}\cdot\mathbf{G}^{\mathbf{sp}}_{\alpha,\alpha}\cdot\boldsymbol{\delta}\mathbf{G}_{0;\alpha}^{-1})$
(93)

and only retaining diagonal site index terms and neglecting bath-induced
interaction terms between different $\alpha$ sites.

Now, we consider the quadratic terms in $\delta\psi$ with

$\delta\Lambda=\begin{bmatrix}\delta\psi_{cl}(\omega)\\\
\delta\bar{\psi}_{cl}(-\omega)\\\ \delta\psi_{q}(\omega)\\\
\delta\bar{\psi}_{q}(-\omega)\end{bmatrix}$ (94)

the action for fluctuations becomes

$\delta
S_{f}=\int\dfrac{\mathrm{d}\omega}{2\pi}\delta\bar{\Lambda}(\omega)\begin{bmatrix}0&(D^{-1})^{A}\\\
(D^{-1})^{R}&(D^{-1})^{K}\end{bmatrix}\delta\Lambda(\omega)$ (95)

where

$(D^{-1})^{R,A,K}=\begin{bmatrix}K_{1}^{R,A,K}&K_{2}^{R,A,K}\\\
K_{3}^{R,A,K}&K_{4}^{R,A,K}\end{bmatrix}.$ (96)

In the normal state, Nambu ($2\times 2$) structure is redundant so the
distribution function is a diagonal constant matrix $F_{s}=1+n_{s}$, where
$n_{s}$ is the occupation of (bosonic) modes. There are also no anomalous
(off-diagonal in Nambu space) contributions, and so $K$ can be simplified
using Keldysh Green’s functions symmetry relations.

Taking the difference of luminescence and absorption gets the spectral weight

$2\pi
W(\omega,\mathbf{p}=0)=\dfrac{\Im(K^{R}_{1}(\omega))}{\Re(K^{R}_{1}(\omega,\mathbf{p}=0))^{2}+\Im(K^{R}_{1}(\omega))^{2}}$
(97)

where

$\displaystyle K^{R}_{1}(\omega)$
$\displaystyle=\dfrac{1}{2}(\omega-\tilde{\omega}_{0}+i\kappa)+i\dfrac{g^{2}}{4}\int\dfrac{\mathrm{d}\omega^{\prime}}{2\pi}\left(G^{R}_{bb}(\omega^{\prime})G^{K}_{aa}(\omega^{\prime}-\omega)+G^{K}_{bb}(\omega^{\prime})G^{A}_{aa}(\omega^{\prime}-\omega)\right).$
(98)

If the imaginary part of $K^{R}_{1}$ is a smooth function of $\omega$, then
the distribution function is

$F_{s}(\omega)=\dfrac{-iK^{K}_{1}(\omega)}{2\Im(K^{R}_{1}(\omega))}.$ (99)

If $\Im(K^{R}_{1}(\omega))=0$, then $F_{s}$ diverges. However the spectral
weight $W=0$ so that number of photons does not diverge. Hence in order to
find $\mu^{\textrm{eff}}$, it is necessary to solve for

$\Im(K^{R}_{1}(\mu^{\textrm{eff}}))=0.$ (100)

(a)

(b)

Figure B2: Normal state Green’s function $\Im(K^{R}_{1}(\omega))$ as a
function of frequency $\omega$ with and without incoherent drive. To determine
$\mu^{\textrm{eff}}$, it is necessary to find where the graph crosses the
horizontal axis.

As a simple example, let $\mathbf{G}=\mathbf{G_{0}}$. When
$\Im(K^{R}_{1}(\omega))$ is sketched as a function of frequency (Fig. B2a), we
find where the graph crosses the horizontal axis in order to determine
$\mu^{\textrm{eff}}$. As the dephasing bath becomes stronger, the
$\Im(K^{R}_{1}(\omega))$ graph is shifted upwards. This is observed as
$\mu_{B}$ decreases to negative values: physically it means that increasingly
more spins are in the ground state. At some critical $\mu_{B}$, the graph no
longer crosses the horizontal axis, hence $\mu^{\textrm{eff}}$ does not exist.
This problem also exists when we apply incoherent driving and define
$K^{R}_{1}(\omega)$ in terms of
$\mathbf{G}=\mathbf{G_{0}}+\mathbf{G_{0}}\cdot\boldsymbol{\Sigma}\cdot\mathbf{G}$
(Fig. B2b). The graph shapes can change dramatically when the dephasing
strength increases ($\mu_{B}$ decreases), compared to the simpler case without
bosonic driving. This makes it more difficult to determine whether
$\mu^{\textrm{eff}}$ exists, and if it is physical.

## References

  * (1) J. Kasprzak, M. Richard, S. Kundermann, A. Baas, P. Jeambrun, J. M. J. Keeling, F. M. Marchetti, M. H. Szymańska, R. Andrè, J. L. Staehli, et al., Nature 443, 409 (2006).
  * (2) M. H. Szymańska, J. Keeling, and P. B. Littlewood, Phys. Rev. B 75, 195331 (2007).
  * (3) D. Kremp, D. Semkat, and K. Henneberger, Phys. Rev. B 78, 125315 (2008).
  * (4) V. Zapf, M. Marcelo, and C. D. Batista, Rev. Mod. Phys. 86, 563 (2014).
  * (5) D. Marcos, A. Tomadin, S. Diehl, and P. Rabl, New J. Phys. 14, 055005 (2012).
  * (6) C. Wu, C. Guo, Y. Wang, G. Wang, X.-L. Feng, and J.-L. Chen, Phys. Rev. A 95, 013845 (2017).
  * (7) J. M. Fink, R. Bianchetti, M. Baur, M. Göppl, L. Steffen, S. Filipp, P. J. Leek, A. Blais, and A. Wallraff, Phys. Rev. Lett. 103, 083601 (2009).
  * (8) A. M. Zagoskin, Rev. Phys. 3, 1 (2018).
  * (9) J. Cirac and P. Zoller, Nat. Phys. 8, 264 (2012).
  * (10) Z.-L. Xiang, S. Ashhab, J. Q. You, and F. Nori, Rev. Mod. Phys. 85, 623 (2013).
  * (11) I. M. Georgescu, S. Ashhab, and F. Nori, Rev. Mod. Phys. 86, 153 (2014).
  * (12) E. Altman, K. R. Brown, G. Carleo, L. D. Carr, E. Demler, C. Chin, B. DeMarco, S. E. Economou, M. A. Eriksson, K.-M. C. Fu, et al., PRX Quantum 2, 017003 (2021).
  * (13) G. Kurizki, P. Bertet, Y. Kubo, K. Mølmer, D. Petrosyan, P. Rabl, and J. Schmiedmayer, Proc. Natl. Acad. Sci. 112, 3866 (2015).
  * (14) M. Fitzpatrick, N. M. Sundaresan, A. C. Y. Li, J. Koch, and A. A. Houck, Phys. Rev. X 7, 011016 (2017).
  * (15) A. Blais, A. L. Grimsmo, S. M. Girvin, and A. Wallraff, Rev. Mod. Phys. 93, 025005 (2021).
  * (16) C. Bonizzoni, A. Ghirri, and M. Affronte, Adv. Phys. X 3, 1435305 (2018).
  * (17) X. Mi, M. Benito, S. Putz1, D. M. Zajac, J. M. Taylor, G. Burkard, and J. R. Petta, Nature 555, 599 (2018).
  * (18) A. A. Clerk, K. W. Lehnert, P. Bertet, J. R. Petta, and Y. Nakamura, Nat. Phys. 16, 257 (2020).
  * (19) D. I. Schuster, A. P. Sears, E. Ginossar, L. DiCarlo, L. Frunzio, J. J. L. Morton, H. Wu, G. A. D. Briggs, B. B. Buckley, and D. D. Awschalom, Phys. Rev. Lett. 105, 140501 (2010).
  * (20) Y. Kubo, F. R. Ong, P. Bertet, D. Vion, V. Jacques, D. Zheng, A. Dréau, J.-F. Roch, A. Auffeves, and F. Jelezko, Phys. Rev. Lett. 105, 140502 (2010).
  * (21) R. Amsüss, S. Saito, and W. J. Munro, in Quantum information processing with diamond: Principles and applications, eds. S. Prawer and I. Aharonovich (Elsevier, 2014), pp. 264-289.
  * (22) Y. Liu, J. You, and Q. Hou, Sci. Rep. 6, 21775 (2016).
  * (23) G. Kirchmair, B. Vlastakis, Z. Leghtas, S. E. Nigg, H. Paik, E. Ginossar, M. Mirrahimi, L. Frunzio, S. M. Girvin, and R. J. Schoelkopf, Nature 495, 205 (2013).
  * (24) T. Weiss and O. Romero-Isart, Phys. Rev. Research 1, 033157 (2019).
  * (25) L. M. Sieberer, M. Buchhold, and S. Diehl, Rep. Prog. Phys. 79, 096001 (2016).
  * (26) H. Deng, H. Haug, and Y. Yamamoto, Rev. Mod. Phys. 82, 1489 (2010).
  * (27) T. Byrnes, N. Y. Kim, and Y. Yamamoto, Nat. Phys. 10, 803 (2014).
  * (28) B. Deveaud, Annu. Rev. Condens. Matter Phys. 6, 155 (2015).
  * (29) J. Keeling, P. R. Eastham, M. H. Szymańska, and P. B. Littlewood, Phys. Rev. B 72, 115320 (2005).
  * (30) J. Keeling, M. H. Szymańska, and P. B. Littlewood, in Optical Generation and Control of Quantum Coherence in Semiconductor Nanostructures, eds. G. Slavcheva and P. Roussignol (Springer, Berlin, 2010), pp. 293-329.
  * (31) M. H. Szymańska, J. Keeling, and P. B. Littlewood, in Quantum Gases: Finite Temperature and Non-Equilibrium Dynamics, eds. N. Proukakis, S. Gardiner, M. Davis, and M. H. Szymańska (Imperial College Press, London, 2013), pp. 447-459.
  * (32) A. Kamenev, Field Theory of Non-equilibrium Systems (Cambridge University Press, 2011).
  * (33) A. Altland and B. Simons, Condensed Matter Field Theory (Cambridge University Press, 2010), 2nd ed.
  * (34) M. H. Szymańska, PhD thesis, University of Cambridge, 2001\.
  * (35) P. Kirton and J. Keeling, New J. Phys. 20, 015009 (2018).
  * (36) P. Kirton, M. M. Roses, J. Keeling, and E. G. Dalla Torre, Adv. Quantum Technol. 2, 1800043 (2019).
  * (37) N. Nagaosa, Quantum Field Theory in Condensed Matter Physics (Springer, 1995).
  * (38) C. W. Gardiner and P. Zoller, Quantum Noise (Springer, 2000).
  * (39) K. Dunnett, and M. H. Szymańska, Phys. Rev. B 93, 195306 (2016).
  * (40) M. H. Szymańska, J. Keeling, and P. B. Littlewood, Phys. Rev. Lett. 96, 230602 (2006).
  * (41) See Appendix A for technical details of the calculation.
  * (42) The Tavis-Cummings behaviour $\mathbf{G_{0}}$ has been defined and studied in the context of microcavity polaritons Szymanska2007 , and we reproduce the $\mathbf{G_{0}}$ functions in Appendix A.
  * (43) F. M. Marchetti, J. Keeling, M. H. Szymańska, and P. B. Littlewood, Phys. Rev. Lett. 96, 066405 (2006).
  * (44) J. Nocedal, and S. J. Wright, Numerical Optimization (Springer, 2006).
  * (45) See Appendix B for a brief overview of the kinds of calculations involved.
  * (46) G. D. Mahan, Many-Particle Physics (Plenum Press, 1990).


# Test of the TMD gluon density in a proton with the longitudinal structure
function $F_{L}(x,Q^{2})$

A.V. Lipatov1,2, G.I. Lykasov2, M.A. Malyshev1,2

1Skobeltsyn Institute of Nuclear Physics, Lomonosov Moscow State University,
119991, Moscow, Russia
2Joint Institute for Nuclear Research, 141980, Dubna, Moscow region, Russia

Abstract

We investigate the dependence of the deep inelastic proton structure function
$F_{L}(x,Q^{2})$ on different forms of the transverse momentum dependent (TMD,
or unintegrated) gluon distribution. We present a comparison of theoretical
results with the latest experimental data taken by the H1 and ZEUS
Collaborations at HERA. We demonstrate that these data, despite of having
quite large uncertainties, can test different kinds of the gluon TMD in a
proton. Moreover, different phenomenological models at low $x$ could be tested
by future experiments on deep inelastic scattering.

Keywords: small-$x$ physics, high-energy factorization, deep inelastic
scattering, CCFM evolution, TMD gluon density

It is well known that the differential cross section of deep inelastic
electron-proton scattering (DIS) can be expressed in terms of proton structure
functions (SFs) $F_{2}(x,Q^{2})$ and $F_{L}(x,Q^{2})$. The longitudinal
structure function $F_{L}(x,Q^{2})$ is a very sensitive characteristic of
Quantum Chromodynamics (QCD) since it is directly related to the gluon content
of the proton. In fact, it is equal to zero in the naive parton model with
spin $1/2$ partons (quarks) and has got nonzero values in the framework of
perturbative QCD. Indeed, knowledge about the quark and, especially, the gluon
content of the proton is necessary for any theoretical study of high energy
processes performed within the QCD. Moreover, it is necessary for future
experiments on deep inelastic scattering planned at the Large Hadron electron
Collider (LHeC) and Future Circular hadron-electron Collider (FCC-he), where
facilities for using electron-proton center-of-mass energies $\sqrt{s}=1.3$
TeV and $3.5$ TeV are proposed. Usually, the number of collider data is
employed to constrain the gluon and quark density functions in a proton. In
this sense, an appropriate analysis of available HERA data on the structure
function $F_{L}(x,Q^{2})$ taken by the H1[1] and ZEUS[2] Collaborations at
$\sqrt{s}=318$ GeV is a rather important and interesting task.

However, calculated next-to-leading (NLO) pQCD corrections to the leading
order (LO) predictions for $F_{L}(x,Q^{2})$ are large and negative at small
$x$ (see, for example[3, 4]), that could even lead to negative
$F_{L}(x,Q^{2})$ values at small $x$ and low $Q^{2}$. The latter demonstrate
limitations of the applicability of perturbative QCD expansion and the
necessity of a resummation procedure (see[5, 6, 7, 8, 9, 10]). Such
resummation can be performed in the framework of so-called high energy QCD
factorization[11] (or $k_{T}$-factorization approach[12]) using the Catani-
Ciafaloni-Fiorani-Marchesini (CCFM) evolution equation[13] for transverse
momentum dependent (TMD, or unintegrated) gluon density in a proton. This
equation resumes large logarithmic terms proportional to
$\alpha_{s}^{n}\ln^{n}1/x$ and $\alpha_{s}^{n}\ln^{n}1/(1-x)$ and therefore
valid at both small and large $x$. The $k_{T}$-factorization QCD approach
supplemented with the CCFM gluon evolution dynamics is becoming an essential
tool which allows one to make theoretical predictions for processes studied at
modern and future colliders. More information can be found in review[14].

At present, a special interest is connected with investigation of the TMD
gluon density in a proton. In fact, while a great amount of knowledge about
the conventional parton density functions (PDFs), obeying Dokshitzer-Gribov-
Lipatov-Altarelli-Parisi (DGLAP) equations[15] has been accumulated over past
years, the TMD gluon and quark densities are still poorly known quantities.
Similar to DGLAP-based PDFs, the TMDs are usually parameterized in a general
form and then fitted to some experimental data. So, the latest precision HERA
measurements of the DIS structure functions $F_{2}(x,Q^{2})$ and
$F_{2}^{c}(x,Q^{2})$ were used[16] to determine the CCFM-evolved gluon
densities in a proton, where the empirical expression for TMD gluon density at
some starting scale $\mu_{0}^{2}$ (which is of the order of the hadron scale)
with factorized Gauss smearing in transverse momentum was applied111There are
other popular approaches to evaluate the TMDs, for example, the Kimber-Martin-
Ryskin (KMR) prescription[17, 18] or the Parton Branching (PB) approach[19,
20].. In contrast, in our previous study[21] the initial TMD gluon
distribution was chosen in the unfactorised form, which depends on $x$ and
$k_{T}$ simultaneously. This starting gluon density LLM’2022, like the GBW
gluon one [22, 23], results in the saturation of the dipole cross section at
low $Q^{2}$. Moreover, the form of the LLM’2022 input was chosen in order to
describe the LHC data on the charged hadron production at low transverse
momenta $p_{T}\sim 1$ GeV within the modified soft quark-gluon string model
(QGSM) [24, 25]. Then, the CCFM equation was applied to extend the proposed
TMD gluon density (LLM’2022) in the whole kinematical region. Several
parameters important at moderate and large $x$ have been fitted from the LHC
data on inclusive $b$-jet and Higgs boson production as well as the latest
HERA data on proton structure functions $F_{2}^{c}(x,Q^{2})$,
$F_{2}^{b}(x,Q^{2})$ and reduced cross sections $\sigma_{\rm
red}^{c}(x,Q^{2})$ and $\sigma_{\rm red}^{b}(x,Q^{2})$, thus moving forward to
the global fit of TMD gluon density from the collider data. However, the
measurements[1, 2] of the longitudinal structure function $F_{L}(x,Q^{2})$
performed by the H1 and ZEUS Collaborations at HERA were not taken into
account, since they have quite large experimental uncertainties as compared to
other data. The main goal of our present study is to test the LLM’2022 gluon
density with available data on $F_{L}(x,Q^{2})$ in a wide region of $x$ and
$Q^{2}$.

Let us briefly recall our main steps. We started from quite a general
expression for the initial TMD gluon distribution[21]:

$\displaystyle
f_{g}(x,{\mathbf{k}}_{T}^{2})=c_{g}(1-x)^{b_{g}}\sum_{n=1}^{3}\left(c_{n}R_{0}(x)|{\mathbf{k}}_{T}|\right)^{n}e^{-R_{0}(x)|{\mathbf{k}}_{T}|},$
(1)

where

$\displaystyle R_{0}(x)={1\over Q_{0}}\left({x\over
x_{0}}\right)^{\lambda/2},\quad
b_{g}=b_{g}(0)+{4C_{A}\over\beta_{0}}\ln{\alpha_{s}(Q_{0}^{2})\over\alpha_{s}({\mathbf{k}}_{T}^{2})},$
(2)

with $C_{A}=N_{C}$, $\beta_{0}=11-2N_{f}/3$, $x_{0}=4.1\cdot 10^{-5}$ and
$\lambda=0.22$. From the best description of recent LHC data on charged soft
hadron production in $pp$ collisions at low $p_{T}$ and the mid-rapidity
region (collected at $\sqrt{s}=0.9$, $2.36$, $7$ and $13$ TeV) within the
modified QGSM the parameters $c_{1}=5$, $c_{2}=3$, $c_{3}=2$ and $Q_{0}=1.233$
GeV, essential at small $x$, have been extracted. The values $b_{g}(0)=5.854$
and $c_{g}=0.173$ were determined from the simultaneous fit of HERA and LHC
data on some hard QCD processes performed in the framework of
$k_{T}$-factorization approach, as it was mentioned above. All the details of
calculations and fitting procedure can be found in ref.[21]. Here we only note
that a rather good value $\chi^{2}/d.o.f.=2.2$ is achieved. The obtained TMD
gluon density in a proton is shown in Figs. 1 and 2 as a function of the
proton longitudinal momentum fraction $x$ and gluon transverse momentum
${\mathbf{k}}_{T}^{2}$ for several values of hard scale $\mu^{2}$. The shaded
bands represent uncertainties of the fitting procedure, which become important
at large $x\geq 10^{-1}$. For comparison, we also show the CCFM-evolved TMD
gluon distribution JH’2013 set 2[16], which is often used in phenomenological
applications. In contrast with our approach, the $x$-dependence of JH’2013 set
2 input has quite a general form proportional to $x^{a}(1-x)^{b}$ with
parameters derived from the high-precision HERA data on the proton structure
functions $F_{2}(x,Q^{2})$ and $F_{2}^{c}(x,Q^{2})$ at $x<5\cdot 10^{-3}$ and
$Q^{2}>3.5$ GeV2. One can see that both TMD gluon densities have a remarkably
different $x$ and ${\mathbf{k}}_{T}^{2}$ behaviour, especially in the region
of small ${\mathbf{k}}_{T}^{2}$. As a consequence, the LLM’2022 gluon density
provides a better description of the considered hard QCD processes (see[21]).
Moreover, it saturates at a scale lower than the ones predicted by the popular
dipole Golec-Biernat-Wüsthoff (GBW) model[22, 23]. Currently, the LLM’2022
gluon density is implemented into the Monte-Carlo event generator pegasus[26]
and included in the tmdlib package[27], which is used by the cascade[28] and
KaTie[29] Monte-Carlo generators.

Figure 1: The TMD gluon densities in the proton
$f_{g}(x,{\mathbf{k}}_{T}^{2},\mu^{2})$ calculated as a function of the
longitudinal momentum fraction $x$ at different values of the transverse
momentum ${\mathbf{k}}_{T}^{2}$ and hard scale $\mu^{2}$. Shaded bands
represent the uncertainties of the $b_{g}(0)$ fitting procedure. Note that the
gluon densities calculated at $\mu^{2}=10^{4}$ GeV2 are multiplied by a factor
of $100$.

Figure 2: The TMD gluon densities in a proton
$f_{g}(x,{\mathbf{k}}_{T}^{2},\mu^{2})$ calculated as a function of the
transverse momentum ${\mathbf{k}}_{T}^{2}$ at different values of the
longitudinal momentum fraction $x$ and hard scale $\mu^{2}$. The notations are
the same as in Fig. 1.

In this paper we focus on investigation of the longitudinal SF
$F_{L}(x,Q^{2})$, which is related to the DIS differential cross section as
follows:

$\displaystyle\frac{d\sigma}{dxdy}=\frac{2\pi\alpha^{2}}{xQ^{4}}\left[\left(1-y+\frac{y^{2}}{2}\right)F_{2}(x,Q^{2})-\frac{y^{2}}{2}F_{L}(x,Q^{2})\right],$
(3)

where $Q^{2}$ is the virtuality of the probing photon, $x$ is the Bjorken
variable and $y$ is the inelasticity. Our calculations below are based on the
formulas[30] and here we strictly follow our previous consideration in all
aspects. So, according to the $k_{T}$-factorization prescription, the proton
structure function $F_{L}(x,Q^{2})$ can be calculated as a convolution

$\displaystyle F_{L}(x,Q^{2})=\int\limits_{x}^{1}\frac{dz}{z}\int
d{\mathbf{k}}_{T}^{2}\sum
e_{f}^{2}\hat{C}^{g}_{L}(x/z,Q^{2},m_{f}^{2},{\mathbf{k}}_{T}^{2})f_{g}(z,{\mathbf{k}}_{T}^{2},\mu^{2}),$
(4)

where the summation is performed over flavors $f$ of quarks with electric
charges $e_{f}$ and masses $m_{f}$. The hard coefficient function
$C^{g}_{L}(x/z,Q^{2},m_{f}^{2},{\mathbf{k}}_{T}^{2})$ corresponds to the
quark-box diagram for off-shell (dependent on the incoming gluon virtuality)
photon-gluon fusion subprocess and was calculated earlier[30]. Numerically, we
set the charm and beauty masses to $m_{c}=1.67$ GeV and $m_{b}=4.75$ GeV and
use the massless limit to evaluate the corresponding contributions from the
light quarks. Also, we apply the $2$-loop formula for the strong coupling
constant $\alpha_{s}$ with $N_{f}=4$ quark flavours at $\Lambda_{\rm QCD}=200$
MeV.

The results of our calculations are shown in Figs. 3 and 4. Here we plot the
evaluated $F_{L}(x,Q^{2})$ as a function of $x$ for different values of
$Q^{2}$ in comparison with available HERA data, taken by the H1[1] and ZEUS[2]
Collaborations. In addition, we also show the results obtained with the
JH’2013 set 2 gluon density. The shaded bands correspond to theoretical
uncertainties of our calculations connected with the choice of hard scales.
Note that in the case of the JH’2013 set 2 gluon density we use auxiliary
’$+$’ and ’$-$’ distributions corresponding to the variated renormalization
scales (see[16]). So, we find that the LLM’2022 gluon density describes the
HERA data quite well within the estimated theoretical and experimental
uncertainties. It tends to slightly overestimate the H1 measurements[1] at
large $Q^{2}$, but the predictions are still compatible with the data at $\sim
2\sigma$. We obtained, also, that the JH’2013 set 2 gluon density provides a
bit worse description of the HERA data, especially at low $Q^{2}$. The better
agreement achieved with the LLM’2022 distribution is an immediate consequence
of using the physically motivated expression (1) for the starting TMD gluon
density. In fact, at low $Q^{2}$ the relatively small gluon
${\mathbf{k}}_{T}^{2}$ are probed, where the difference between the considered
TMD gluon distributions becomes large, as it is shown in Fig. 2. This big
difference can be due to the saturation of the dipole cross section at low
$Q^{2}$ [22, 23], which is taken into account by the LLM’2022 gluon. So, our
calculations demonstrate that available experimental data for the proton
structure function $F_{L}(x,Q^{2})$ are sensitive to the TMD gluon densities,
and the relatively low $Q^{2}$ region could help to clearly distinguish the
latter.

To conclude, we have tested the recently proposed TMD gluon density in a
proton (LLM’2022) with the longitudinal structure function $F_{L}(x,Q^{2})$ in
a wide region of $x$ and $Q^{2}$. We have compared corresponding predictions
with the experimental data taken by the H1 and ZEUS Collaborations at HERA as
well as with results obtained with another CCFM-based gluon distribution,
JH’2013 set 2. We have found that the predictions with the LLM’2022 gluon
density do not contradict available data in a wide $Q^{2}$ range. Moreover,
they show a better description of the data than the JH’2013 set 2 distribution
due to more physically motivated expression for the initial TMD gluon density.
Thus, the measurements of $F_{L}(x,Q^{2})$ could verify different TMDs,
especially at not very large $Q^{2}$, where the further effects of QCD
evolution do not play an essential role. It could be important for experiments
at future colliders, such as LHeC or FCC-he, and for improvement of different
small-$x$ phenomenological models.

Figure 3:  The longitudinal SF $F_{2}(x,Q^{2})$ calculated at different
$Q^{2}$. The green band correspond to the result obtained with the LLM’2022
gluon density with scale uncertainties, while the red one represents results
calculated with the JH’2013 set 2 one. The experimental data are from H1[1].
Figure 4:  The longitudinal SF $F_{2}(x,Q^{2})$ calculated at different
$Q^{2}$. The notations are the same as in Fig. 3. The experimental data are
from ZEUS[2].

Acknowledgements. We thank S.P. Baranov, A.V. Kotikov, H. Jung, S. Taheri
Monfared for their important comments and remarks. We also thank G. Pontecorvo
for careful reading of our manuscript. Studies related with derivation of the
LLM’2022 gluon density were supported by the Russian Science Foundation under
grant 22-22-00387. All other performed calculations were supported by the
Russian Science Foundation, grant 22-22-00119.

## References

  * [1] H1 Collaboration, Eur. Phys. J. C 71, 1579 (2011).
  * [2] ZEUS Collaboration, Phys. Lett. B 682, 8 (2009).
  * [3] C.D. White, R.S. Thorne, Phys. Rev. D 74, 014002 (2006).
  * [4] C.D. White, R.S. Thorne, Phys. Rev. D 75, 034005 (2007).
  * [5] M. Ciafaloni, D. Colferai, G.P. Salam, Phys. Rev. D 60, 114036 (1999).
  * [6] M. Ciafaloni, D. Colferai, G.P. Salam, JHEP 07, 054 (2000).
  * [7] R.S. Thorne, Phys. Lett. B 474, 372 (2000).
  * [8] R.S. Thorne, Phys. Rev. D 60, 054031 (1999).
  * [9] R.S. Thorne, Phys. Rev. D 64, 074005 (2001).
  * [10] G. Altarelli, R.D. Ball, S. Forte, Nucl. Phys. B 621, 359 (2002).
  * [11] S. Catani, M. Ciafaloni, F. Hautmann, Nucl. Phys. B 366, 135 (1991);
J.C. Collins, R.K. Ellis, Nucl. Phys. B 360, 3 (1991).

  * [12] L.V. Gribov, E.M. Levin, M.G. Ryskin, Phys. Rep. 100, 1 (1983);
E.M. Levin, M.G. Ryskin, Yu.M. Shabelsky, A.G. Shuvaev, Sov. J. Nucl. Phys.
53, 657 (1991).

  * [13] M. Ciafaloni, Nucl. Phys. B 296, 49 (1988);
S. Catani, F. Fiorani, G. Marchesini, Phys. Lett. B 234, 339 (1990);
S. Catani, F. Fiorani, G. Marchesini, Nucl. Phys. B 336, 18 (1990);
G. Marchesini, Nucl. Phys. B 445, 49 (1995).

  * [14] R. Angeles-Martinez, A. Bacchetta, I.I. Balitsky, D. Boer, M. Boglione, R. Boussarie, F.A. Ceccopieri, I.O. Cherednikov, P. Connor, M.G. Echevarria, G. Ferrera, J. Grados Luyando, F. Hautmann, H. Jung, T. Kasemets, K. Kutak, J.P. Lansberg, A. Lelek, G.I. Lykasov, J.D. Madrigal Martinez, P.J. Mulders, E.R. Nocera, E. Petreska, C. Pisano, R. Placakyte, V. Radescu, M. Radici, G. Schnell, I. Scimemi, A. Signori, L. Szymanowski, S. Taheri Monfared, F.F. Van der Veken, H.J. van Haevermaet, P. Van Mechelen, A.A. Vladimirov, S. Wallon, Acta Phys. Polon. B 46, 2501 (2015).
  * [15] V.N. Gribov, L.N. Lipatov, Sov. J. Nucl. Phys. 15, 438 (1972);
L.N. Lipatov, Sov. J. Nucl. Phys. 20, 94 (1975);
G. Altarelli, G. Parisi, Nucl. Phys. B 126, 298 (1977);
Yu.L. Dokshitzer, Sov. Phys. JETP 46, 641 (1977).

  * [16] F. Hautmann, H. Jung, Nucl. Phys. B 883, 1 (2014).
  * [17] M.A. Kimber, A.D. Martin, M.G. Ryskin, Phys. Rev. D 63, 114027 (2001);
A.D. Martin, M.G. Ryskin, G. Watt, Eur. Phys. J. C 31, 73 (2003).

  * [18] A.D. Martin, M.G. Ryskin, G. Watt, Eur. Phys. J. C 66, 163 (2010).
  * [19] F. Hautmann, H. Jung, A. Lelek, V. Radescu, R. Zlebcik, Phys. Lett. B 772, 446 (2017).
  * [20] F. Hautmann, H. Jung, A. Lelek, V. Radescu, R. Zlebcik, JHEP 1801, 070 (2018).
  * [21] A.V. Lipatov, G.I. Lykasov, M.A. Malyshev, Phys. Rev. D 107, 014022 (2023).
  * [22] K. Golec-Biernat, M. Wüsthoff, Phys. Rev. D 59, 014017 (1998).
  * [23] K. Golec-Biernat, M. Wüsthoff, Phys. Rev. D 60, 114023 (1999).
  * [24] V.A. Bednyakov, G.I. Lykasov, V.V. Lyubushkin, Europhys. Lett. 92, 31001 (2010).
  * [25] V.A. Bednyakov, A.A. Grinyuk, G.I. Lykasov, M. Poghosyan, Int. J. Mod. Phys. A 27, 1250042 (2012).
  * [26] A.V. Lipatov, M.A. Malyshev and S.P. Baranov, Eur. Phys. J. C 80, 330 (2020).
  * [27] N.A. Abdulov, A. Bacchetta, S.P. Baranov, A. Bermudez Martinez, V. Bertone, C. Bissolotti, V. Candelise, L.I. Estevez Banos, M. Bury, P.L.S. Connor, L. Favart, F. Guzman, F. Hautmann, M. Hentschinski, H. Jung, L. Keersmaekers, A.V. Kotikov, A. Kusina, K. Kutak, A. Lelek, J. Lidrych, A.V. Lipatov, G.I. Lykasov, M.A. Malyshev, M. Mendizabal, S. Prestel, S. Sadeghi Barzani, S. Sapeta, M. Schmitz, A. Signori, G. Sorrentino, S. Taheri Monfared, A. van Hameren, A.M. van Kampen, M. Vanden Bemden, A. Vladimirov, Q. Wang, H. Yang, Eur. Phys. J. C 81, 752 (2021).
  * [28] H.Jung, S.P. Baranov, A. Bermudez Martinez, L.I. Estevez Banos, F. Guzman, F. Hautmann, A. Lelek, J. Lidrych, A.V. Lipatov, M.A. Malyshev, M. Mendizabal, S. Taheri Monfared, A.M. van Kampen, Q. Wang, H. Yang, Eur. Phys. J. C 81, 425 (2021).
  * [29] A. van Hameren, Comput. Phys. Commun. 224, 371 (2018).
  * [30] A.V. Kotikov, A.V. Lipatov, G. Parente, N.P. Zotov, Eur. Phys. J. C 26, 51 (2002).