#  Cavity QED materials: Comparison and validation of two linear response
theories at arbitrary light-matter coupling strengths

Juan Román-Roche  Instituto de Nanociencia y Materiales de Aragón (INMA),
CSIC-Universidad de Zaragoza, Zaragoza 50009, Spain  Departamento de Física de
la Materia Condensada, Universidad de Zaragoza, Zaragoza 50009, Spain  Álvaro
Gómez-León  Institute of Fundamental Physics IFF-CSIC, Calle Serrano 113b,
28006 Madrid, Spain  Fernando Luis  Instituto de Nanociencia y Materiales de
Aragón (INMA), CSIC-Universidad de Zaragoza, Zaragoza 50009, Spain
Departamento de Física de la Materia Condensada, Universidad de Zaragoza,
Zaragoza 50009, Spain  David Zueco  Instituto de Nanociencia y Materiales de
Aragón (INMA), CSIC-Universidad de Zaragoza, Zaragoza 50009, Spain
Departamento de Física de la Materia Condensada, Universidad de Zaragoza,
Zaragoza 50009, Spain

###### Abstract

We develop a linear response theory for materials collectively coupled to a
cavity that is valid in all regimes of light-matter coupling, including
symmetry-broken phases. We present and compare two different approaches.
First, using a coherent path integral formulation for the partition function
to obtain thermal Green functions. This approach relies on a saddle point
expansion for the action, that can be truncated in the thermodynamic limit.
Second, by formulating the equations of motion for the retarded Green
functions and solving them. We use a mean-field decoupling of high-order Green
functions in order to obtain a closed, solvable system of equations. Both
approaches yield identical results in the calculation of response functions
for the cavity and material. These are obtained in terms of the bare cavity
and material responses. In combination, the two techniques clarify the
validity of a mean-field decoupling in correlated light-matter systems and
provide complementary means to compute finite-size corrections to the
thermodynamic limit. The theory is formulated for a general model that
encompasses most of the systems typically considered in the field of cavity
QED materials, within a long-wavelength approximation. Finally, we provide a
detailed application of the theory to the Quantum Hall effect and to a
collection of magnetic models. We validate our predictions against analytical
and finite-size exact-diagonalization results.

##  I Introduction

Linear Response Theory (LRT) is an elegant tool for understanding how systems
at equilibrium respond when they are perturbed. The linearity refers to the
condition that the perturbation is sufficiently weak, allowing the system’s
dynamics to be treated to first order in the perturbation. This framework
yields various relations that determine response properties, regardless of the
specific situation or model considered. Fundamentally, it relates equilibrium
fluctuations and dissipation through the fluctuation-dissipation theorem
[1][Chap. 7].

In this paper, we investigate a particular scenario where a system containing
$N$ local degrees of freedom is coupled to a single-mode cavity. The coupling
is assumed to be collective, and our theory yields exact results in the
thermodynamic limit for the matter, $N\to\infty$. This is the typical limit
considered in cavity QED materials, on which we will elaborate below [2, 3,
4]. Beyond the assumption of collective coupling and the thermodynamic limit,
our approach is general. We obtain formulas for the response of both the
material and the cavity in terms of the bare photonic and material response
functions, for any light-matter coupling strength. At equilibrium, it has been
shown that the $N\to\infty$ limit simplifies the problem, allowing for exact
results, particularly famous is the existence (or lack thereof) of the
superradiant phase transition. Various techniques have been employed,
including thermodynamic inequalities [5], scaling arguments [6], the stiffness
theorem [7], and path integral methods [8]. In this work, we are interested in
perturbing the equilibrium to build an equally rigorous result for LRT. To
achieve this, we use two independent but complementary approaches: one based
on a path integral formulation of the partition function and the other on the
equations of motion for the correlators. The main idea underpinning our
derivation is the fact that collective interactions can be exactly decoupled
by an auxiliary mean-field variable. This is a common resource in large-$N$
theories in quantum field theory which is used in combination with a saddle-
point approximation for the action [9]. The complementary formulation in terms
of the equations of motion for the correlators renounces the provable
exactness afforded by the path integral in exchange for a simpler derivation
where the mean-field decoupling is explicit.

Cavity QED materials are actively discussed in the literature, particularly
the possibility of modifying the properties of matter using quantum light [2,
3]. This is pertinent in the regime of strong light-matter correlations, which
is why a non-perturbative approach is relevant. The problem is far from being
a theoretical curiosity, as initial experiments have reported modifications in
conductivity [10, 11], magnetism [12], and metal-to-insulator transition [13].
Remarkably, one way to observe these modifications is through the hybrid
light-matter response, which is where our work fits, providing a general
theory for such situations. In this context, our work presents the LRT for
cavity materials using both a path-integral and an equations-of-motion
approach. We discuss their relationship and the application of our theory. We
demonstrate that our theory recovers previous results on the cavity-modified
quantum Hall effect and explore the response of magnetic materials in
cavities, discussing dipolar Ising systems with longitudinal and transverse
coupling to the cavity, Heisenberg ferromagnetism, and the Dicke-Ising model.
The latter involves a one-dimensional Ising chain coupled to a cavity,
allowing for finite-size numerical simulations to validate our general theory.
It should be noted that this paper is quite technical and is accompanied by a
letter focusing on the presentation of the main formulas and their use [14].
In the letter, we also discuss the emergence of a new type of bound polariton
state in the Dicke-Ising model. Here, we delve into all the theoretical
details and provide the aforementioned examples not fully covered in the
letter.

The paper is organized as follows. In Sec. II we present the Hamiltonian and
discuss its validity to describe the typical models of cavity QED materials.
In Sec. III we develop and compare two distinct linear response theories. We
provide expressions for different response functions of the hybrid light-
matter system in terms of the bare response functions of the cavity and the
material. Section IV is dedicated to the application of the LRT to the quantum
Hall effect. More precisely to a 2D electron gas subject a homogeneous
perpendicular classical magnetic field. In Sec. V we apply the theory to
several magnetic models: the Dicke model, the LMG model with longitudinal
coupling to a cavity, the Dicke-LMG model, the Dicke-Ising model and the
Heisenberg model. We end the paper with some Conclusions and relegate some
technical details to the Appendices.

##  II Models in cavity QED materials

We will consider a Hamiltonian of the form

$H=H_{\rm m}+\Omega
a^{\dagger}a+g\left(a+a^{\dagger}\right)C_{x}+\zeta\frac{g^{2}}{\Omega}C_{x}^{2}\,.$
(1)

Here the bare material system is described by $H_{\rm m}$ and remains
unspecified. Likewise $C_{x}$ is an unspecified Hermitian matter coupling
operator. The light matter coupling constant $g=\lambda/\sqrt{N}$ represents
the coupling per particle, with $\lambda$ the collective coupling. The last
term is an optional $P^{2}$ term, toggled by $\zeta\in\\{0,1\\}$.

Depending on the choice of $H_{\rm m}$, $C_{x}$, $\Omega$, $g$ and $\zeta$
this Hamiltonian can describe several distinct microscopic systems. The Pauli-
Fierz Hamiltonian can be brought to the form of Hamiltonian (1) by invoking
the long-wavelength approximation and diagonalizing the photonic sector with a
Bogoliubov transformation to eliminate the $A^{2}$ term [15, 7, 8]. It would
correspond to setting $\zeta=0$ and associating the electronic collective
momenta to $C_{x}$ and renormalized frequency and couplings, stemming from the
Bogoliubov transform, to $\Omega$ and $g$ (See Sec. IV for an example). Fermi-
Hubbard models cannot be cast to the form of Hamiltonian (1) due to the non-
linear nature of the Peierls coupling, unless it is linearized by considering
the weak coupling limit or by moving to the dipole gauge [16, 17, 18].
Localized electric dipoles can be described by setting $\zeta=1$ and
associating the collective dipole operator to $C_{x}$ [19, 8]. Since the
inclusion of the $P^{2}$ term to ensure gauge invariance is a relatively
recent mandate and many common models do not incorporate it, toggling off
$\zeta$ will clarify its effects within linear response theory [20].
Furthermore, we consider $\zeta$ a binary parameter here, but promoting it to
take arbitrary values also allows to study a model where the $P^{2}$ term
amplitude is not fixed by the light-matter coupling and cavity frequency. To
describe magnetic dipoles instead, one can associate a collective spin
operator to $C_{x}$ and toggle off the $P^{2}$ term by setting $\zeta=1$ [21,
22, 23]. Finally, a system of coupled oscillators is also compatible with
Hamiltonian (1), although solving it with the theory described in this paper
would be overkill.

##  III Linear response theory

Throughout this text, we will compute linear response functions, i.e. retarded
Green functions, for the light and matter subsystems. In our derivation, we
will not specify the matter subsystem; instead, we will derive relationships
between the Green functions for the hybrid light-matter system with the bare
matter and photonic Green functions. The bare photonic Green function is a
well know quantity that is trivial to compute. A full solution will thus be
contingent on being able to compute the Green function of the bare matter,
possibly subject to some additional mean fields, as we will gather from the
derivation. The retarded Green function for operators $A$ and $B$ is defined
as

$G^{\rm
r}_{A,B}(t,t^{\prime})=-i\theta(t-t^{\prime})\langle[A(t),B(t^{\prime})]\rangle\,.$
(2)

We will be particularly concerned with the photonic propagator

$D(t)=G^{\rm r}_{a,a^{\dagger}}(t,0)$ (3)

and matter response functions

$\chi_{ab}(t)=-\frac{1}{N}G^{\rm r}_{C_{a},C_{b}}(t,0)\,,$ (4)

where $C_{a}$ might be $C_{x}$ or any other collective matter operator.

###  III.1 Imaginary-time path integral: Thermal Green function

One way to compute retarded Green functions is to obtain them by analytic
continuation from the thermal Green function, which can be computed from an
imaginary-time path integral formulation of the partition function of the
system. The thermal Green function for operators $A$ and $B$ reads

$G^{\rm
t}_{A,B}(\tau-\tau^{\prime})=-\langle\mathcal{T}_{\tau}A(\tau)B(\tau^{\prime})\rangle$
(5)

Within the imaginary-time formalism, we will obtain the thermal Green
functions in terms of imaginary time, $\tau$, or Matsubara frequencies,
$\omega_{m}$. The retarded Green function is then obtained with the
replacement

$G^{\rm r}_{A,B}(\omega)=G^{\rm
t}_{A,B}(\omega_{m})|_{i\omega_{m}\to\omega_{+}=\omega+i0^{+}}\;,$ (6)

for all positive Matsubara frequencies [1][Chap. 7].

####  III.1.1 Effective action

We formulate the partition function as a path integral over coherent states

$Z=\oint_{a,\bar{a},c}e^{-S}\,,$ (7)

with

$\begin{split}S=&S_{\rm
m}+\int_{\tau}\bar{a}(\tau)(\partial_{\tau}-\Omega)a(\tau)\\\
&+g\int_{\tau}(a(\tau)+\bar{a}(\tau))C_{x}(\tau)+\zeta\frac{g^{2}}{\Omega}\int_{\tau}C_{x}^{2}(\tau)\,.\end{split}$
(8)

Here $\oint_{a,\bar{a},c}\equiv\int\mathcal{D}[a,\bar{a}]\mathcal{D}[c]$ is
the functional integral over coherent-state paths and
$\int_{\tau}\equiv\int_{0}^{\beta}d\tau$ is the integral over imaginary time.
A partial integration over the photonic fields, which is but a collection of
Gaussian integrals, yields an effective action $S_{\rm eff}=S_{\rm m}+S_{\rm
ind}$ with an induced retarded collective interaction [24, 25]

$S_{\rm
ind}=g^{2}\int_{\tau,\tau^{\prime}}C_{x}(\tau)D_{0}(\tau-\tau^{\prime})C_{x}(\tau^{\prime})+\zeta\frac{g^{2}}{\Omega}\int_{\tau}C_{x}^{2}(\tau)\,.$
(9)

Here $D_{0}$ is the free photon propagator

$D_{0}(\tau)=-\frac{e^{-\tau\Omega}}{1-e^{-\beta\Omega}}\,,$ (10)

or in Matsubara frequency space

$D_{0}(\omega_{m})=\frac{1}{i\omega_{m}-\Omega}\,.$ (11)

It is convenient to write

$S_{\rm ind}=\frac{1}{2}\int_{\tau,\tau^{\prime}}C_{x}(\tau)\frac{1}{N}V_{\rm
ind}(\tau-\tau^{\prime})C_{x}(\tau^{\prime})\,,$ (12)

with

$V_{\rm
ind}(\tau)=\lambda^{2}(D_{0}(\tau)+D_{0}(\beta-\tau))+2\zeta\frac{\lambda^{2}}{\Omega}\delta(\tau)\,,$
(13)

or in Matsubara frequency space

$V_{\rm
ind}(\omega_{m})=2\lambda^{2}\frac{\Omega^{2}(\zeta-1)+\zeta\omega_{m}^{2}}{\Omega(\omega_{m}^{2}+\Omega^{2})}\,.$
(14)

Note that we symmetrize $V_{\rm ind}$ with respect to $\tau$. It can be seen
that the odd part of $V_{\rm ind}$ does not contribute to the action. This
will be justified later, when we relate $V_{\rm ind}$ to the propagator of an
auxiliary real scalar field, which must be even.

####  III.1.2 Auxiliary field decoupling

The induced interaction can be decoupled with a Hubbard-Stratonovich
transformation that introduces an auxiliary scalar field, $\varphi$:

$\begin{split}e^{-S_{\mathrm{ind}}}=\frac{1}{Z_{\varphi}}\oint_{\varphi}\exp\Bigl{(}&-\frac{1}{2}N\int_{\tau,\tau^{\prime}}\varphi(\tau)V_{\mathrm{ind}}^{-1}(\tau-\tau^{\prime})\varphi(\tau^{\prime})\\\
&-i\int_{\tau}\varphi(\tau)C_{x}(\tau)\Bigr{)}\,.\end{split}$ (15)

We define the propagator of the auxiliary field as

$W(\tau)=\langle\varphi(\tau)\varphi(0)\rangle^{\rm c}\,.$ (16)

The generating functional for bare connected correlation functions of the
matter coupling operator, $C_{x}$, is $\mathcal{G}_{\rm
m}^{0}[\xi]=-N^{-1}\log Z_{\rm m}[\xi]$ with

$Z_{\rm m}[\xi]=\oint_{c}e^{-(S_{\rm m}+i\int_{\tau}\xi(\tau)C_{x}(\tau))}\,,$
(17)

such that

$\begin{split}\frac{\delta}{\delta\xi(\tau_{1})}\cdots&\frac{\delta}{\delta\xi(\tau_{n})}\mathcal{G}_{\rm
m}^{0}[\xi]\Bigr{|}_{\xi=0}=\\\ &=-\frac{(-i)^{n}}{N}\langle
C_{x}(\tau_{1})\cdots C_{x}(\tau_{n})\rangle^{\rm c}_{\rm m}\\\
&\equiv\chi_{xx,0}^{(n)}(\tau_{1},\ldots,\tau_{n})\,.\end{split}$ (18)

In the following we will denote $\chi_{xx,0}^{(2)}\equiv\chi_{xx,0}$. With
this, after partial integration over the matter degrees of freedom, we can
write [26]

$Z=\oint_{\varphi}e^{-Nf[\varphi]}\,,$ (19)

with

$f[\varphi]=\frac{1}{2}\int_{\tau,\tau^{\prime}}\varphi(\tau)V_{\mathrm{ind}}^{-1}(\tau-\tau^{\prime})\varphi(\tau^{\prime})+\mathcal{G}_{\rm
m}^{0}[\varphi]\,.$ (20)

The exponent is proportional to $N$, suggesting that it can be treated with
the saddle-point method for large $N$. The number of terms required for a good
approximation will be determined by the value of $N$. In the thermodynamic
limit, $N\to\infty$, it will be justified to truncate the expansion at the
quadratic term, which will reveal what is the free propagator of the auxiliary
field $\varphi$. The condition that the functional derivative of $f$ with
respect to $\varphi$ vanishes

$\begin{split}\int_{\tau}^{\prime}V_{\rm
ind}^{-1}(\tau-\tau^{\prime})\varphi_{\rm
sp}(\tau^{\prime})+\frac{\delta\mathcal{G}_{\rm m}^{0}[\varphi_{\rm
sp}]}{\delta\varphi_{\rm sp}(\tau)}&=\\\ \int_{\tau}^{\prime}V_{\rm
ind}^{-1}(\tau-\tau^{\prime})\varphi_{\rm
sp}(\tau^{\prime})+\frac{i}{N}\langle C_{x}(\tau)\rangle_{\varphi_{\rm
sp}}&=0\end{split}$ (21)

or in Matsubara frequency space

$\varphi_{\rm sp}(\omega_{m})=-\frac{i}{N}V_{\rm ind}(\omega_{m})\langle
C_{x}(\omega_{m})\rangle_{\varphi_{\rm sp}}$ (22)

formally defines $\varphi_{\rm sp}$. Let us postpone the issue of finding
$\varphi_{\rm sp}$ for now and continue on by expanding $f$ in powers of
$\varphi$ around $\varphi_{\rm sp}$. It is convenient to note that, by
definition,

$\mathcal{G}_{\rm
m}^{0}[\xi]=\sum_{n_{\rm{even}}}\frac{1}{n!}\int_{\tau_{1},\ldots,\tau_{n}}\xi\left(\tau_{1}\right)\cdots\xi\left(\tau_{n}\right)\chi_{xx,0}^{(n)}\left(\tau_{1},\ldots,\tau_{n}\right)\,.$
(23)

Odd terms are absent because the model is invariant under a sign flip of
$C_{x}$. Therefore

$\begin{split}f[\varphi]-f[\varphi_{\rm
sp}]=&\frac{1}{2}\int_{\tau,\tau^{\prime}}\delta\varphi(\tau)(NW_{0})^{-1}(\tau-\tau^{\prime})\delta\varphi(\tau^{\prime})\\\
&+\mathcal{O}(\varphi^{4})\,.\end{split}$ (24)

Here $\delta\varphi=\varphi-\varphi_{\rm sp}$ and

$(NW_{0})^{-1}=V_{\rm ind}^{-1}+\tilde{\chi}_{xx,0}\,,$ (25)

with

$\begin{split}\tilde{\chi}_{xx,0}(\tau,\tau^{\prime})=&\frac{\delta}{\delta\varphi_{\rm
sp}(\tau)}\frac{\delta}{\delta\varphi_{\rm sp}(\tau^{\prime})}\mathcal{G}_{\rm
m}^{0}[\varphi_{\rm sp}]\\\ =&\frac{1}{N}\langle
C_{x}(\tau)C_{x}(\tau^{\prime})\rangle_{\varphi_{\rm{sp}}}^{\rm
c}\,.\end{split}$ (26)

Note that if $\varphi_{\rm sp}=0$ then $\tilde{\chi}_{xx,0}=\chi_{xx,0}$. At
this point we note that $W_{0}$, which is the free propagator of $\varphi$, a
real scalar field, is a function of $V_{\rm ind}$. It now clear why we
symmetrize $V_{\rm ind}$, $W_{0}$ must be even in time. The importance of Eq.
(25) stems from the fact that for $N\to\infty$ one can safely truncate the
saddle-point expansion to second order, which implies that $W=W_{0}$. We can
now obtain relations between the auxiliary field propagator, $W$, the photonic
propagator, $D$, and matter response functions.

####  III.1.3 Dressed response functions in terms of bare response functions

In the following we will define generating functionals for connected
correlators in imaginary time. These differ from the thermal Green function as
defined in Eq. (5) in symmetry broken phases when the expectation value of the
operator for which one is computing the connected correlator acquires a finite
value. However, this difference only affects the zero frequency component of
the correlator and is thus unimportant for the analytic continuation to real
frequencies.

Let us define the generating functional for photonic connected correlators by
introducing a complex field in Eq. (7): $\mathcal{G}_{\rm
ph}[\eta,\bar{\eta}]=-\log Z[\eta,\bar{\eta}]$, with

$Z[\eta,\bar{\eta}]=\oint_{a,\bar{a},c}e^{-(S+\int_{\tau}(a(\tau)\bar{\eta}(\tau)+\bar{a}(\tau)\eta(\tau))}\,.$
(27)

After partial integration of the cavity fields, this yields

$\begin{split}Z[\eta,\bar{\eta}]=\oint_{c}\exp-\Bigl{(}&S_{\rm
m}+\zeta\frac{g^{2}}{\Omega}\int_{\tau}C_{x}^{2}(\tau)\\\
&+\int_{\tau,\tau^{\prime}}\bar{m}(\tau)D_{0}(\tau-\tau^{\prime})m(\tau)\Bigr{)}\,,\end{split}$
(28)

with $m(\tau)=\eta(\tau)+gC_{x}(\tau)$. With this

$\begin{split}D(\tau-\tau^{\prime})=&\left.\frac{\delta}{\delta\bar{\eta}(\tau)}\frac{\delta}{\delta\eta(\tau^{\prime})}\mathcal{G}_{\rm
ph}[\eta,\bar{\eta}]\right|_{\eta=\bar{\eta}=0}\\\
=&D_{0}(\tau-\tau^{\prime})\\\
&-\lambda^{2}\int_{u,v}D_{0}(\tau-u)\chi_{xx}(u-v)D_{0}(v-\tau^{\prime})\,,\end{split}$
(29)

or in Matsubara frequency space

$D(\omega_{m})=D_{0}(\omega_{m})-\lambda^{2}D_{0}(\omega_{m})\chi_{xx}(\omega_{m})D_{0}(\omega_{m})\,.$
(30)

Likewise, we can define the generating functional for matter connected
correlators. Since we are interested in response functions, i.e. two point-
correlators, we can rather generally consider correlations of three distinct
matter operators, the matter coupling operator, $C_{x}$, and two others, that
we label $C_{y}$ and $C_{z}$. Note that $x,y$ and $z$ are simply indices here
and they do not necessarily indicate spatial direction. The generating
functional is $\mathcal{G}_{\rm m}[\xi_{x},\xi_{y},\xi_{z}]=-N^{-1}\log
Z[\xi_{x},\xi_{y},\xi_{z}]$, with

$Z[\xi_{x},\xi_{y},\xi_{z}]=\oint_{c,\varphi}e^{-(S_{\rm
eff}+i\int_{\tau}\sum_{a}\xi_{a}(\tau)C_{a}(\tau))}\,.$ (31)

After partial integration over the cavity fields, this yields

$Z[\xi_{x},\xi_{y},\xi_{z}]=\oint_{\varphi}e^{-Nf[\xi_{x}+\varphi,\xi_{y},\xi_{z}]}\,,$
(32)

with

$\begin{split}f[\xi_{x}+\varphi,\xi_{y},\xi_{z}]=&\frac{1}{2}\int_{\tau,\tau^{\prime}}\varphi(\tau)V_{\mathrm{ind}}^{-1}(\tau-\tau^{\prime})\varphi(\tau^{\prime})\\\
&+\mathcal{G}_{\rm m}^{0}[\xi_{x}+\varphi,\xi_{y},\xi_{z}]\,,\end{split}$ (33)

Here $\mathcal{G}_{\rm m}^{0}[\xi_{x},\xi_{y},\xi_{z}]=-N^{-1}\log Z_{\rm
m}[\xi_{x},\xi_{y},\xi_{z}]$ is the generalization of $\mathcal{G}_{\rm
m}^{0}[\xi]$ (17), with

$Z_{\rm m}[\xi_{x},\xi_{y},\xi_{z}]=\oint_{c}e^{-(S_{\rm
m}+i\int_{\tau}\sum_{a}\xi_{a}(\tau)C_{a}(\tau))}\,.$ (34)

Then a second-order expansion of $f[\xi_{x}+\varphi,\xi_{y},\xi_{z}]$ around
$\varphi_{\rm sp}$ yields

$\begin{split}f[\xi_{x}+\varphi,\xi_{y},\xi_{z}]-f[\varphi_{\rm
sp},0,0]=\frac{1}{2}\int_{\tau,\tau^{\prime}}\Bigl{(}&\delta\varphi(\tau)(NW)^{-1}(\tau-\tau^{\prime})\delta\varphi(\tau^{\prime})+\sum_{a,b}\xi_{a}(\tau)\tilde{\chi}_{ab,0}(\tau-\tau^{\prime})\xi_{b}(\tau^{\prime})\\\
&+\,\delta\varphi(\tau)\sum_{a}\tilde{\chi}_{xa,0}(\tau-\tau^{\prime})\xi_{a}(\tau^{\prime})+\sum_{a}\xi_{a}(\tau)\tilde{\chi}_{ax,0}(\tau-\tau^{\prime})\delta\varphi(\tau)\Bigr{)}\,.\end{split}$
(35)

The functional integral over the auxiliary-field displacements $\delta\varphi$
is just a Gaussian integral that we can perform, yielding

$\begin{split}\mathcal{G}_{\rm
m}[\xi_{x},\xi_{y},\xi_{z}]=\rm{cst.}+\frac{1}{2}\sum_{a,b}\int_{\tau,\tau^{\prime}}\xi_{a}(\tau)\Biggl{(}&\tilde{\chi}_{ab,0}(\tau-\tau^{\prime})-\int_{u,v}\tilde{\chi}_{ax,0}(\tau-u)NW(u-v)\tilde{\chi}_{xb,0}(u-\tau^{\prime})\Biggr{)}\xi_{b}(\tau^{\prime})\,.\end{split}$
(36)

With this

$\begin{split}\chi_{ab}(\tau-\tau^{\prime})=\left.\frac{\delta}{\delta\xi_{a}(\tau)}\frac{\delta}{\delta\xi_{b}(\tau^{\prime})}\mathcal{G}_{\rm
m}[\xi_{x},\xi_{y},\xi_{z}]\right|_{\xi_{a}=0\;\forall
a}=\tilde{\chi}_{ab,0}(\tau-\tau^{\prime})-\int_{u,v}\tilde{\chi}_{ax,0}(\tau-u)NW(u-v)\tilde{\chi}_{xb,0}(u-\tau^{\prime})\,,\end{split}$
(37)

or in Matsubara frequency space

$\chi_{ab}(\omega_{m})=\tilde{\chi}_{ab,0}(\omega_{m})-\chi_{ax,0}(\omega_{m})NW(\omega_{m})\chi_{xa,0}(\omega_{m})\,.$
(38)

Putting together Eqs. (25) and (38) we arrive at the relations

$\displaystyle\chi_{xa}(\omega_{m})=\frac{\tilde{\chi}_{xa,0}(\omega_{m})}{1+V_{\rm
ind}(\omega_{m})\tilde{\chi}_{xx,0}(\omega_{m})}\,,$ (39)
$\displaystyle\chi_{ax}(\omega_{m})=\frac{\tilde{\chi}_{ax,0}(\omega_{m})}{1+V_{\rm
ind}(\omega_{m})\tilde{\chi}_{xx,0}(\omega_{m})}\,,$ (40)
$\displaystyle\begin{split}&\left.\chi_{ab}\right|_{a,b\neq
x}(\omega_{m})=\tilde{\chi}_{ab,0}(\omega_{m})\\\
&\qquad\qquad\qquad\quad\;-\frac{\tilde{\chi}_{ax,0}(\omega_{m})V_{\rm
ind}(\omega_{m})\tilde{\chi}_{xb,0}(\omega_{m})}{1+V_{\rm
ind}(\omega_{m})\tilde{\chi}_{xx,0}(\omega_{m})}\,,\end{split}$ (41)

which together with Eq. (30) give us the dependence of the photonic propagator
and matter response functions on the bare matter response function and the
free photon propagator. Note from Eqs. (159) and (34) that
$\tilde{\chi}_{ab,0}$ is simply the response function of the bare matter
subject to an external field $\varphi_{\rm sp}$, the saddle point. Equations
(39), (39) and (41) exhibit a structure similar to the response functions
commonly derived using the random phase approximation, such as those found in
[26]. It is important to highlight that these equations remain valid for all
parameter ranges and are applicable to both normal and superradiant phases.

####  III.1.4 Computing the saddle point $\varphi_{\rm sp}$.

We should now tackle the issue of finding $\varphi_{\rm sp}$. Without
specifying the matter subsystem it is reasonable to assume a constant
solution, i.e. $\varphi_{\rm sp}(\tau)=\varphi_{\rm sp}=\varphi_{\rm
sp}(\omega_{m}=0)$ (all other frequency components being zero). Then, from Eq.
22 we find

$\varphi_{\rm sp}=-\frac{i}{N}V_{\rm ind}(\omega_{m}=0)\langle
C_{x}\rangle_{\varphi_{\rm sp}}\,.$ (42)

This tells us that $\varphi_{\rm sp}$ is self-consistently proportional to
$\langle C_{x}\rangle_{\varphi_{\rm sp}}$, i.e. to the expectation value of
$C_{x}$ for the bare matter subject to a field
$i\varphi_{\rm{sp}}=\frac{1}{N}V_{\rm ind}(\omega_{m}=0)\langle C_{x}\rangle$.
This is precisely the self-consistent condition that arises from computing
$\langle C_{m}\rangle$ from the mean-field Hamiltonian

$H_{\rm eff}^{\rm MF}=H_{\rm m}+\frac{2\lambda^{2}(\zeta-1)}{N\Omega}\langle
C_{x}\rangle C_{x}-\frac{\lambda^{2}(\zeta-1)}{N\Omega}\langle
C_{x}\rangle^{2}\,.$ (43)

This is the mean-field theory of the effective Hamiltonian that arises from
taking the static limit in $V_{\rm ind}$ in the effective action [8],
justifying (42). Thus, if we assume $\varphi_{\rm sp}$ constant,
$\tilde{\chi}_{ab,0}$ are just the response functions of the mean field matter
Hamiltonian of Eq. (43). Thus, $\tilde{\chi}_{ab,0}$ are obtained as functions
of $\langle C_{x}\rangle$ which can be computed by solving $H_{\rm eff}^{\rm
MF}$ variationaly with respect to $\langle C_{x}\rangle$. Note that in the
absence of symmetry breaking, i.e. when $\langle C_{x}\rangle=0$, or if
$\zeta=1$, the effective matter Hamiltonian is just the bare matter
Hamiltonian and thus $\tilde{\chi}_{ab,0}=\chi_{ab,0}$ is just the bare matter
response.

###  III.2 Equations of motion for the retarded Green functions

An alternative way to compute retarded Green functions is by formulating their
equations of motion and, under some assumptions, solving them. Defining
$G^{\rm r}_{A,B}(t)\equiv G^{\rm r}_{A,B}(t,0)$ we have from Eq. (2) that

$i\partial tG^{\rm r}_{A,B}(t)=\delta(t)\langle[A,B]\rangle-G^{\rm
r}_{[H,A],B}(t)$ (44)

or in frequency space

$\omega_{+}G^{\rm r}_{A,B}(\omega)=\langle[A,B]\rangle-G^{\rm
r}_{[H,A],B}(\omega)\,,$ (45)

with $\omega_{+}=\omega+i0^{+}$.

####  III.2.1 Solving the equations of motion

Let $\\{|\alpha\rangle\\}$ be a basis of the Hilbert space of the matter
subsystem and let us define the Hubbard operators
$X^{\alpha\beta}=|\alpha\rangle\langle\beta|$. We start by computing the
matter correlation function $\chi_{ab}$. Note that

$\chi_{ab}=-\frac{1}{N}G^{\rm
r}_{C_{a},C_{b}}=-\frac{1}{N}\sum_{\alpha,\beta}C_{a}^{\alpha\beta}G^{\rm
r}_{X^{\alpha\beta},C_{b}}$ (46)

where $C_{a}^{\alpha\beta}=\langle\alpha|C_{a}|\beta\rangle$. Then

$\begin{split}\omega_{+}G^{\rm
r}_{X^{\alpha\beta},C_{b}}(\omega)=&\langle[X^{\alpha\beta},C_{b}]\rangle-G^{\rm
r}_{[H,X^{\alpha\beta}],C_{b}}(\omega)\\\ =&-\Gamma_{b}^{\alpha\beta}-G^{\rm
r}_{[H_{\rm m},X^{\alpha\beta}],C_{b}}(\omega)\\\ &-gG^{\rm
r}_{(a+a^{\dagger})[C_{x},X^{\alpha\beta}],C_{b}}(\omega)\\\
&-\zeta\frac{g^{2}}{\Omega}G^{\rm
r}_{[C_{x}^{2},X^{\alpha\beta}],C_{b}}(\omega)\,.\end{split}$ (47)

Here
$\Gamma_{b}^{\alpha\beta}=\langle[C_{b},X^{\alpha\beta}]\rangle=\sum_{\mu}(C_{b}^{\mu\alpha}\langle
X^{\mu\beta}\rangle-C_{b}^{\beta\mu}\langle X^{\alpha\mu}\rangle)$. The last
two terms of Eq. (47) correspond to three-point Green functions. One could now
write the equations of motion for such three-point Green functions, which
would in turn depend on four-point Green functions, and so on and so forth.
This infinite hierarchy of $n$-point correlators involves light-matter
interactions, as discussed in the context of lasing physics [27, Chapter 9] or
light emission [28]. At some point, the hierarchy needs to be truncated. The
crudest approximation is to restrict the theory to two-point Green functions,
assuming that $G^{\rm r}_{AB,C}\approx\langle A\rangle G^{\rm r}_{B,C}+\langle
B\rangle G^{\rm r}_{A,C}$, which amounts to neglecting correlations of order
higher than two for all terms stemming from light-matter interaction. In terms
of fluctuations, this means that only linear fluctuations between the two
systems are relevant, and that correlated fluctuations are largely suppressed.
This approach is reminiscent of RPA or mean-field approximations, which need
to be justified. This has indeed been discussed within the Dicke model, where
the traditional educated guess that quantum fluctuations vanish in the
$N\to\infty$ limit [27, Chapter 9] [Cf. with the equilibrium argument in [6]]
has been rigorously demonstrated recently [29]. In that work, the proof was
restricted to the Dicke model and within the Lindblad framework, _i.e._ ,
adding dissipation in a Markovian manner. This result is relevant to us since
LRT falls within the conditions of the theorem demonstrated in that reference,
which requires the dynamics to start with an uncorrelated initial condition
between light and matter. This can be safely assumed here since the response
is independent of the excitation and can be obtained from excitation functions
where the initial condition is assumed to be at equilibrium, which can be well
approximated by uncorrelated light-matter states [7, 23]. It is true, however,
that in our case we are not considering dissipative Lindblad dynamics and we
include potential matter-matter interactions. In any case, we believe that
restricting to two-point Green functions can be proven to be exact, as it
agrees with the exactness of the saddle-point solution for $N\to\infty$
discussed in the previous section. Additionally, this section concludes by
obtaining the same response functions as within the path integral formalism.
This suggests that generalizing Ref. [29] might be possible in a future work.
With this,

$\begin{split}\omega_{+}G^{\rm
r}_{X^{\alpha\beta},C_{b}}(\omega)=&-\Gamma_{b}^{\alpha\beta}-G^{\rm
r}_{[H_{\rm m},X^{\alpha\beta}],C_{b}}(\omega)\\\ &-g\langle
a+a^{\dagger}\rangle G^{\rm r}_{[C_{x},X^{\alpha\beta}],C_{b}}(\omega)\\\
&-\zeta\frac{2g^{2}}{\Omega}\langle C_{x}\rangle G^{\rm
r}_{[C_{x},X^{\alpha\beta}],C_{b}}(\omega)\\\
&-g\Gamma_{x}^{\alpha\beta}\left(G^{\rm r}_{a,C_{b}}(\omega)+G^{\rm
r}_{a^{\dagger},C_{b}}(\omega)\right)\\\
&-\zeta\frac{2g^{2}}{\Omega}\Gamma_{x}^{\alpha\beta}G^{\rm
r}_{C_{x},C_{b}}(\omega)\,.\end{split}$ (48)

We can put together the second, third and fourth elements in the right-hand
side to build a mean field matter Hamiltonian

$H_{\rm m}^{\rm MF}=H_{\rm m}-g\langle a+a^{\dagger}\rangle
C_{x}-\zeta\frac{2g^{2}}{\Omega}\langle C_{x}\rangle C_{x}\,.$ (49)

Note that the constant terms in this mean-field Hamiltonian are omitted here.
We will show later that this Hamiltonian is equivalent to the mean-field
effective Hamiltonian (43) obtained with the path integral approach under a
mean-field decoupling of full light-matter Hamiltonian. At this point we can
assume that the basis that we have been using, $\\{|\alpha\rangle\\}$, is
precisely the eigenbasis of $H_{\rm m}^{\rm MF}$, with eigenenergies
$\\{E_{\alpha}\\}$, obtaining

$\begin{split}\omega_{+}G^{\rm
r}_{X^{\alpha\beta},C_{b}}(\omega)=&-\Gamma_{b}^{\alpha\beta}-(E_{\alpha}-E_{\beta})G^{\rm
r}_{X^{\alpha\beta},C_{b}}(\omega)\\\ &-g\Gamma_{x}^{\alpha\beta}\left(G^{\rm
r}_{a,C_{b}}(\omega)+G^{\rm r}_{a^{\dagger},C_{b}}(\omega)\right)\\\
&-\zeta\frac{2g^{2}}{\Omega}\Gamma_{x}^{\alpha\beta}G^{\rm
r}_{C_{x},C_{b}}(\omega)\,.\end{split}$ (50)

We can compute the photon-matter Green functions that appear on the second
line obtaining

$\displaystyle(\omega_{+}-\Omega)G^{\rm r}_{a,C_{b}}(\omega)$
$\displaystyle=gG^{\rm r}_{C_{x},C_{b}}(\omega)\,,$ (51)
$\displaystyle(\omega_{+}+\Omega)G^{\rm r}_{a^{\dagger},C_{b}}(\omega)$
$\displaystyle=-gG^{\rm r}_{C_{x},C_{b}}(\omega)\,.$ (52)

Putting together Eqs. (50), (51) and (52) we get

$(\omega_{+}+E_{\alpha}-E_{\beta})G^{\rm
r}_{X^{\alpha\beta},C_{b}}(\omega)=-\Gamma_{b}^{\alpha\beta}-\Gamma_{x}^{\alpha\beta}V_{\rm
ind}(\omega)\frac{1}{N}G^{\rm r}_{C_{x},C_{b}}(\omega)\,,$ (53)

with $V_{\rm ind}(\omega)=V_{\rm
ind}(\omega_{m})|_{i\omega_{m}\to\omega+i0^{+}}$. Together with Eq. (46) this
gives

$\begin{split}\chi_{ab}(\omega)=\frac{1}{N}\sum_{\alpha\beta}C_{a}^{\alpha\beta}\frac{\Gamma_{b}^{\alpha\beta}}{\omega_{+}+E_{\alpha}-E_{\beta}}-\frac{1}{N}\sum_{\alpha\beta}C_{a}^{\alpha\beta}\frac{\Gamma_{x}^{\alpha\beta}}{\omega_{+}+E_{\alpha}-E_{\beta}}V_{\rm
ind}(\omega)\chi_{xb}(\omega)\,.\end{split}$ (54)

From the definition of $\Gamma_{b}^{\alpha\beta}$ and noting that $\langle
X^{\alpha\beta}\rangle=Z^{-1}\delta_{\alpha\beta}\exp(-\beta E_{\alpha})$ we
find that

$\frac{1}{N}\sum_{\alpha\beta}C_{a}^{\alpha\beta}\frac{\Gamma_{b}^{\alpha\beta}}{\omega_{+}+E_{\alpha}-E_{\beta}}=\frac{1}{N}\frac{1}{Z}\sum_{\alpha\beta}\frac{e^{-\beta
E_{\beta}}-e^{-\beta
E_{\alpha}}}{\omega_{+}+E_{\alpha}-E_{\beta}}\langle\alpha|C_{a}|\beta\rangle\langle\beta|C_{b}|\alpha\rangle=\tilde{\chi}_{ab,0}(\omega)\,.$
(55)

Note here the abuse of notation by using $\beta$ both as an index and as
inverse temperature. With this and Eq. (54), we finally obtain

$\displaystyle\chi_{xa}(\omega)=\frac{\tilde{\chi}_{xa,0}(\omega)}{1+V_{\rm
ind}(\omega)\tilde{\chi}_{xx,0}(\omega)}\,,$ (56)
$\displaystyle\chi_{ax}(\omega)=\frac{\tilde{\chi}_{ax,0}(\omega)}{1+V_{\rm
ind}(\omega)\tilde{\chi}_{xx,0}(\omega)}\,,$ (57)
$\displaystyle\begin{split}&\left.\chi_{aa}\right|_{a,b\neq
x}(\omega)=\tilde{\chi}_{ab,0}(\omega)\\\
&\qquad\qquad\qquad\;\;\;-\frac{\tilde{\chi}_{ax,0}(\omega)V_{\rm
ind}(\omega)\tilde{\chi}_{xb,0}(\omega)}{1+V_{\rm
ind}(\omega)\tilde{\chi}_{xx,0}(\omega)}\,,\end{split}$ (58)

which coincide with the analytic continuation of Eqs. (39), (40) and (41).

Proceeding in a similar fashion, we can now compute $D$ and
$D_{+}=G_{a^{\dagger},a^{\dagger}}^{\rm r}$. First we have

$\displaystyle(\omega_{+}-\Omega)D(\omega)$ $\displaystyle=1+gG^{\rm
r}_{C_{x},a^{\dagger}}(\omega)\,,$ (59)
$\displaystyle(\omega_{+}+\Omega)D_{+}(\omega)$ $\displaystyle=-gG^{\rm
r}_{C_{x},a^{\dagger}}(\omega)\,.$ (60)

Again we use the decomposition $G^{\rm
r}_{C_{x},a^{\dagger}}=\sum_{\alpha\beta}C_{x}^{\alpha\beta}G^{\rm
r}_{X^{\alpha\beta},a^{\dagger}}$ and compute

$\begin{split}(\omega_{+}+E_{\alpha}-E_{\beta})G^{\rm
r}_{X^{\alpha\beta},a^{\dagger}}=&-g\Gamma_{b}^{\alpha\beta}\left(D(\omega)+D_{+}(\omega)\right)\\\
&-\zeta\frac{2g^{2}}{\Omega}\Gamma_{b}^{\alpha\beta}G^{\rm
r}_{C_{x},a^{\dagger}}(\omega)\,,\end{split}$ (61)

such that

$G^{\rm
r}_{C_{x},a^{\dagger}}=\frac{-Ng\tilde{\chi}_{xx,0}(\omega)}{1+\zeta\frac{2\lambda^{2}}{\Omega}\tilde{\chi}_{xx,0}(\omega)}\left(D(\omega)+D_{+}(\omega)\right)\,.$
(62)

Plugging this into Eqs. (59) and (60) we get

$\displaystyle\begin{split}(\omega-\Omega)D(\omega)&=1\\\
&\quad-\frac{\lambda^{2}\tilde{\chi}_{xx,0}(\omega)}{1+\zeta\frac{2\lambda^{2}}{\Omega}\tilde{\chi}_{xx,0}(\omega)}\left(D(\omega)+D_{+}(\omega)\right)\,,\end{split}$
(63) $\displaystyle(\omega+\Omega)D_{+}(\omega)$
$\displaystyle=\frac{\lambda^{2}\tilde{\chi}_{xx,0}(\omega)}{1+\zeta\frac{2\lambda^{2}}{\Omega}\tilde{\chi}_{xx,0}(\omega)}\left(D(\omega)+D_{+}(\omega)\right)\,.$
(64)

We can now solve for $D$ and $D_{+}$ to get

$\displaystyle D(\omega)$
$\displaystyle=\frac{\omega_{+}+\Omega+\lambda^{2}\frac{2\zeta(\omega_{+}+\Omega)-\Omega}{\Omega}\tilde{\chi}_{xx,0}(\omega)}{\omega_{+}^{2}-\Omega^{2}+2\lambda^{2}\frac{\zeta(\omega_{+}^{2}-\Omega^{2})+\Omega^{2}}{\Omega}\tilde{\chi}_{xx,0}(\omega)}\,,$
(65) $\displaystyle D_{+}(\omega)$
$\displaystyle=\frac{\lambda^{2}\tilde{\chi}_{xx,0}(\omega)}{\omega_{+}^{2}-\Omega^{2}+2\lambda^{2}\frac{\zeta(\omega_{+}^{2}-\Omega^{2})+\Omega^{2}}{\Omega}\tilde{\chi}_{xx,0}(\omega)}\,.$
(66)

Although not immediately obvious, Eq. (65) is precisely the analytic
continuation of Eq. (30).

####  III.2.2 Mean field decoupling of the light-matter Hamiltonian

Let us consider the Hamiltonian in Eq. (1) with a mean-field decoupling of the
light-matter interaction and $P^{2}$ terms. We expand $C_{x}=\langle
C_{x}\rangle+\delta C_{x}$ and $a+a^{\dagger}=\langle
a+a^{\dagger}\rangle+\delta(a+a^{\dagger})$ and neglect terms quadratic in
fluctuations, obtaining two decoupled Hamiltonians and some constants:
$H=H_{\rm m}^{\rm{MF}}+H_{\rm{ph}}^{\rm{MF}}+\rm{cst.}$ with

$\displaystyle H_{\rm m}^{\rm{MF}}=H_{\rm m}+g\langle a+a^{\dagger}\rangle
C_{x}-\zeta\frac{2g^{2}}{\Omega}\langle C_{x}\rangle C_{x}\,,$ (67)
$\displaystyle H_{\rm{ph}}^{\rm{MF}}=\Omega
a^{\dagger}a+g(a+a^{\dagger})\langle C_{x}\rangle\,,$ (68) $\displaystyle{\rm
cst}.=-g\langle a+a^{\dagger}\rangle\langle
C_{x}\rangle-\zeta\frac{2g^{2}}{\Omega}\langle C_{x}\rangle^{2}\,.$ (69)

The mean-field photonic Hamiltonian $H_{\rm{ph}}^{\rm{MF}}$ can be
diagonalized with a displacement, $a=b-g/\Omega\langle C_{x}\rangle$, yielding
$H_{\rm{ph}}^{\rm{MF}}=\Omega b^{\dagger}b$, a new constant
$-g^{2}/\Omega\langle C_{x}\rangle^{2}$ and the relation $\langle
a+a^{\dagger}\rangle=-2g/\Omega\langle C_{x}\rangle$. Plucking this relation
into the mean-field matter Hamiltonian, $H_{\rm m}^{\rm{MF}}$, and all the
constants gathered in the different steps yields precisely the mean-field
effective Hamiltonian, $H_{\rm eff}^{\rm MF}$, of Eq. (43). The constants are
important for a potential variational solution.

###  III.3 Comparing the two theories

In this section we have shown the equivalence between two different approaches
to address the linear response regime of cavity QED materials. Although
seemingly very different, the path-integral in imaginary time and the equation
of motion in real time provide exact analytical results in the thermodynamic
limit, $N\to\infty$. The effect of the thermodynamic limit is clearly
identified in the path integral approach when the saddle point method is
invoked. Instead, with the equation of motion method it requires an analysis
of the scaling of the correlated parts of higher-order correlation functions
with $N$. On the other hand, the equations-of-motion approach does not require
the auxiliary field of the Hubbard-Stratonovich transformation, which allows
us to interpret the steps used to calculate its propagator $W$. We can see
that the auxiliary field $\varphi_{\text{sp}}$ plays the same role as the
mean-field background in the equation of motion, which couples to the
fluctuations. Together, the two approaches cement the understanding that mean-
field approaches are exact in the thermodynamic limit for Dicke-like models,
even if matter interactions are considered. This has been studied and
exploited before at equilibrium [6, 30, 7, 23, 8] and in the dynamics [29] and
in Refs. [26, 31] and now here for the LRT. Beyond the LRT, the equations-of-
motion approach opens the possibility of studying real-time dynamics.

Finally, the comparison between the two approaches is not just an academic
exercise. In many cases it can be interesting to go beyond the $N\to\infty$
limit. This would be the case in topological systems, where boundaries play an
important role. In that case the two approaches offer complementary methods
incorporate finite-size corrections. In the path integral approach,
corrections are computed by extending the saddle-point expansion to higher
orders, and thus relaxing $W=W_{0}$ for the weaker $W^{-1}=W_{0}^{-1}+\Pi$,
where $\Pi$ can be computed with standard diagrammatic techniques in terms of
fourth- and higher-order matter response functions [26]. In the equations-of-
motion approach, corrections are introduced by expanding the theory from two-
to three- and higher-order Green functions before invoking the mean-field
factoring of the correlators.

##  IV Integer quantum Hall effect in a cavity

###  IV.1 The Hamiltonian

To demonstrate the applicability of our theory, we will use it to study the
modification of the integer quantum Hall effect. The matter system under
consideration is a two-dimensional electron gas subject to a perpendicular
classical magnetic field. This model has been used to explain the breakdown of
topological protection that underlies the observed modifications of the
quantum Hall effect by coupling to a cavity [32, 11, 33, 34]. The Hamiltonian
of the electron gas coupled to a single-mode cavity reads

$\begin{split}H=&\sum_{j}\frac{(\boldsymbol{p}_{j}-e\boldsymbol{A}_{\rm
ext}(\boldsymbol{r}_{j})-e\boldsymbol{A}_{0}(b+b^{\dagger}))^{2}}{2m}\\\
&+\sum_{i>j}V(\boldsymbol{r}_{i}-\boldsymbol{r}_{j})+\Omega
b^{\dagger}b\,.\end{split}$ (70)

Here $\boldsymbol{r}_{i}$, $\boldsymbol{p}_{i}$ are respectively the position
and momentum operators of the $i$-th electron. Both the external classical
vector potential and the cavity’s vector potential point along the
$x$-direction: $\boldsymbol{A}_{\rm
ext}(\boldsymbol{r}_{j})=-By_{j}\boldsymbol{e}_{x}$ and
$\boldsymbol{A}_{0}=\sqrt{1/(2\epsilon_{0}V\Omega)}\boldsymbol{e}_{x}$. The
$V$ and $\epsilon_{0}$ are respectively the cavity mode volume and the
dielectric constant. $V(\boldsymbol{r}_{i}-\boldsymbol{r}_{j})$ is the Coulomb
interaction between the electrons. After some manipulation, we can write the
Hamiltonian as

$\begin{split}H=&H_{\rm m}-\frac{\omega_{\rm
p}}{\sqrt{N}}\left(\sqrt{\frac{\omega_{\rm
p}}{\Omega}}\bar{p}_{x}+\sqrt{\frac{\omega_{\rm
c}}{\Omega}}\bar{y}\right)\left(b+b^{\dagger}\right)\\\ &+\Omega
b^{\dagger}b+\frac{\omega_{\rm
p}^{2}}{4\Omega}\left(b+b^{\dagger}\right)^{2}\,,\end{split}$ (71)

where $H_{\rm m}$ is the bare (without cavity) electron gas Hamiltonian

$H_{\rm m}=\sum_{j}\frac{(\boldsymbol{p}_{j}-e\boldsymbol{A}_{\rm
ext}(\boldsymbol{r}_{j}))^{2}}{2m}+\sum_{i>j}V(\boldsymbol{r}_{i}-\boldsymbol{r}_{j})$
(72)

and $\bar{p}_{x}$ and $\bar{y}$ are adimensionalized collective operators

$\displaystyle\bar{p}_{x}=\frac{1}{\sqrt{2m\omega_{\rm p}}}\sum_{j}p_{x,j}\,,$
(73) $\displaystyle\bar{y}=\sqrt{\frac{m\omega_{\rm c}}{2}}\sum_{j}y_{j}\,.$
(74)

We have defined the plasma frequency, $\omega_{\rm
p}=\sqrt{e^{2}\rho/(m\epsilon_{0})}$ with $\rho=N/V$ the density of electrons
is the cavity, and the cyclotron frequency $\omega_{\rm c}=eB/m$. The photonic
terms of the Hamiltonian can be diagonalized with a Bogoliubov transformation,
which yields

$H=H_{\rm m}-\frac{\omega_{\rm p}}{\sqrt{N}}\left(\sqrt{\frac{\omega_{\rm
p}}{\tilde{\omega}}}\bar{p}_{x}+\sqrt{\frac{\omega_{\rm
c}}{\tilde{\omega}}}\bar{y}\right)\left(a+a^{\dagger}\right)+\tilde{\omega}a^{\dagger}a\,,$
(75)

with $\tilde{\omega}^{2}=\Omega^{2}+\omega_{\rm p}^{2}$. By identifying
$g=-\omega_{\rm p}/\sqrt{N}$ and $C_{x}=\sqrt{\omega_{\rm
p}/\tilde{\omega}}\bar{p}_{x}+\sqrt{\omega_{\rm c}/\tilde{\omega}}\bar{y}$ the
Hamiltonian takes the form of Eq. (1). . Note that in this case there is no
$P^{2}$ term, i.e. $\zeta=0$, as we are working in the Coulomb gauge. In
exchange, the elimination of the $A^{2}$ term through the Bogoliubov transform
renormalizes the cavity frequency and light-matter coupling.

###  IV.2 Current and optical conductivity

In the quantum Hall effect [35, 36], one is typically interested in computing
the longitudinal and transverse conductivities. The optical conductivity
tensor can be computed with the Kubo formalism in terms of the current
response functions

$\sigma_{ab}(\omega)=\frac{i}{\omega_{+}}\left[\frac{e^{2}\rho_{\rm{2D}}}{m}\delta_{ab}+\frac{G^{\rm
r}_{J_{a},J_{b}}(\omega)}{A}\right]\,,$ (76)

where $\omega_{+}=\omega+i\delta$, $A$ and $\rho_{\rm{2D}}$ are the area and
surface density of the two-dimensional electron gas, $\delta$ is the
broadening parameter, $\delta_{ab}$ is the Kronecker delta and $J_{a,b}$ are
the current operators in the $a$ and $b$ directions, with $a,b\in\\{x,y\\}$
[37, 38]. We will now compute $G^{\rm r}_{J_{a},J_{b}}$ using the theory
outlined in Sec. III, i.e. we will express $G^{\rm r}_{J_{a},J_{b}}$ in terms
of the response functions of the bare electron gas, $G^{r0}_{J_{a},J_{b}}$.

The current operator is given by

$\begin{split}\boldsymbol{J}&=e\sum_{j}\boldsymbol{v}_{j}\\\
&=\frac{e}{m}\sum_{j}\left(\boldsymbol{p}_{j}-e\boldsymbol{A}_{\rm
ext}(\boldsymbol{r}_{j})-e\boldsymbol{A}_{0}(b+b^{\dagger})\right)\,.\end{split}$
(77)

Then, after some manipulation, we can express

$\displaystyle
J_{x}=\sqrt{\frac{2e^{2}\tilde{\omega}}{m}}\left(C_{x}+\gamma(a+a^{\dagger})\right)\,,$
(78) $\displaystyle J_{y}=\sqrt{\frac{2e^{2}\tilde{\omega}}{m}}C_{y}\,,$ (79)

with $\gamma=-\sqrt{N}\omega_{\rm p}/(2\tilde{\omega})$ and
$C_{y}=\sqrt{\omega_{\rm p}/\tilde{\omega}}\bar{p}_{y}$. Here
$\sqrt{2m\omega_{\rm p}}\bar{p}_{y}=\sum_{j}p_{y,j}$. Meaning that we need to
compute

$\displaystyle G^{\rm r}_{J_{x},J_{x}}=\frac{2e^{2}\tilde{\omega}}{m}G^{\rm
r}_{C_{x}+\gamma(a+a^{\dagger}),C_{x}+\gamma(a+a^{\dagger})}\,,$ (80)
$\displaystyle G^{\rm r}_{J_{x},J_{y}}=\frac{2e^{2}\tilde{\omega}}{m}G^{\rm
r}_{C_{x}+\gamma(a+a^{\dagger}),C_{y}}\,,$ (81) $\displaystyle G^{\rm
r}_{J_{y},J_{y}}=\frac{2e^{2}\tilde{\omega}}{m}G^{\rm r}_{C_{y},C_{y}}\,.$
(82)

###  IV.3 Implicit expressions for the current response functions

####  IV.3.1 Computing mixed light-matter correlations

To compute these alternative Green functions following our section III, we
formulate the following generating functional for mixed connected correlators:
$\mathcal{G}_{\rm mix}[\eta,\xi_{x},\xi_{y}]=-\log Z[\eta,\xi_{x},\xi_{y}]$,
with

$\begin{split}Z[\eta,\xi_{x},\xi_{y}]=\oint_{a,\bar{a},c}\exp-\Bigl{(}&S+\int_{\tau}\eta(\tau)(a(\tau)+\bar{a}(\tau))\\\
&+\sum_{a}\int_{\tau}\xi_{a}(\tau)C_{a}(\tau)\Bigr{)}\,.\end{split}$ (83)

After partial integration of the cavity fields, this yields

$\begin{split}Z[\eta,\xi_{x},\xi_{y}]=\oint_{c}\exp-\Bigl{(}&S_{\rm
m}+\int_{\tau,\tau^{\prime}}m(\tau)D_{0}(\tau-\tau^{\prime})m(\tau)\\\
&+\sum_{a}\int_{\tau}\xi_{a}(\tau)C_{a}(\tau)\Bigr{)}\,,\end{split}$ (84)

with $m(\tau)=\eta(\tau)+gC_{x}(\tau)$. With this

$\begin{split}G^{\rm
t}_{\gamma(a+a^{\dagger}),C_{a}}(\tau-\tau^{\prime})&=-\gamma\frac{\delta}{\delta\eta(\tau)}\frac{\delta}{\delta\xi_{a}(\tau^{\prime})}\mathcal{G}_{\rm
mix}[\eta,\xi_{x},\xi_{y}]\\\ &=\gamma\int_{u}D_{0}^{\rm
s}(\tau-u)Ng\chi_{xa}(u-\tau^{\prime})\,,\end{split}$ (85)

or in Matsubara frequency space

$G^{\rm t}_{\gamma(a+a^{\dagger}),C_{a}}(\omega_{m})=-\gamma D_{0}^{\rm
s}(\omega_{m})Ng\chi_{xa}(\omega_{m})\,.$ (86)

Here $D_{0}^{s}(\tau)=D_{0}(\tau)+D_{0}(-\tau)$ is the symmetrized free photon
propagator

$D_{0}^{\rm
s}(\omega_{m})=\frac{1}{i\omega_{m}-\tilde{\omega}}-\frac{1}{i\omega_{m}+\tilde{\omega}}\,.$
(87)

Likewise

$G^{\rm t}_{C_{x},\gamma(a+a^{\dagger})}(\omega_{m})=-\gamma D_{0}^{\rm
s}(\omega_{m})Ng\chi_{xx}(\omega_{m})$ (88)

and

$\begin{split}G^{\rm
t}_{\gamma(a+a^{\dagger}),\gamma(a+a^{\dagger})}(\omega_{m})=&\gamma^{2}D_{0}^{\rm
s}(\omega_{m})\\\ &-\gamma^{2}D_{0}^{\rm
s}(\omega_{m})\lambda^{2}\chi_{xx}(\omega_{m})D_{0}^{\rm
s}(\omega_{m})\,.\end{split}$ (89)

With these we obtain, after analytic continuation $i\omega_{m}\to\omega_{+}$,

$\begin{split}G^{\rm
r}_{J_{x},J_{y}}(\omega)&=-N\frac{2e^{2}\tilde{\omega}}{m}\left(1+\gamma
gD_{0}^{\rm s}(\omega)\right)\chi_{xy}(\omega)\,,\\\
&=-N\frac{2e^{2}\tilde{\omega}}{m}\left(1+\left(\frac{\omega_{\rm
p}}{\tilde{\omega}}\right)^{2}\frac{\tilde{\omega}}{2}D_{0}^{\rm
s}(\omega)\right)\chi_{xy}(\omega)\,,\end{split}$ (90)

and similarly

$\displaystyle G^{\rm
r}_{J_{x},J_{x}}(\omega)=-N\frac{2e^{2}\tilde{\omega}}{m}\left(\left(1+\left(\frac{\omega_{\rm
p}}{\tilde{\omega}}\right)^{2}\frac{\tilde{\omega}}{2}D_{0}^{\rm
s}(\omega)\right)^{2}\chi_{xx}(\omega)-\frac{1}{2\tilde{\omega}}\left(\frac{\omega_{\rm
p}}{\tilde{\omega}}\right)^{2}\frac{\tilde{\omega}}{2}D_{0}^{\rm
s}(\omega)\right)\,,$ (91) $\displaystyle G^{\rm
r}_{J_{y},J_{y}}(\omega)=-N\frac{2e^{2}\tilde{\omega}}{m}\chi_{yy}(\omega)\,.$
(92)

We can relate $\chi_{xx}$, $\chi_{xy}$ and $\chi_{yy}$ to
$\tilde{\chi}_{xx,0}$, $\tilde{\chi}_{xy,0}$ and $\tilde{\chi}_{yy,0}$ through
Eqs. (39) and (41), obtaining

$G^{\rm r}_{J_{x},J_{y}}(\omega)=\left(1+\left(\frac{\omega_{\rm
p}}{\tilde{\omega}}\right)^{2}\frac{\tilde{\omega}}{2}D_{0}^{\rm
s}(\omega)\right)\frac{\tilde{G}^{\rm r,0}_{J_{x},J_{y}}(\omega)}{1-V_{\rm
ind}(\omega)\frac{m}{2e^{2}\tilde{\omega}N}\tilde{G}^{\rm
r,0}_{J_{x},J_{x}}(\omega)}\,,$ (93)

with $\tilde{G}^{\rm r,0}_{J_{a},J_{b}}=-N\tilde{\chi}_{ab,0}$. Likewise

$\displaystyle G^{\rm r}_{J_{x},J_{x}}(\omega)=\left(1+\left(\frac{\omega_{\rm
p}}{\tilde{\omega}}\right)^{2}\frac{\tilde{\omega}}{2}D_{0}^{\rm
s}(\omega)\right)^{2}\frac{\tilde{G}^{\rm r,0}_{J_{x},J_{x}}(\omega)}{1-V_{\rm
ind}(\omega)\frac{m}{2e^{2}\tilde{\omega}N}\tilde{G}^{\rm
r,0}_{J_{x},J_{x}}(\omega)}+\frac{Ne^{2}}{m}\left(\frac{\omega_{\rm
p}}{\tilde{\omega}}\right)^{2}\frac{\tilde{\omega}}{2}D_{0}^{\rm
s}(\omega)\,,$ (94) $\displaystyle G^{\rm
r}_{J_{y},J_{y}}(\omega)=\frac{\tilde{G}^{\rm
r,0}_{J_{y},J_{y}}(\omega)+V_{\rm
ind}(\omega)\frac{m}{2e^{2}\tilde{\omega}N}\left(\tilde{G}^{\rm
r,0}_{J_{x},J_{y}}(\omega)\tilde{G}^{\rm
r,0}_{J_{y},J_{x}}(\omega)-\tilde{G}^{\rm
r,0}_{J_{x},J_{x}}(\omega)\tilde{G}^{\rm
r,0}_{J_{y},J_{y}}(\omega)\right)}{1-V_{\rm
ind}(\omega)\frac{m}{2e^{2}\tilde{\omega}N}\tilde{G}^{\rm
r,0}_{J_{x},J_{x}}(\omega)}\,.$ (95)

We will now show that $\tilde{\chi}_{ab,0}=\chi_{ab,0}$ and thus
$\tilde{G}^{\rm r,0}_{J_{a},J_{b}}=G^{\rm r,0}_{J_{a},J_{b}}$.

####  IV.3.2 Proof that $\tilde{\chi}_{ab,0}=\chi_{ab,0}$

From Eqs. (43) and (75) we see that the mean-field effective Hamiltonian of
the electron gas is given by

$\begin{split}H_{\rm eff}^{\rm
MF}=&\sum_{j}\frac{\left(\boldsymbol{p}_{j}+m\omega_{\rm
c}y_{j}\boldsymbol{e}_{x}\right)^{2}}{2m}+\sum_{i>j}V(\boldsymbol{r}_{i}-\boldsymbol{r}_{j})\\\
&-2\frac{\omega_{\rm p}^{2}}{\tilde{\omega}^{2}}\langle p_{x}+m\omega_{\rm
c}y\rangle\sum_{j}\frac{p_{j,x}+m\omega_{\rm c}y_{j}}{2m}\\\
&+\rm{cst.}\,,\end{split}$ (96)

with $\langle p_{x}+m\omega_{\rm c}y\rangle=\langle p_{j,x}+m\omega_{\rm
c}y_{j}\rangle\,\forall j$. Completing squares yields

$\begin{split}H_{\rm eff}^{\rm
MF}&=\sum_{j}\frac{1}{2m}\left(p_{j,x}+m\omega_{\rm c}y_{j}-\frac{\omega_{\rm
p}^{2}}{\tilde{\omega}^{2}}\langle p_{x}+m\omega_{\rm c}y\rangle\right)^{2}\\\
&+\sum_{j}\frac{p_{j,y}^{2}}{2m}+\sum_{i>j}V(\boldsymbol{r}_{i}-\boldsymbol{r}_{j})+\rm{cst.}\,.\end{split}$
(97)

We can see now that

$p_{j,x}+m\omega_{\rm c}y_{j}-\frac{\omega_{\rm
p}^{2}}{\tilde{\omega}^{2}}\langle p_{x}+m\omega_{\rm
c}y\rangle\propto[x_{j},H_{\rm eff}^{\rm MF}]$ (98)

and thus

$\langle p_{x}+m\omega_{\rm c}y\rangle\left(1-\frac{\omega_{\rm
p}^{2}}{\Omega^{2}+\omega_{\rm p}^{2}}\right)=0\,,$ (99)

which implies that $\langle p_{x}+m\omega_{\rm c}y\rangle=0$ and thus $H_{\rm
eff}^{\rm MF}=H_{\rm m}$ and $\tilde{\chi}_{ab,0}=\chi_{ab,0}$. This is a no
go theorem, implying that the equilibrium properties of the electron gas are
not modified by the cavity. Despite this, the linear response properties, are
modified, alternative: as we now discuss.

This implies that $\tilde{G}^{\rm r,0}_{J_{a},J_{b}}=G^{\rm
r,0}_{J_{a},J_{b}}$. With this, Equations (93), (94) and (95) provide us with
an implicit expression of the current response functions of the electron gas
coupled to a cavity in terms of the current response functions of the bare
electron gas. This culminates the application of our linear response theory.
Obtaining explicit expressions for the current response functions is now just
a matter of computing the current response functions of the bare electron gas.
We refer to the later as bare current response functions and we compute them
in the next section.

###  IV.4 Computing the bare current response functions

To study the bare 2D electron gas we lean on the reasoning laid out in Ref.
[34]. There, they study the same 2D electron gas subject to a perpendicular
classical magnetic field and coupled to a cavity (70). They show that the
cavity only couples to the center-of-mass coordinates. Furthermore, they show
that the center-of-mass and relative coordinate sectors of the Hamiltonian
commute. This essentially splits the problem into two, the problem of the
relative coordinates unmodified by the cavity, and the problem of the center
or mass coupled to the cavity. Note from Eqs. (78) and (79) that the currents
also only depend on center-of-mass coordinates. Most importantly, the Coulomb
interaction between electrons only depends on the relative coordinates. In
summary, the full problem of a 2D electron gas coupled to a cavity can be
broken into the problem of the relative coordinates, subject to Coulomb
interactions but decoupled from the cavity, and the problem of an independent
center of mass coordinate, which couples to the cavity and is solely
responsible for the conduction properties of the electron gas. Consequently,
in order to compute bare current response functions, we only need to study the
center-of-mass sector of the matter Hamiltonian

$H_{\rm
m}^{\rm{cm}}=\frac{\left(\boldsymbol{P}-e\boldsymbol{A}_{\rm{ext}}(\boldsymbol{R})\right)^{2}}{2m}\,,$
(100)

with $\sqrt{N}\boldsymbol{R}=\sum_{j}\boldsymbol{r}_{j}$ and
$\sqrt{N}\boldsymbol{P}=\sum_{j}\boldsymbol{p}_{j}$ the center-of-mass
position and momentum operators. This is just the Hamiltonian of a single
electron under a classical magnetic field, the prototipical example of Landau
quantization [39][Chap. 1]. The eigenstates of the system are separable, into
plane waves in the $X$ direction and eigenstates of a displaced harmonic
oscillator of frequency $\omega_{\rm c}$ along the $Y$ direction

$\begin{split}H_{\rm m}^{\rm{cm}}|K_{X},n\rangle&=\omega_{\rm
c}\left(c^{\dagger}c+\frac{1}{2}\right)|K_{X},n\rangle\\\ &=\omega_{\rm
c}\left(n+\frac{1}{2}\right)|K_{X},n\rangle\,,\end{split}$ (101)

With $P_{X}|K_{X},n\rangle=K_{X}|K_{X},n\rangle$. Here

$\displaystyle V=\sqrt{\frac{1}{2\omega_{\rm
c}}}\left(c+c^{\dagger}\right)\,,$ (102) $\displaystyle
P_{V}=i\sqrt{\frac{\omega_{\rm c}}{2}}\left(c^{\dagger}-c\right)\,,$ (103)

and $V=\sqrt{m}(Y+K_{X}l_{c}^{2})$, $P_{V}=P_{Y}/\sqrt{m}$, with
$l_{c}=\sqrt{1/(eB)}$.

From the definition in Eq. (77) (ignoring the contribution from the cavity to
get the bare current) we can see that the current only depends on center-of-
mass coordinates

$\boldsymbol{J}=\sqrt{N}\frac{e}{m}\left(\boldsymbol{P}-e\boldsymbol{A}_{\rm
ext}(\boldsymbol{R})\right)\,.$ (104)

After some manipulation, we can express

$\displaystyle J_{x}=\frac{e\omega_{\rm
c}}{\sqrt{2m}}\left(\frac{1}{\omega_{\rm
c}\sqrt{m}}(P_{X}-K_{X})+\frac{1}{\sqrt{\omega_{\rm
c}}}(c+c^{\dagger})\right)\,,$ (105) $\displaystyle
J_{y}=i\frac{e\sqrt{\omega_{\rm c}}}{\sqrt{2m}}\left(c^{\dagger}-c\right)\,,$
(106)

and the matrix elements

$\displaystyle\begin{split}\langle
K_{X}^{\prime},n|J_{x}|K_{X},m\rangle&=\delta_{K_{X},K_{X}^{\prime}}\frac{e\sqrt{\omega_{\rm
c}}}{\sqrt{2m}}(\sqrt{m+1}\delta_{n,m+1}\\\
&\qquad\qquad\qquad\quad\quad+\sqrt{m}\delta_{n,m-1})\,,\end{split}$ (107)
$\displaystyle\begin{split}\langle
K_{X}^{\prime},n|J_{y}|K_{X},m\rangle&=\delta_{K_{X},K_{X}^{\prime}}i\frac{e\sqrt{\omega_{\rm
c}}}{\sqrt{2m}}(\sqrt{m+1}\delta_{n,m+1}\\\
&\qquad\qquad\qquad\quad\quad-\sqrt{m}\delta_{n,m-1})\,.\end{split}$ (108)

This shows that the current operators are diagonal with respect to the plane
waves along the $X$ direction, such that the spectral decomposition of the
current response functions can be written simply as a sum over eigenstates of
the harmonic oscillator along the $Y$ direction. At zero temperature, this
reads

$G^{\rm r,0}_{J_{a},J_{b}}(\omega)=\sum_{n}\left(\frac{\langle
0|J_{a}|n\rangle\langle
n|J_{b}|0\rangle}{\omega_{+}+E_{0}-E_{n}}-\frac{\langle
n|J_{a}|0\rangle\langle 0|J_{b}|n\rangle}{\omega_{+}+E_{n}-E_{0}}\right)\,.$
(109)

From Eqs. (107), (108) and (109) we obtain the bare current response functions

$\displaystyle\begin{split}&G^{\rm r,0}_{J_{x},J_{x}}(\omega)=G^{\rm
r,0}_{J_{y},J_{y}}(\omega)\\\ &\qquad\qquad\;=-\frac{Ne^{2}\omega_{\rm
c}}{m}\frac{1}{2}\left(\frac{1}{\omega_{+}+\omega_{\rm
c}}-\frac{1}{\omega_{+}-\omega_{\rm c}}\right)\end{split}\,,$ (110)
$\displaystyle G^{\rm r,0}_{J_{x},J_{y}}(\omega)=\frac{Ne^{2}\omega_{\rm
c}}{m}\frac{i}{2}\left(\frac{1}{\omega_{+}+\omega_{\rm
c}}+\frac{1}{\omega_{+}-\omega_{\rm c}}\right)\,,$ (111) $\displaystyle G^{\rm
r,0}_{J_{y},J_{x}}(\omega)=-\frac{Ne^{2}\omega_{\rm
c}}{m}\frac{i}{2}\left(\frac{1}{\omega_{+}+\omega_{\rm
c}}+\frac{1}{\omega_{+}-\omega_{\rm c}}\right)\,.$ (112)

###  IV.5 Explicit expressions for the current response functions

Substituting Eqs. (110), (111), (112) and $V_{\rm ind}(\omega)=\omega_{\rm
p}^{2}D_{0}^{\rm s}(\omega)$, with $D_{0}^{\rm s}(\omega)$ the analytic
continuation of $D_{0}^{\rm s}(\omega_{m})|_{i\omega_{m}\to\omega_{+}}$ (87),
into Eqs. (93), (94) and (95) yields, after some manipulation,

$\displaystyle G^{\rm
r}_{J_{x},J_{y}}(\omega)=i\frac{Ne^{2}}{m}\frac{\omega_{\rm
c}\omega_{+}\left(\omega_{+}^{2}-\Omega^{2}\right)}{\left(\omega_{+}^{2}-\tilde{\omega}^{2}\right)\left(\omega_{+}^{2}-\omega_{\rm
c}^{2}\right)-\omega_{\rm p}^{2}\omega_{\rm c}^{2}}\,,$ (113) $\displaystyle
G^{\rm r}_{J_{x},J_{x}}(\omega)=\frac{Ne^{2}}{m}\frac{(\omega_{\rm
c}^{2}+\omega_{\rm p}^{2})\left(\omega_{+}^{2}-\frac{\Omega^{2}\omega_{\rm
c}^{2}}{\omega_{\rm c}^{2}+\omega_{\rm
p}^{2}}\right)}{\left(\omega_{+}^{2}-\tilde{\omega}^{2}\right)\left(\omega_{+}^{2}-\omega_{\rm
c}^{2}\right)-\omega_{\rm p}^{2}\omega_{\rm c}^{2}}\,,$ (114) $\displaystyle
G^{\rm r}_{J_{y},J_{y}}(\omega)=\frac{Ne^{2}}{m}\frac{\omega_{\rm
c}^{2}\left(\omega_{+}^{2}-\Omega^{2}\right)}{\left(\omega_{+}^{2}-\tilde{\omega}^{2}\right)\left(\omega_{+}^{2}-\omega_{\rm
c}^{2}\right)-\omega_{\rm p}^{2}\omega_{\rm c}^{2}}\,.$ (115)

These current response functions have poles at $\omega_{+}=\pm\Omega_{\pm}$,
with $\Omega_{\pm}$ the frequencies of the Landau polaritons

$2\Omega_{\pm}=\tilde{\omega}^{2}+\omega_{\rm
c}^{2}\pm\sqrt{(\tilde{\omega}^{2}-\omega_{\rm c}^{2})^{2}+4\omega_{\rm
c}^{2}\omega_{\rm p}^{2}}\,.$ (116)

This is all in agreement with the Landau polaritons and current response
functions computed in Ref. [34]. There, they solve the system exactly,
computing the energies and wavefunctions and subsequently the current response
functions from the spectral decomposition formula. Therefore, the agreement
serves as validation that our effective theory correctly predicts the response
functions of a material coupled to a cavity. We also recover the dc
conductivities in the long-wavelength limit $\Omega\to 0$, $\delta\to 0$ of
Ref. [33].

For completeness, let us write down the dc ($\omega\to 0$) conductivities,
which are tipically of interest when considering the quantum Hall effect,

$\displaystyle\sigma_{xy}=\frac{e^{2}\nu}{h}\frac{\omega_{\rm
c}^{2}(\Omega^{2}+\delta^{2})}{(\Omega_{-}^{2}+\delta^{2})(\Omega_{+}^{2}+\delta^{2})}\,,$
(117) $\displaystyle\sigma_{xx}=\sigma_{\rm D}\left(1-\frac{(\omega_{\rm
c}^{2}+\omega_{\rm p}^{2})\left(\frac{\Omega^{2}\omega_{\rm
c}^{2}}{\omega_{\rm c}^{2}+\omega_{\rm
p}^{2}}+\delta^{2}\right)}{(\Omega_{-}^{2}+\delta^{2})(\Omega_{+}^{2}+\delta^{2})}\right)\,,$
(118) $\displaystyle\sigma_{yy}=\sigma_{\rm D}\left(1-\frac{\omega_{\rm
c}^{2}(\Omega^{2}+\delta^{2})}{(\Omega_{-}^{2}+\delta^{2})(\Omega_{+}^{2}+\delta^{2})}\right)\,.$
(119)

Note that $\sigma_{\rm D}=e^{2}\rho_{\rm{2D}}/(m\delta)$ is the Drude dc
conductivity and that in $\sigma_{xy}$ we have introduced the Landau level
filling factor $\nu=\rho_{\rm{2D}}h/(eB)$.

It is also interesting to consider the limit of vanishing classical magnetic
field $\omega_{\rm c}\to 0$, where the system is just the bare electron gas
coupled to a cavity. In this case, we can see from Eqs. (113), (114) and (115)
that only $G^{\rm r}_{J_{x},J_{x}}$ is non-zero. This stems from the fact that
all the bare current response functions vanish in this limit, Cf. Eqs (110),
(111) and (112). Despite this, $G^{\rm r}_{J_{x},J_{x}}$ does not vanish
because it contains a term that depends solely on the cavity, i.e. the cavity
acts as a current channel in the direction of the light-matter coupling
independent of whether the material is a conductor or not. This is evidenced
by Eq. (91) where the last term solely depends on the free photonic
propagator. With this, the optical conductivity reads

$\sigma_{xx}(\omega)=\sigma_{\rm
D}\frac{i\delta}{\omega_{+}}\left(1+\frac{\omega_{\rm
p}^{2}}{\omega_{+}^{2}-\tilde{\omega}^{2}}\right)\,.$ (120)

This is in agreement with the optical conductivity of the 2D free electron gas
computed in Ref. [15]. Once again, there, they solve the system exactly,
computing the energies and wavefunctions and subsequenctly the current
response functions.

##  V Spin models in a cavity

###  V.1 Dicke model

We apply now our LRT to the paradigmatic Dicke model [40], for which $H_{\rm
m}=\frac{\omega_{z}}{2}\sum_{j}\sigma_{j}^{z}$ and
$C_{x}=\sum_{j}\sigma_{j}^{x}$. Here the matter subsystem is just a collection
of independent emitters and the only interactions are those mediated by the
cavity. Accordingly

$H_{\rm{eff}}^{\rm{MF}}=\frac{\omega_{z}}{2}\sum_{j}\sigma_{j}^{z}+\frac{2\lambda^{2}(\zeta-1)}{\Omega}m_{x}\sum_{j}\sigma_{j}^{x}-\frac{\lambda^{2}(\zeta-1)}{\Omega}m_{x}^{2}\,$
(121)

with $m_{x}=N^{-1}\langle\sum_{j}\sigma_{j}^{x}\rangle$. Note that for
$\zeta=1$, i.e. in the presence of a $P^{2}$ term, the effective term vanishes
and $H_{\rm{eff}}^{\rm{MF}}=H_{\rm m}$. Otherwise, Hamiltonian (121) can be
solved variationaly with respect to $m_{x}$ (See App. B for details)
obtaining, at zero temperature,

$m_{x}=\begin{cases}0&{\rm
if}\quad\lambda^{2}(1-\zeta)\leq\frac{1}{4}\omega_{z}\Omega\,,\\\
\sqrt{1-\left(\frac{\Omega\omega_{z}}{4\lambda^{2}(1-\zeta)}\right)^{2}}&{\rm
if}\quad\lambda^{2}(1-\zeta)>\frac{1}{4}\omega_{z}\Omega\,.\end{cases}$ (122)

For $\zeta=0$ the model exhibits a phase transition between a paramagnetic
phase with $m_{x}=0$ and a ferromagnetic phase (also termed superradiant when
the emphasis is put on the cavity) with $m_{x}\neq 0$. The response function
reads

$\tilde{\chi}_{xx,0}(\omega)=-\frac{\omega_{z}^{2}}{\varepsilon^{2}}\frac{2\varepsilon}{\omega_{+}^{2}-\varepsilon^{2}}\,,$
(123)

with
$\varepsilon^{2}=\omega_{z}^{2}+\left(\frac{4\lambda^{2}(1-\zeta)}{\Omega}m_{x}\right)^{2}$.
It depends on the value of $m_{x}$ and thus on the phase.

The photonic propagator, $D(\omega)$, which exhibits poles at the resonant
frequencies of the hybrid system [14]. The exact polaritons can be obtained
with a two-oscillator solution of the Dicke model in the thermodynamic limit.
In the case of $\zeta=0$, which corresponds to the standard Dicke model, these
read [41]

$2\Omega_{\pm}^{2}=\begin{cases}\omega_{z}^{2}+\Omega^{2}\pm\sqrt{\left(\omega_{z}^{2}-\Omega^{2}\right)^{2}+16\lambda^{2}\omega_{z}\Omega}&{\rm
if}\quad\lambda^{2}\leq\frac{1}{4}\omega_{z}\Omega\,,\\\
\omega_{z}^{2}/\mu^{2}+\Omega^{2}\pm\sqrt{\left(\omega_{z}^{2}/\mu^{2}-\Omega^{2}\right)^{2}+4\omega_{z}^{2}\Omega^{2}}&{\rm
if}\quad\lambda^{2}>\frac{1}{4}\omega_{z}\Omega\,,\end{cases}$ (124)

with $\mu=\omega_{z}\Omega/(4\lambda^{2})$. In the case of $\zeta=1$, the same
analysis is possible, except in this case there is only a paramagnetic phase
due to the presence of the $P^{2}$ term. In the two-oscillator solution, the
$P^{2}$ term affecting the spin becomes an $A^{2}$ term of the corresponding
oscillator. It can be eliminated with a Bogoliubov transform, yielding a
renormalized frequency
$\tilde{\omega}_{z}^{2}=\omega_{z}(\omega_{z}+4\lambda^{2}/\Omega)$ and
coupling
$\tilde{\lambda}=\lambda\left(1+4\lambda^{2}/(\omega_{z}\Omega)\right)^{-1/4}$.
Then, the polaritons are just those of the standard Dicke model in the normal
phase (top case of Eq. (124)), with the substitution
$\omega_{z}\to\tilde{\omega}_{z}$ and $\lambda\to\tilde{\lambda}$. For both
$\zeta=0$ and $\zeta=1$ we find perfect agreement between the exact Dicke
polaritons and the poles of the propagator computed with our linear response
theory [14].

###  V.2 Lipskin-Meshkov-Glick model with longitudinal coupling to the cavity

We consider now the Lipskin-Meshkov-Glick (LMG) model with longitudinal
coupling to the cavity as a marginal generalization of the Dicke model. In
contrast with the Dicke model, this one incorporates intrinsic all-to-all
interactions in the matter subsystem.

The full Hamiltonian reads

$\begin{split}H=&\frac{\omega_{z}}{2}\sum_{j}\sigma_{j}^{z}-\frac{J}{N}\sum_{ij}\sigma_{i}^{x}\sigma_{j}^{x}\\\
&+\Omega a^{\dagger}a+g(a+a^{\dagger})\sum_{j}\sigma_{j}^{x}\,.\end{split}$
(125)

Interestingly, the collective nature of the intrinsic interaction allows two
alternative treatments within the linear response theory. The first option is
to treat it like we would any other intrinsic interaction: bundle it together
with the field term to form the matter Hamiltonian, $H_{\rm
m}=\frac{\omega_{z}}{2}\sum_{j}\sigma_{j}^{z}-\frac{J}{N}\sum_{ij}\sigma_{i}^{x}\sigma_{j}^{x}$,
and then follow the linear response theory as described in Sec. III. The
second option exploits the collective nature of the intrinsic interaction and
treats this term not within $H_{\rm m}$ but explicitly within the linear
response theory, analogously to the $P^{2}$ term, as it is, in fact, just a
negative $P^{2}$ term [Cf. Eqs. (1) and (125)]. We will follow this second
option for its simplicity, since, as we show below, the resulting
$H_{\rm{eff}}^{\rm{MF}}$ corresponds to free spins and is identical in
structure to the one from the Dicke model (121). Accordingly, we have

$V_{\rm
ind}(\omega)=\frac{2\lambda^{2}}{\Omega}\frac{\Omega^{2}}{\omega^{2}-\Omega^{2}}-2J\,,$
(126)

and

$\begin{split}H_{\rm{eff}}^{\rm{MF}}=&\frac{\omega_{z}}{2}\sum_{j}\sigma_{j}^{z}-2J_{\rm
eff}m_{x}\sum_{j}\sigma_{j}^{x}\\\ &+J_{\rm eff}m_{x}^{2}\,,\end{split}$ (127)

with $m_{x}=N^{-1}\langle\sum_{j}\sigma_{j}^{x}\rangle$ and $J_{\rm
eff}=\frac{\lambda^{2}}{\Omega}+J$. Solving the Hamiltonian variationally with
respect to $m_{x}$ (See App. B for details) yields, at zero temperature,

$m_{x}=\begin{cases}0&{\rm if}\quad\omega_{z}\geq 4J_{\rm eff}\,,\\\
\sqrt{1-\left(\frac{\omega_{z}}{4J_{\rm eff}}\right)^{2}}&{\rm
if}\quad\omega_{z}<4J_{\rm eff}\,.\end{cases}$ (128)

Like the Dicke model (for $\zeta=0$), the LMG model exhibits a phase
transition between a paramagnetic phase with $m_{x}=0$ and a ferromagnetic
phase with $m_{x}\neq 0$. The response function reads [26]

$\tilde{\chi}_{xx,0}(\omega)=-\frac{\omega_{z}^{2}}{\varepsilon^{2}}\frac{2\varepsilon}{\omega_{+}^{2}-\varepsilon^{2}}\,,$
(129)

with $\varepsilon^{2}=\omega_{z}^{2}+\left(4J_{\rm eff}m_{x}\right)^{2}$. It
depends on the value of $m_{x}$ and thus on the phase.

Figure 1: Cavity response, $D$, of the LMG model with longitudinal coupling to
the cavity as a function of the collective coupling, $\lambda$, for different
values of the intrinsic interaction, $J$. The dashed lines correspond to a fit
of the polaritons with a two-oscillator model. The top right insets show the
magnetization. The transverse field is set to $\omega_{z}=\Omega$.

In Fig. 1 we show the photonic propagator, $D(\omega)$, which exhibits poles
at the resonant frequencies of the hybrid system, as a function of the
collective coupling. For $J<\omega_{z}/4$ the lower mode goes to zero at
finite coupling, signaling the phase transition. The critical point is
displaced to a lower coupling, $\lambda$, with respect to the Dicke model due
to the synergy between the intrinsic and effective interactions. For
$J>\omega_{z}/4$ the system is always on the ordered phase independently of
the value of $\lambda$. In both cases there is a splitting at zero coupling,
$\lambda=0$, due to the intrinsic interactions. In Fig. 1 we also overlay the
exact polaritons obtained with a two-oscillator solution of the LMG model in
the thermodynamic limit. It is is a generalization of the two-oscillator
solution of the Dicke model, developed in Ref. [41], to account for the
interaction term. The resulting LMG polaritons are

$2\Omega_{\pm}^{2}=\begin{cases}\tilde{\omega}_{z}^{2}+\Omega^{2}\pm\sqrt{\left(\tilde{\omega}_{z}^{2}-\Omega^{2}\right)^{2}+16\tilde{\lambda}^{2}\omega_{z}\Omega}&{\rm
if}\quad\omega_{z}\geq 4J_{\rm eff}\,,\\\
\omega_{z}^{2}/\tilde{\mu}^{2}-4\tilde{\mu}J\omega_{z}+\Omega^{2}\pm\sqrt{\left(\omega_{z}^{2}/\tilde{\mu}^{2}-4\tilde{\mu}J\omega_{z}+\Omega^{2}\right)^{2}+16\lambda^{2}\tilde{\mu}\omega_{z}\Omega}&{\rm
if}\quad\omega_{z}<4J_{\rm eff}\,,\end{cases}$ (130)

with $\tilde{\omega}_{z}^{2}=\omega_{z}(\omega_{z}-4J)$,
$\tilde{\lambda}=\lambda(1-4J/\omega_{z})^{-1/4}$ and
$\tilde{\mu}=\omega_{z}/(4J_{\rm eff})$. We find perfect agreement between the
exact LMG polaritons and the poles of the propagator computed with our linear
response theory.

This model describes recent experiments in rare earth quantum Ising systems,
such as LiHoF4 [21, 22]. These experiments aim to observe how the lowest
excitation mode softens as it approaches the ferromagnetic transition. We can
understand that Eq. (125) is equivalent to a Dicke model where $\lambda$
becomes $\sqrt{\lambda^{2}+J\Omega}$. Therefore, the softening of the mode can
be viewed as the behavior of the lower polariton in the Dicke model.

###  V.3 Lipskin-Meshkov-Glick model with transverse coupling to the cavity
(Dicke-LMG model)

Thus far we have only considered models where the intrinsic and cavity-
mediated interactions are synergistic, here we consider the Dicke-LMG model:
an LMG model with transverse coupling to the cavity, named in analogy to the
Dicke-Ising model that we will study in the next section. The full Hamiltonian
reads

$\begin{split}H=&\frac{\omega_{x}}{2}\sum_{j}\sigma_{j}^{x}+\frac{\omega_{z}}{2}\sum_{j}\sigma_{j}^{z}-\frac{J}{N}\sum_{ij}\sigma_{i}^{z}\sigma_{j}^{z}\\\
&+\Omega
a^{\dagger}a+g\left(a+a^{\dagger}\right)\sum_{j}\sigma_{j}^{x}\,.\end{split}$
(131)

For vanishing longitudinal and transverse fields the Dicke-LMG model has a
$\mathbb{Z}_{2}\times\mathbb{Z}_{2}$ symmetry. The first symmetry corresponds
to a spin flip, $\sigma_{j}^{z}\to-\sigma_{j}^{z}$, and in the bare LMG model
it is spontaneously broken in a second order phase transition from a
paramagnetic to a ferromagnetic phase. The second symmetry corresponds to a
simultaneous cavity-field and spin flip, $a\to-a$ and
$\sigma_{j}^{x}\to-\sigma_{j}^{x}$, and in the bare Dicke model it is
spontaneously broken in a second order phase transition from a normal to a
superradiant phase. As we show in the following, we find that the combination
of the two symmetries gives rise to a first-order phase transition in the
Dicke-LMG model between two symmetry-broken phases: a ferromagnetic normal
phase for large $J$ and a paramagnetic superradiant phase for large
$\lambda^{2}/\Omega$, this is akin to what occurs in the Dicke-Ising model
(See Sec. V.4). Switching on the external field has the effect of explicitely
breaking the corresponding symmetry, preventing the spontaneous symmetry
breaking that characterizes a phase transition. We find that if only one
external field is switched on, the phase transition is demoted from first to
second order. If the two external fields are switched on, the phase transition
is eliminated completely.

In order to solve the model, we will consider the intrinsic collective
interactions explicitly within the linear response theory. However, unlike in
the case of the LMG model with longitudinal coupling to the cavity, in this
case the cavity-mediated and intrinsic interactions act on different axes,
i.e. they couple different collective operators:
$C_{x}=\sum_{j}\sigma_{j}^{x}$ and $C_{z}=\sum_{j}\sigma_{j}^{z}$
respectively. Accordingly, a slight generalization of the linear response
theory presented in Sec. III is required as there are now two induced
interaction terms within the effective action

$\displaystyle
V_{\rm{ind},x}(\omega)=\frac{2\lambda^{2}}{\Omega}\frac{\Omega^{2}}{\omega_{+}^{2}-\Omega^{2}}\,,$
(132) $\displaystyle V_{\rm{ind},z}(\omega)=-2J\,.$ (133)

Thus, the system can be solved using a multimode generalization of the linear
response theory, developed in App. A. The expression for the photonic
propagator remains unchanged, see Eq. 30, as it solely determined by the
operator that couples to the cavity. In contrast, matter correlators have new
contributions afforded by the new interaction channel, $V_{\rm{ind},z}$. We
find

$\chi_{xx}=\frac{\tilde{\chi}_{xx,0}+V_{\rm{ind},z}\det(\tilde{\chi}_{0})}{1+V_{\rm{ind},x}\tilde{\chi}_{xx,0}+V_{\rm{ind},z}\tilde{\chi}_{zz,0}+V_{\rm{ind},x}V_{\rm{ind},z}\det(\tilde{\chi}_{0})}\,,$
(134)

with $(\tilde{\chi}_{0})_{ab}=\tilde{\chi}_{ab,0}$. The mean-field effective
Hamiltonian reads

$H_{\rm{eff}}^{\rm{MF}}=\frac{\tilde{\omega}_{x}}{2}\sum_{j}\sigma_{j}^{x}+\frac{\tilde{\omega}_{z}}{2}\sum_{j}\sigma_{j}^{z}+\frac{N\lambda^{2}}{\Omega}m_{x}^{2}+NJm_{z}^{2}\,,$
(135)

with $\tilde{\omega}_{x}=\omega_{x}-4\frac{\lambda^{2}}{\Omega}m_{x}$ and
$\tilde{\omega}_{z}=\omega_{z}-4Jm_{z}$. Solving variationaly with respect to
$m_{x}$ and $m_{z}$ (See App. B for details) allows us to compute the
equilibrium values of $m_{x}$ and $m_{z}$ numerically and subsequently the
response functions $\tilde{\chi}_{xx,0}$, $\tilde{\chi}_{xz,0}$,
$\tilde{\chi}_{zx,0}$ and $\tilde{\chi}_{zz,0}$, which depend on the values of
$m_{x}$ and $m_{z}$.

In Fig. 2 we show the photonic propagator, $D(\omega)$, the matter response
function, $\chi_{zz}(\omega)$, (left inset) and the magnetizations, $m_{x}$
and $m_{z}$, (right inset) as functions of the collective coupling, $\lambda$,
in four different scenarios of external fields values. The intrinsic
interaction is set so that at zero coupling the cavity and the spins are
resonant: $4J=\Omega$. In the case of vanishing $\omega_{x}$ and $\omega_{z}$,
Fig. 2(a), we observe a first order phase transition, evidenced by the
discontinuous behaviour in the order parameters of the two ordered phases. The
photonic propagator has a pole that goes to zero at the critical point,
signaling the gap closing. We include $\chi_{zz}$ in the left inset for
completeness because the photonic propagator depends on $\chi_{xx}$ [Cf. Eq.
(30)] and the system becomes unresponsive to probing along the $x$ direction
in the superradiant phase. This is caused by the fact that the intrinsic
interaction is the only term in the Hamiltonian containing $\sigma^{z}$
operators. Due to its collective nature, its effect is switched off (mean-
field behaviour) in the superradiant phase where $m_{z}=0$ [Cf. Eq. (135)].
The system is then fully ordered along $x$ and with no off-diagonal terms it
becomes unresponsive to probing with $\sigma^{x}$. Accordingly, in the
superradiant phase the photonic propagator only shows a pole at the cavity
frequency. To witness the other excitation of the system we look at
$\chi_{zz}$ which shows a pole that emerges from zero at the critical point,
signaling the gap reopening. The combination of $D$ and $\chi_{zz}$ provides a
complete picture of the excitations of the system. In Fig. 2(b) we see that
switching on the transverse field, $\omega_{x}$, eliminates the superradiant
symmetry. The phase transition becomes of second order, with $m_{z}$ the sole
order parameter. Like in the previous case, the absence of non mean-field
terms containing $\sigma^{z}$ in the Hamiltonian makes the system unresponsive
to $\sigma_{x}$ in the superradiant phase. Thus, we also include $\chi_{zz}$
to witness the gap reopening. Alternatively, in Fig. 2(c) we see that
switching on the longitudinal field, $\omega_{z}$, eliminates the
ferromagnetic symmetry. The phase transition becomes of second order, with
$m_{x}$ as the sole order parameter. In this case, the presence of the
longitudinal field endows the photonic propagator with visibility of all the
excitations, showcasing the full gap closing and reopening at the critical
point. Finally in Fig. 2(d) we switch on both the longitudinal and transverse
fields. All symmetries are now explicitely broken and the phase transition
gives way to a smooth crossover of $m_{z}$ and $m_{x}$ dominant regimes for
small and large $\lambda$ respectively. The gap remains finite at all times.

Figure 2: Cavity response, $D$, of the Dicke-LMG model (LMG model with
transverse coupling to the cavity) as a function of the collective coupling,
$\lambda$, for different values of the longitudinal and transverse fields,
$\omega_{z}$ and $\omega_{x}$ respectively. The top left insets show the
matter response $\chi_{zz}$. The top right insets show the magnetizations. The
intrinsic interaction is set to $4J=\Omega$.

###  V.4 Ising model with transverse coupling to the cavity (Dicke-Ising
model)

Thus far we have only considered models with collective intrinsic
interactions, such that despite being intrinsic, these collective interactions
were handled on a par with the light-matter interactions using the linear
response theory. As a result, the mean-field effective Hamiltonian
corresponded to free spins and was trivially solved. Here we consider the
Dicke-Ising model: a Ising model with transverse coupling to the cavity [30,
42, 43, 44]. The intrinsic interactions are between nearest-neighbours, i.e.
not collective, and as such they cannot be handled within the linear response
theory like a $P^{2}$ term. The full Hamiltonian reads

$\begin{split}H=&\frac{\omega_{x}}{2}\sum_{j}\sigma_{j}^{x}-J\sum_{j}\sigma_{j}^{z}\sigma_{j+1}^{z}\\\
&+\Omega
a^{\dagger}a+g\left(a+a^{\dagger}\right)\sum_{j}\sigma_{j}^{x}\,.\end{split}$
(136)

Like the LMG-Dicke model, and for the same reasons (See Sec. V.2), for
vanishing transverse field the Dicke-Ising model has a
$\mathbb{Z}_{2}\times\mathbb{Z}_{2}$ symmetry that gives rise to a first-order
phase transition between two symmetry-broken phases: a ferromagnetic normal
phase for large $J$ and a paramagnetic superradiant phase for large
$\lambda^{2}/\Omega$. However, unlike in the Dicke-LMG model, the phase
transition remains of first order after switching on the transverse field. We
attribute this robustness to the fact that the intrinsic interactions are not
collective in this case. As we will see below, we do not consider a
longitudinal field because the resulting mean-field effective Hamiltonian
would not be analytically solvable.

The corresponding mean-field effective Hamiltonian is that of a Ising chain in
transverse field

$H_{\rm{eff}}^{\rm{MF}}=\frac{\tilde{\omega}_{x}}{2}\sum_{j}\sigma_{j}^{x}-J\sum_{j}\sigma_{j}^{z}\sigma_{j+1}^{z}+\frac{N\lambda^{2}}{\Omega}m_{x}^{2}\,,$
(137)

with $\tilde{\omega}_{x}=\omega_{x}-4\lambda^{2}/\Omega m_{x}$ and
$m_{x}=N^{-1}\sum_{j}\langle\sigma_{j}^{x}\rangle$. The transverse field is a
combination of the external field and the mean-field cavity-induced
interaction. It is clear now that adding a longitudinal field would make
$H_{\rm{eff}}^{\rm{MF}}$ analytically intractable. In the thermodynamic limit,
$N\to\infty$, the ground-state energy per spin is given by

$e_{0}(m_{x})=\frac{\lambda^{2}}{\Omega}m_{x}^{2}-\frac{1}{2}\int_{-\pi}^{\pi}\frac{dk}{2\pi}\epsilon_{k}\,,$
(138)

with

$\epsilon_{k}=\sqrt{(2J)^{2}+\tilde{\omega}_{x}^{2}-4J\tilde{\omega}_{x}\cos
k}\,.$ (139)

Solving variationaly allows us to compute the equilibrium value of $m_{x}$
numerically and subsequently the response function $\tilde{\chi}_{xx,0}$,
which depends on the value of $m_{x}$. At zero temperature and in the
continuum limit, $\tilde{\chi}_{xx,0}$ is given by (See App. C for details)

$\tilde{\chi}_{xx,0}(\omega)=-32J^{2}\int_{-\pi}^{\pi}\frac{dk}{2\pi}\frac{\sin^{2}k}{\epsilon_{k}(\omega^{2}-4\epsilon_{k}^{2})}\,.$
(140)

Interestingly, we find that $\tilde{\chi}_{xx,0}$ has poles at
$2\epsilon_{k}$. This stems from the fact that

$\begin{split}C_{x}=\sum_{j}\sigma_{j}^{x}=N-2\sum_{k}\bigl{(}&v_{k}^{2}+(u_{k}^{2}-v_{k}^{2})\gamma_{k}^{\dagger}\gamma_{k}\\\
&+iu_{k}v_{k}(\gamma_{k}^{\dagger}\gamma_{-k}^{\dagger}-\gamma_{-k}\gamma_{k})\bigr{)}\,,\end{split}$
(141)

where $\gamma_{k}$ and $\gamma_{k}^{\dagger}$ are the annihilation and
creation operators of the Bogoliuvov fermions that constitute the elementary
excitations of the Ising model after a Jordan-Wigner fermionization, with
$u_{k}$ and $v_{k}$ the Bogoliubov coefficients [45][Chap. 10]. The physical
consequences of this coupling where discussed in Ref. [14]. Here we present a
comparison between our LRT and finite-size exact-diagonalization results in
Fig. 3. This allows us to discuss how quickly finite-size effects are washed
out as we increase the system size and provides another check of our formulas
in this non-trivial case.

Figure 3: Cavity response, $D$, of the Dicke-Ising model as a function of the
collective coupling, $\lambda$, computed analytically in the thermodynamic
limit, $N\to\infty$, (left) and with exact diagonalization, $N=14$, (right)
for different values of the classical field, $\omega_{x}$. The yellow dashed
lines correspond to a fit of the polaritons with a two-oscillator model. The
top right insets show the magnetizations. The bottom right inset in the right
plots shows a vertical cut at the black dashed line for two finite sizes,
$N=4$ and $N=14$. The dotted lines mark the edges of the band of the mean-
field effective Hamiltonian (137). The parameters are $\omega_{x}=0$ and
$4J=\Omega$. In the exact-diagonalization results, the Fock basis for the
photonic Hilbert space is truncated at 40 photons.

We compare our LRT in Figs. 3(a), (c) and (e), valid in the thermodynamic
limit, with exact-diagonalization results for system sizes up to $N=14$ in
Figs. 3(b), (d) and (f). Let us begin by comparing Figs. 3 (a) and (b), which
correspond to the case of vanishing classical field. In this case we observe
the formation of polaritons in the normal phase and the opening of the two-
excitation band and the hardening of the lower polartion in the superradiant
phase in Fig. 3(a). The same features are observed in Fig. 3(b), although the
two-excitation band is not fully formed and instead we can distinguish a
collection of discrete levels. Additionally, there are some finite-size
effects. Most prominently, there is a pole corresponding to the single-
excitation band in the normal phase, at $\omega=2J=\Omega/2$ for $\lambda\to
0$. This is explained by noting that the coupling operator $C_{x}$ can create
domain walls at the edges of the chain, something which is impossible in the
thermodynamic limit, an effect that becomes negligible as $N\to\infty$. This
interpretation is confirmed by the bottom right inset in Fig. 3(b), which
shows that the intensity of this pole decreases with size, unlike the poles
corresponding to the polaritons. Figures 3(c), (d), (e) and (f) feature a
finite classical field $\omega_{x}$ and thus a finite bandwidth in the normal
phase. In Fig. 3(c) the narrow bandwidth allows the formation of bound
polariton states (BPS) with well defined energies below and above the band
[14]. This is validated in Fig. 3(d). Again, we observe additional features
that we attribute to finite size-effects. In particular the single-excitation
band, at $\omega=2J=\Omega/2$ for $\lambda\to 0$. Its visibility is shown to
decrease with size in the bottom right inset of Fig. 3(d). In Fig. 3(e) the
large bandwidth prevents the formation for BPS. This is confirmed in Fig. 3(f)
where we observe a collection of closely packed levels of equal visibility
that are expected to form the band in the thermodynamic limit. There are not
any two levels that stand out as polaritons, in agreement with our LRT. The
levels that fall outside the would-be band are shown to be finite-size
artifacts in the bottom right inset of Fig. 3(f).

###  V.5 Heisenberg ferromagnet coupled to a cavity

We finish the section by commenting on the paradigmatic Heisenberg model. In
particular we consider a ferromagnetic Heisenberg model on a $3D$ lattice,
such that the model is ordered below a critical temperature,

$H_{\rm m}=-J\sum_{\langle
ij\rangle}\boldsymbol{\sigma}_{i}\boldsymbol{\sigma}_{j}+\frac{\omega_{z}}{2}\sum_{j}\sigma_{j}^{z}\,.$
(142)

The ground state of the model corresponds to a symmetry broken ferromagnetic
state. For vanishing field, $\omega_{z}=0$, the spins will be magnetized in a
random direction of space, spontaneously breaking the $SO(3)$ symmetry of the
model. The addition of a field demotes the symmetry to $SO(2)$ in the plane
perpendicular to the field, setting the magnetization direction along the
field. It also opens a gap between the spin-wave band and the ground state,
such that the zero momentum spin-wave has an energy gap with respect to the
ground state of $\omega_{z}$. The coupling to a cavity in a direction
perpendicular to the classical field will demote the remaining $SO(2)$
symmetry to $\mathbb{Z}_{2}$ in the cavity field direction. The cavity
coupling, $C_{x}=\sum_{j}\sigma_{j}^{x}$, will compete with the classical
field in setting the magnetization direction. The corresponding effective
Hamiltonian reads

$\begin{split}H_{\rm{eff}}^{\rm{MF}}=&\frac{\omega_{z}}{2}\sum_{j}\sigma_{j}^{z}-J\sum_{\langle
ij\rangle}\boldsymbol{\sigma}_{i}\boldsymbol{\sigma}_{j}\\\
&-\frac{2\lambda^{2}}{\Omega}m_{x}\sum_{j}\sigma_{j}^{x}+\frac{\lambda^{2}}{\Omega}m_{x}^{2}\,\end{split}$
(143)

If we assume that the model is ordered, the interaction term contributes an
energy $-Jz$ with $z$ the coordination number, independently of the value of
$m_{x}$. Thus, the variational Hamiltonian is equivalent to that of a Dicke
model [Cf. Sec. V.1.]. The behaviour of the magnetization, $m_{x}$, is also
inherited from the Dicke model. For a subcritical coupling,
$\lambda<\sqrt{\omega_{z}\Omega}/2$, the model is ordered along $z$. When the
coupling reaches the critical point, the magnetization starts to turn towards
the $x$ axis in a second order phase transition.

The response function $\tilde{\chi}_{xx,0}$ corresponds to the excitation of
the zero-momentum magnon, with a pole at
$\omega^{2}=\omega_{z}^{2}+(4\lambda^{2}m_{x}/\Omega)^{2}$. Therefore, our LRT
predicts the formation of polaritons of this zero-momentum magnon with
frequencies given by Eq. 124, following the physics of the Dicke model.

##  VI Conclusions

In this paper, we have presented a linear response theory for materials
collectively coupled to a cavity. We have employed two different approaches.
The first relies on a path integral formulation of the partition function and
a saddle point expansion of the partition function that is truncated exactly
in the thermodynamic limit. The second is more direct and relies on
formulating the equations of motion for the response functions. These are
truncated to second order in light-matter fluctuations and solved, yielding
equivalent expressions for the response functions to the first approach. This
provides a posteriori validation for the truncation of the equations of
motion, which confirms the validity of a mean-field decoupling of the light-
matter interaction in cavity QED materials. We obtain exact expressions for
the response functions of the hybrid system in terms of the bare response
functions of the cavity and the material.

The theory has been demonstrated in several systems, showcasing its
applicability. In the Quantum Hall effect we have recovered the optical
response and Landau polaritons computed previously in Refs. [33, 34]. In
magnetic systems, we have started from the non-interacting Dicke model and
progressively complicated the model by adding collective intrinsic
interactions to yield the LMG model with longitudinal coupling to the cavity,
then changing the cavity coupling direction to induce a competition between
intrinsic and light-matter interactions in the Dicke-LMG model and finally
making the intrinsic interactions non-collective in the Dicke-Ising model.
Finally, we have explored the effect of the cavity in an ordered Heisenberg
ferromagnet. Our predictions have been validated against analytical and exact-
diagonalization results.

Even if the focus of the paper has been the linear response theory, the study
of the different spin models has provided some interesting independent
physical insights. In particular that first order phase transitions arise in
hybrid systems with competing orders, such as the Dicke-LMG and Dicke-Ising
models, when the competing constituents feature second-order phase transitions
by themselves. In this sense, our results also point to the fact that the
robustness of this first order phase transition depends on the nature of the
intrinsic material interactions. A transverse field demotes the phase
transition to second order for the Dicke-LMG model (collective interactions)
but not for for the Dicke-Ising model (nearest-neighbor interactions). We
underline that this has been witnessed through the lens of the photonic
response, which is relevant for its experimental implications. Cavity
transmission is a common probing mechanism for magnetic materials inductively
coupled to a cavity [46, 47, 48, 49, 50]. Finally, the finite-size exact-
diagonalization results allow us to gauge what the commonly-used
“thermodynamic limit” means in terms of actual system sizes. As shown in Sec.
V.4, all the features (bar the blending of discrete levels to form the band)
expected in the thermodynamic limit are already present for $N=14$.

Although not explored here, the computation of finite-size corrections is
possible. In the path-integral approach this would be achieved by considering
higher order terms in the saddle-point expansion. Alternatively, one could
truncate the equations of motion at higher order, to include three- or higher-
order correlators [51, 52, 53]. This would be relevant in the study of
topological systems, where boundary effects are crucial [54, 55, 56]. Also,
real time dynamics could be easily studied with the framework of the equations
of motion [57]. In conclusion, our work lays the foundations of a linear
response theory for cavity QED materials.

## Acknowledgements

We acknowledge discussions with Katharina Lenk, Martin Eckstein, Andreas Wipf,
Yuto Ashida, Luis Martín-Moreno, Fernando Falceto and José G. Esteve. Despite
the enormous impediments in spending the budget, the authors must acknowledge
funding from the grant TED2021-131447B-C21 funded by
MCIN/AEI/10.13039/501100011033 and the EU ‘NextGenerationEU’/PRTR. We also
acknoledge the Gobierno de Aragón (Grant E09-17R Q-MAD), Quantum Spain and the
CSIC Quantum Technologies Platform PTI-001. J. R-R acknowledges support from
the Ministry of Universities of the Spanish Government through the grant
FPU2020-07231.

##  Appendix A Linear response theory for a multimode cavity

Consider the Hamiltonian for a material coupled collectively to a multimode
cavity

$\begin{split}H=&H_{\rm m}+\sum_{k}\Omega_{k}a_{k}^{\dagger}a_{k}\\\
&+\sum_{k}g_{k}\left(a_{k}+a_{k}^{\dagger}\right)C_{k}+\zeta\sum_{k}\frac{g_{k}^{2}}{\Omega_{k}}C_{k}^{2}\,.\end{split}$
(144)

The corresponding action reads

$\begin{split}S=&S_{\rm
m}+\sum_{k}\int_{\tau}\bar{a}_{k}(\tau)(\partial_{\tau}-\Omega_{k})a_{k}(\tau)\\\
&+\sum_{k}g_{k}\int_{\tau}(a_{k}(\tau)+\bar{a}_{k}(\tau))C_{k}(\tau)\\\
&+\zeta\sum_{k}\frac{g_{k}^{2}}{\Omega_{k}}\int_{\tau}C_{k}^{2}(\tau)\,.\end{split}$
(145)

After partial integration over the cavity field, the induced action reads

$S_{\rm
ind}=\frac{1}{2}\sum_{k}\int_{\tau,\tau^{\prime}}C_{k}(\tau)\frac{1}{N}V_{\rm
ind,k}(\tau-\tau^{\prime})C_{k}(\tau^{\prime})\,,$ (146)

with

$V_{\rm
ind,k}(\omega_{m})=2\lambda_{k}^{2}\frac{\Omega_{k}^{2}(\zeta-1)+\zeta\omega_{m}^{2}}{\Omega_{k}(\omega_{m}^{2}+\Omega_{k}^{2})}\,.$
(147)

Here $D_{k,0}$ is the free photon propagator of the $k$-th mode

$D_{k,0}(\omega_{m})=\frac{1}{i\omega_{m}-\Omega_{k}}\,.$ (148)

The induced interaction can be decoupled with a Hubbard-Stratonovich
transformation that introduces an auxiliary scalar field for each mode,
$\varphi_{k}$:

$e^{-S_{\mathrm{ind}}}=\frac{1}{Z_{\varphi}}\oint_{\varphi}e^{-\frac{1}{2}N\sum_{k}\int_{\tau,\tau^{\prime}}\varphi_{k}(\tau)V_{\mathrm{ind},k}^{-1}(\tau-\tau^{\prime})\varphi_{k}(\tau^{\prime})-i\sum_{k}\int_{\tau}\varphi_{k}(\tau)C_{k}(\tau)}\,.$
(149)

We define the propagator of the auxiliary fields as

$W_{kp}(\tau)=\langle\varphi_{k}(\tau)\varphi_{p}(0)\rangle^{\rm c}\,.$ (150)

The generating functional for bare connected correlation functions of the
matter coupling operators, $C_{k}$, is $\mathcal{G}_{\rm
m}^{0}[\xi_{k}]=-N^{-1}\log Z_{\rm m}[\xi_{k}]$ with

$Z_{\rm m}[\xi_{k}]=\oint_{c}e^{-(S_{\rm
m}+i\sum_{k}\int_{\tau}\xi_{k}(\tau)C_{k}(\tau))}\,,$ (151)

such that

$\begin{split}\frac{\delta}{\delta\xi_{k_{1}}(\tau_{1})}\cdots&\frac{\delta}{\delta\xi_{k_{n}}(\tau_{n})}\mathcal{G}_{\rm
m}^{0}[\xi_{k}]\Bigr{|}_{\xi_{k}=0}=\\\ &=-\frac{(-i)^{n}}{N}\langle
C_{k_{1}}(\tau_{1})\cdots C_{k_{n}}(\tau_{n})\rangle^{\rm c}_{\rm m}\\\
&\equiv\chi_{k_{1}\ldots
k_{n},0}^{(n)}(\tau_{1},\ldots,\tau_{n})\,.\end{split}$ (152)

In the following we will denote $\chi_{kp,0}^{(2)}\equiv\chi_{kp,0}$. With
this, after partial integration over the matter degrees of freedom, we can
write

$Z=\oint_{\varphi_{k}}e^{-Nf[\varphi_{k}]}\,,$ (153)

with

$f[\varphi_{k}]=\frac{1}{2}\sum_{k}\int_{\tau,\tau^{\prime}}\varphi_{k}(\tau)V_{\mathrm{ind},k}^{-1}(\tau-\tau^{\prime})\varphi_{k}(\tau^{\prime})+\mathcal{G}_{\rm
m}^{0}[\varphi_{k}]\,.$ (154)

Like in the single-mode case, we apply a saddle-point approximation for large
$N$. This requires that the number of modes, $M$, be finite, i.e. that
$\lim_{N\to\infty}M/N\to 0$. Otherwise higher order terms in the saddle-point
expansion might not actually be negligible. The condition that the functional
derivatives of $f$ with respect to $\varphi_{k}$ vanish

$\varphi_{k,\rm sp}(\omega_{m})=-\frac{i}{N}V_{\rm ind,k}(\omega_{m})\langle
C_{k}(\omega_{m})\rangle_{\varphi_{k,\rm sp}}$ (155)

formally defines $\varphi_{k,\rm sp}$. Making the simplifying assumption that
$\varphi_{k,\rm sp}$ be constant we find

$\varphi_{k,\rm sp}=-\frac{i}{N}V_{\rm{ind},k}(\omega_{m}=0)\langle
C_{k}\rangle_{\varphi_{k,\rm sp}}\,.$ (156)

This tells us that $\varphi_{k,\rm sp}$ is self-consistently proportional to
$\langle C_{x}\rangle_{\varphi_{k,\rm sp}}$, i.e. to the expectation value of
$C_{k}$ for the bare matter subject to fields
$i\varphi_{k,\rm{sp}}=\frac{1}{N}V_{\rm ind,k}(\omega_{m}=0)\langle
C_{k}\rangle$. This is precisely the self-consistent condition that arises
from computing $\langle C_{k}\rangle$ from the mean-field Hamiltonian

$\begin{split}H_{\rm eff}^{\rm MF}=&H_{\rm
m}+\sum_{k}\frac{2\lambda_{k}^{2}(\zeta-1)}{N\Omega_{k}}\langle C_{k}\rangle
C_{k}\\\ &-\sum_{k}\frac{\lambda_{k}^{2}(\zeta-1)}{N\Omega_{k}}\langle
C_{k}\rangle^{2}\,.\end{split}$ (157)

This is the mean-field theory of the effective Hamiltonian that arises from
taking the static limit in $V_{\rm ind,k}$ in the effective action [8].

Expanding $f$ around the saddle-point yields

$(NW_{0})_{kp}^{-1}=\delta_{kp}V_{\rm{ind},k}^{-1}+\tilde{\chi}_{kp,0}\,,$
(158)

with

$\begin{split}\tilde{\chi}_{kp,0}(\tau_{1},\tau_{2})=&\frac{\delta}{\delta\varphi_{k,\rm
sp}(\tau_{1})}\frac{\delta}{\delta\varphi_{p,\rm
sp}(\tau_{2})}\mathcal{G}_{\rm m}^{0}\\\ =&\frac{1}{N}\langle
C_{k}(\tau)C_{p}(\tau^{\prime})\rangle_{\varphi_{k,\rm{sp}}}^{\rm
c}\,.\end{split}$ (159)

For $N\to\infty$ one can safely truncate the saddle-point expansion to second
order, which implies that $W=W_{0}$.

We can now obtain relations between the auxiliary field propagator, $W_{kp}$,
the photonic propagator, $D_{k}$, and matter response functions. Let us define
the generating functional for photonic connected correlators by introducing a
complex fields in Eq. (145): $\mathcal{G}_{\rm
ph}[\eta_{k},\bar{\eta}_{k}]=-\log Z[\eta_{k},\bar{\eta}_{k}]$, with

$Z[\eta_{k},\bar{\eta}_{k}]=\oint_{a_{k},\bar{a}_{k},c}e^{-(S+\sum_{k}\int_{\tau}(a_{k}(\tau)\bar{\eta}_{k}(\tau)+\bar{a}_{k}(\tau)\eta_{k}(\tau))}\,.$
(160)

After partial integration of the cavity fields, this yields

$\begin{split}Z[\eta_{k},\bar{\eta}_{k}]=\oint_{c}\exp&-\Bigl{(}S_{\rm
m}+\zeta\sum_{k}\frac{g_{k}^{2}}{\Omega_{k}}\int_{\tau}C_{k}^{2}(\tau)\\\
&+\sum_{k}\int_{\tau,\tau^{\prime}}\bar{m}_{k}(\tau)D_{k}(\tau-\tau^{\prime})m_{k}(\tau)\Bigr{)}\,,\end{split}$
(161)

with $m_{k}(\tau)=\eta_{k}(\tau)+g_{k}C_{k}(\tau)$. With this

$D_{k}(\omega_{m})=D_{k,0}(\omega_{m})-\lambda_{k}^{2}D_{k,0}(\omega_{m})\chi_{kk}(\omega_{m})D_{k,0}(\omega_{m})\,.$
(162)

Likewise, we can define the generating functional for matter connected
correlators: $\mathcal{G}_{\rm m}[\xi_{k}]=-N^{-1}\log Z[\xi_{k}]$ with

$Z[\xi_{k}]=\oint_{c,\varphi_{k}}e^{-(S_{\rm
eff}+i\sum_{k}\int_{\tau}\xi_{k}(\tau)C_{k}(\tau))}\,.$ (163)

After partial integration over the cavity fields, this yields

$Z[\xi_{k}]=\oint_{\varphi_{k}}e^{-Nf[\xi_{k}+\varphi_{k}]}$ (164)

Then, a second order expansion of $f[\xi_{k}+\varphi_{k}]$ around the saddle
point yields

$\begin{split}f[\xi_{k}+\varphi_{k}]-f[\varphi_{k,\rm{sp}}]=\frac{1}{2}\sum_{k,p}\int_{\tau,\tau^{\prime}}\Bigl{(}&\delta\varphi_{k}(\tau)(NW)_{kp}^{-1}(\tau-\tau^{\prime})\delta\varphi_{p}(\tau^{\prime})+\xi_{k}(\tau)\tilde{\chi}_{kp,0}(\tau-\tau^{\prime})\xi_{p}(\tau^{\prime})\\\
&+\,\delta\varphi_{k}(\tau)\tilde{\chi}_{kp,0}(\tau-\tau^{\prime})\xi_{p}(\tau^{\prime})+\xi_{k}(\tau)\tilde{\chi}_{kp,0}(\tau-\tau^{\prime})\delta\varphi_{p}(\tau)\Bigr{)}\,.\end{split}$
(165)

The functional integral over the auxiliary-field displacements
$\delta\varphi_{k}$ is just an $M$-dimensional Gaussian integral that we can
perform, yielding

$G_{\rm m}[\xi_{k}]={\rm
cst.}+\frac{1}{2}\sum_{k,p}\int_{\tau,\tau^{\prime}}\xi_{k}(\tau)\left(\tilde{\chi}_{kp,0}(\tau-\tau^{\prime})-\sum_{k^{\prime},p^{\prime}}\int_{u,v}\tilde{\chi}_{kk^{\prime},0}(\tau-u)NW_{k^{\prime}p^{\prime}}(u-v)\tilde{\chi}_{p^{\prime}p,0}(v-\tau^{\prime})\right)\xi_{p}(\tau^{\prime})\,.$
(166)

With this

$\chi_{kp}(\omega_{m})=\tilde{\chi}_{kp,0}(\omega_{m})-\sum_{k^{\prime},p^{\prime}}\tilde{\chi}_{kk^{\prime},0}(\omega_{m})NW_{k^{\prime}p^{\prime}}(\omega_{m})\tilde{\chi}_{p^{\prime}p,0}(\omega_{m})\,.$
(167)

Note that obtaining final expressions for $\chi_{kp}(\omega_{m})$ requires
inverting the matrix $(NW)^{-1}$ that is defined in Eq. (158).

##  Appendix B A spin subject to variational fields

###  B.1 Spectrum and response functions of a free spin

Let us consider the following free spin model

$H=\frac{\omega_{x}}{2}\sigma_{x}+\frac{\omega_{z}}{2}\sigma_{z}\,.$ (168)

It is diagonalized by a rotation that defines spin operators along new
directions

$\displaystyle\sigma_{z}^{\prime}=\frac{\omega_{z}}{\epsilon}\sigma_{z}-\frac{\omega_{x}}{\epsilon}\sigma_{x}\,,$
(169)
$\displaystyle\sigma_{x}^{\prime}=\frac{\omega_{x}}{\epsilon}\sigma_{z}+\frac{\omega_{z}}{\epsilon}\sigma_{x}\,,$
(170)

where $\epsilon^{2}=\omega_{x}^{2}+\omega_{z}^{2}$. With these, the
Hamiltonian reads

$H=\frac{\epsilon}{2}\sigma_{z}^{\prime}\,,$ (172)

with ground-state energy $E_{0}=-\epsilon/2$.

We compute the response functions from their spectral decomposition, which, at
zero temperature, reads

$\chi_{ab}(\omega)=-\sum_{n}\left(\frac{\langle 0|\sigma_{a}|n\rangle\langle
n|\sigma_{b}|0\rangle}{\omega^{+}+E_{0}-E_{n}}-\frac{\langle
n|\sigma_{a}|0\rangle\langle
0|\sigma_{b}|n\rangle}{\omega^{+}+E_{n}-E_{0}}\right)\,,$ (173)

where the sum is over the eigenstates of the system. With the matrix elements

$\displaystyle\langle n|\sigma_{x}|0\rangle=-\frac{\omega_{x}}{\epsilon}\,,$
(174) $\displaystyle\langle
n|\sigma_{z}|0\rangle=\frac{\omega_{z}}{\epsilon}\,,$ (175)

we obtain

$\displaystyle\chi_{xx}(\omega)=-\frac{\omega_{z}^{2}}{\epsilon^{2}}\frac{2\epsilon}{\omega_{+}^{2}-\epsilon^{2}}\,,$
(176)
$\displaystyle\chi_{xz}(\omega)=\chi_{zx}(\omega)=-\frac{\omega_{x}\omega_{z}}{\epsilon^{2}}\frac{2\epsilon}{\omega_{+}^{2}-\epsilon^{2}}\,,$
(177)
$\displaystyle\chi_{zz}(\omega)=-\frac{\omega_{x}^{2}}{\epsilon^{2}}\frac{2\epsilon}{\omega_{+}^{2}-\epsilon^{2}}\,.$
(178)

###  B.2 Considering variational fields

There will be situations in which we arrive to a free spin model like the one
defined in Eq. (168) from a mean field approximation. As a result we will have
a Hamiltonian of the form

$H(m_{x},m_{z})=\frac{\tilde{\omega}_{x}}{2}\sigma_{x}+\frac{\tilde{\omega}_{z}}{2}\sigma_{z}+E(m_{x},m_{z})\,,$
(179)

where now the fields, $\tilde{\omega}_{x}$ and $\tilde{\omega}_{z}$ depend on
mean-field parameters $m_{a\in\\{x,z\\}}=\langle\sigma_{a}\rangle$ and there
is a constant energy term that also depends on the variational parameters. The
Hamiltonian can be diagonalized as explained in the previous section, and the
ground state energy and response functions are now functions of $m_{x}$ and
$m_{z}$. In particular,
$E_{0}(m_{x},m_{z})=-\epsilon(m_{x},m_{z})/2+E_{0}(m_{x},m_{z})$, with
$\epsilon^{2}=\tilde{\omega}_{x}^{2}+\tilde{\omega}_{z}^{2}$. The values of
$m_{x}$ and $m_{z}$ are determined variationaly from a minimization of
$E_{0}(m_{x},m_{z})$. Depending on the particular problem, this can be done
analytically or numerically.

##  Appendix C Computing $\tilde{\chi}_{xx,0}$ for the Ising model

The spectral decomposition of $\tilde{\chi}_{xx,0}$ reads

$\tilde{\chi}_{xx,0}=-\frac{1}{N}\sum_{n}|\langle
n|C_{x}|0\rangle|^{2}\frac{2(E_{n}-E_{0})}{\omega_{+}^{2}-(E_{n}-E_{0})^{2}}\,,$
(180)

where $|n\rangle$ is the eigenstate of $H_{\rm{eff}}^{\rm{MF}}$ (137) with
eigenvalue $E_{n}$. From Eq. (141) we see that

$\langle
n|C_{x}|0\rangle=\delta_{n0}\left(N+\sum_{k}v_{k}^{2}\right)-2i\sum_{k}u_{k}v_{k}\langle
n|\gamma_{k}^{\dagger}\gamma_{-k}^{\dagger}|0\rangle\,,$ (181)

with [45][Chap. 10]

$\displaystyle u_{k}=\sin(\theta_{k}/2)\,,$ (182) $\displaystyle
v_{k}=\cos(\theta_{k}/2)\,,$ (183) $\displaystyle\tan\theta_{k}=\frac{\sin
k}{\frac{\tilde{\omega}_{x}}{2J}-\cos k}\,.$ (184)

With this, we find

$\tilde{\chi}_{xx,0}=-\frac{16J^{2}}{N}\sum_{k}\frac{\sin^{2}k}{\epsilon_{k}(\omega_{+}^{2}-4\epsilon_{k}^{2})}$
(185)

and in the thermodynamic limit,
$\lim_{N\to\infty}N^{-1}\sum_{k}f_{k}=\int_{-\pi}^{\pi}dk/(2\pi)f_{k}$, we
obtain Eq. (140).

## References

  * Altland and Simons [2010] A. Altland and B. D. Simons,  _Condensed Matter Field Theory_ (Cambridge University Press, 2010).
  * Schlawin  _et al._ [2022] F. Schlawin, D. M. Kennes, and M. A. Sentef, Cavity quantum materials, Appl. Phys. Rev. 9, 011312 (2022).
  * Garcia-Vidal  _et al._ [2021] F. J. Garcia-Vidal, C. Ciuti, and T. W. Ebbesen, Manipulating matter by strong coupling to vacuum fields, Science 373, eabd0336 (2021).
  * Bloch  _et al._ [2022] J. Bloch, A. Cavalleri, V. Galitski, M. Hafezi, and A. Rubio, Strongly correlated electron–photon systems, Nature 606, 41–48 (2022).
  * Hepp and Lieb [1973] K. Hepp and E. H. Lieb, On the superradiant phase transition for molecules in a quantized radiation field: the dicke maser model, Ann. Phys 76, 360 (1973).
  * Wang and Hioe [1973] Y. K. Wang and F. T. Hioe, Phase transition in the dicke model of superradiance, Phys. Rev. A 7, 831 (1973).
  * Andolina  _et al._ [2019] G. M. Andolina, F. M. D. Pellegrino, V. Giovannetti, A. H. MacDonald, and M. Polini, Cavity quantum electrodynamics of strongly correlated electron systems: A no-go theorem for photon condensation, Phys. Rev. B 100, 121109 (2019).
  * Román-Roche and Zueco [2022] J. Román-Roche and D. Zueco, Effective theory for matter in non-perturbative cavity QED, SciPost Phys. Lect. Notes , 50 (2022).
  * Moshe and Zinn-Justin [2003] M. Moshe and J. Zinn-Justin, Quantum field theory in the large n limit: a review, Phys. Rep. 385, 69–228 (2003).
  * Paravicini-Bagliani  _et al._ [2018] G. L. Paravicini-Bagliani, F. Appugliese, E. Richter, F. Valmorra, J. Keller, M. Beck, N. Bartolo, C. Rössler, T. Ihn, K. Ensslin, C. Ciuti, G. Scalari, and J. Faist, Magneto-transport controlled by landau polariton states, Nat. Phys. 15, 186–190 (2018).
  * Appugliese  _et al._ [2022] F. Appugliese, J. Enkner, G. L. Paravicini-Bagliani, M. Beck, C. Reichl, W. Wegscheider, G. Scalari, C. Ciuti, and J. Faist, Breakdown of topological protection by cavity vacuum fields in the integer quantum hall effect, Science 375, 1030–1034 (2022).
  * Thomas  _et al._ [2021] A. Thomas, E. Devaux, K. Nagarajan, G. Rogez, M. Seidel, F. Richard, C. Genet, M. Drillon, and T. W. Ebbesen, Large enhancement of ferromagnetism under a collective strong coupling of ybco nanoparticles, Nano Lett. 21, 4365–4370 (2021).
  * Jarc  _et al._ [2023] G. Jarc, S. Y. Mathengattil, A. Montanaro, F. Giusti, E. M. Rigoni, R. Sergo, F. Fassioli, S. Winnerl, S. Dal Zilio, D. Mihailovic, P. Prelovšek, M. Eckstein, and D. Fausti, Cavity-mediated thermal control of metal-to-insulator transition in 1t-tas2, Nature 622, 487–492 (2023).
  * Román-Roche  _et al._ [2024] J. Román-Roche, Álvaro Gómez-León, F. Luis, and D. Zueco, Linear response theory for cavity qed materials (2024), arXiv:2406.11957 .
  * Rokaj  _et al._ [2022a] V. Rokaj, M. Ruggenthaler, F. G. Eich, and A. Rubio, Free electron gas in cavity quantum electrodynamics, Phys. Rev. Res. 4, 013012 (2022a).
  * Li  _et al._ [2020] J. Li, D. Golez, G. Mazza, A. J. Millis, A. Georges, and M. Eckstein, Electromagnetic coupling in tight-binding models for strongly correlated light and matter, Phys. Rev. B 101, 205140 (2020).
  * Arwas and Ciuti [2023] G. Arwas and C. Ciuti, Quantum electron transport controlled by cavity vacuum fields, Phys. Rev. B 107, 045425 (2023).
  * Dmytruk and Schiró [2021] O. Dmytruk and M. Schiró, Gauge fixing for strongly correlated electrons coupled to quantum light, Phys. Rev. B 103, 075131 (2021).
  * Schuler  _et al._ [2020] M. Schuler, D. D. Bernardis, A. M. Läuchli, and P. Rabl, The vacua of dipolar cavity quantum electrodynamics, SciPost Phys. 9, 066 (2020).
  * Rokaj  _et al._ [2018] V. Rokaj, D. M. Welakuh, M. Ruggenthaler, and A. Rubio, Light–matter interaction in the long-wavelength limit: no ground-state without dipole self-energy, J. Phys. B: At. Mol. Opt. Phys. 51, 034005 (2018).
  * Libersky  _et al._ [2021] M. Libersky, R. D. McKenzie, D. M. Silevitch, P. C. E. Stamp, and T. F. Rosenbaum, Direct observation of collective electronuclear modes about a quantum critical point, Phys. Rev. Lett. 127, 207202 (2021).
  * McKenzie  _et al._ [2022] R. D. McKenzie, M. Libersky, D. M. Silevitch, and T. F. Rosenbaum, Theory of magnon polaritons in quantum ising materials, Phys. Rev. A 106, 043716 (2022).
  * Román-Roche  _et al._ [2021] J. Román-Roche, F. Luis, and D. Zueco, Photon condensation and enhanced magnetism in cavity qed, Phys. Rev. Lett. 127, 167201 (2021).
  * Feynman and Vernon [1963] R. Feynman and F. Vernon, The theory of a general quantum system interacting with a linear dissipative system, Ann. Phys. (N. Y.) 24, 118–173 (1963).
  * Grabert  _et al._ [1988] H. Grabert, P. Schramm, and G.-L. Ingold, Quantum brownian motion: The functional integral approach, Phys. Rep. 168, 115 (1988).
  * Lenk  _et al._ [2022a] K. Lenk, J. Li, P. Werner, and M. Eckstein, Collective theory for an interacting solid in a single-mode cavity (2022a), arXiv:2205.05559 .
  * Gardiner and Zoller [2004] C. Gardiner and P. Zoller,  _Quantum Noise: A Handbook of Markovian and Non-Markovian Quantum Stochastic Methods with Applications to Quantum Optics_, Springer Series in Synergetics (Springer, 2004).
  * del Valle  _et al._ [2009] E. del Valle, F. P. Laussy, and C. Tejedor, Luminescence spectra of quantum dots in microcavities. ii. fermions, Phys. Rev. B 79, 235326 (2009).
  * Carollo and Lesanovsky [2021] F. Carollo and I. Lesanovsky, Exactness of mean-field equations for open dicke models with an application to pattern retrieval dynamics, Phys. Rev. Lett. 126, 230601 (2021).
  * Lee and Johnson [2004] C. F. Lee and N. F. Johnson, First-order superradiant phase transitions in a multiqubit cavity system, Phys. Rev. Lett. 93, 083001 (2004).
  * Lenk  _et al._ [2022b] K. Lenk, J. Li, P. Werner, and M. Eckstein, Dynamical mean-field study of a photon-mediated ferroelectric phase transition, Phys. Rev. B 106, 245124 (2022b).
  * Ciuti [2021] C. Ciuti, Cavity-mediated electron hopping in disordered quantum hall systems, Phys. Rev. B 104, 155307 (2021).
  * Rokaj  _et al._ [2022b] V. Rokaj, M. Penz, M. A. Sentef, M. Ruggenthaler, and A. Rubio, Polaritonic hofstadter butterfly and cavity control of the quantized hall conductance, Phys. Rev. B 105, 205424 (2022b).
  * Rokaj  _et al._ [2023] V. Rokaj, J. Wang, J. Sous, M. Penz, M. Ruggenthaler, and A. Rubio, Weakened topological protection of the quantum hall effect in a cavity, Phys. Rev. Lett. 131, 196602 (2023).
  * von Klitzing [1986] K. von Klitzing, The quantized hall effect, Rev. Mod. Phys. 58, 519 (1986).
  * von Klitzing  _et al._ [2020] K. von Klitzing, T. Chakraborty, P. Kim, V. Madhavan, X. Dai, J. McIver, Y. Tokura, L. Savary, D. Smirnova, A. M. Rey, C. Felser, J. Gooth, and X. Qi, 40 years of the quantum hall effect, Nat. Rev. Phys. 2, 397–401 (2020).
  * Kubo [1957] R. Kubo, Statistical-mechanical theory of irreversible processes. i. general theory and simple applications to magnetic and conduction problems, J. Phys. Soc. Jpn. 12, 570–586 (1957).
  * Giuliani and Vignale [2005] G. Giuliani and G. Vignale,  _Quantum Theory of the Electron Liquid_ (Cambridge University Press, 2005).
  * Tong [2016] D. Tong, Lectures on the quantum hall effect (2016), arXiv:1606.06687 [hep-th] .
  * Dicke [1954] R. H. Dicke, Coherence in spontaneous radiation processes, Phys. Rev 93, 99 (1954).
  * Emary and Brandes [2003] C. Emary and T. Brandes, Chaos and the quantum phase transition in the dicke model, Phys. Rev. E 67, 066203 (2003).
  * Gammelmark and Mølmer [2011] S. Gammelmark and K. Mølmer, Phase transitions and heisenberg limited metrology in an ising chain interacting with a single-mode cavity field, New J. Phys. 13, 053035 (2011).
  * Cortese  _et al._ [2017] E. Cortese, L. Garziano, and S. De Liberato, Polariton spectrum of the dicke-ising model, Phys. Rev. A 96, 053861 (2017).
  * Rohn  _et al._ [2020] J. Rohn, M. Hörmann, C. Genes, and K. P. Schmidt, Ising model in a light-induced quantized transverse field, Phys. Rev. Res. 2, 023131 (2020).
  * Sachdev [2011] S. Sachdev,  _Quantum Phase Transitions_ , 2nd ed. (Cambridge University Press, 2011).
  * Rollano  _et al._ [2022] V. Rollano, M. C. de Ory, C. D. Buch, M. Rubín-Osanz, D. Zueco, C. Sánchez-Azqueta, A. Chiesa, D. Granados, S. Carretta, A. Gomez, S. Piligkos, and F. Luis, High cooperativity coupling to nuclear spins on a circuit quantum electrodynamics architecture, Commun. Phys. 5, 10.1038/s42005-022-01017-8 (2022).
  * MacNeill  _et al._ [2019] D. MacNeill, J. T. Hou, D. R. Klein, P. Zhang, P. Jarillo-Herrero, and L. Liu, Gigahertz frequency antiferromagnetic resonance and strong magnon-magnon coupling in the layered crystal ${\mathrm{crcl}}_{3}$, Phys. Rev. Lett. 123, 047204 (2019).
  * Zollitsch  _et al._ [2023] C. W. Zollitsch, S. Khan, V. T. T. Nam, I. A. Verzhbitskiy, D. Sagkovits, J. O’Sullivan, O. W. Kennedy, M. Strungaru, E. J. G. Santos, J. J. L. Morton, G. Eda, and H. Kurebayashi, Probing spin dynamics of ultra-thin van der waals magnets via photon-magnon coupling, Nat. Commun. 14, 2619 (2023).
  * Mergenthaler  _et al._ [2017] M. Mergenthaler, J. Liu, J. J. Le Roy, N. Ares, A. L. Thompson, L. Bogani, F. Luis, S. J. Blundell, T. Lancaster, A. Ardavan, G. A. D. Briggs, P. J. Leek, and E. A. Laird, Strong coupling of microwave photons to antiferromagnetic fluctuations in an organic magnet, Phys. Rev. Lett. 119, 147701 (2017).
  * Martínez-Losa del Rincón  _et al._ [2023] S. Martínez-Losa del Rincón, I. Gimeno, J. Pérez-Bailón, V. Rollano, F. Luis, D. Zueco, and M. J. Martínez-Pérez, Measuring the magnon-photon coupling in shaped ferromagnets: Tuning of the resonance frequency, Phys. Rev. Appl. 19, 014002 (2023).
  * Queisser and Schützhold [2023] F. Queisser and R. Schützhold, Hierarchy of double-time correlations, J. Stat. Mech.: Theory Exp 2023, 053101 (2023).
  * Gómez-León [2016] A. Gómez-León, Hierarchy of correlations: Application to green’s functions and interacting topological phases, Phys. Rev. B 94, 035144 (2016).
  * Gómez-León [2017] A. Gómez-León, Hierarchy of correlations for the ising model in the majorana representation, Phys. Rev. B 96, 064426 (2017).
  * Pérez-González  _et al._ [2022] B. Pérez-González, Álvaro Gómez-León, and G. Platero, Topology detection in cavity qed, Phys. Chem. Chem. Phys. 24, 15860–15870 (2022).
  * Dmytruk and Schirò [2022] O. Dmytruk and M. Schirò, Controlling topological phases of matter with quantum light, Commun. Phys 5, 271 (2022).
  * Pérez-González  _et al._ [2023] B. Pérez-González, G. Platero, and Álvaro Gómez-León, Light-matter correlations in quantum floquet engineering (2023), arXiv:2302.12290 .
  * Gómez-León [2019] A. Gómez-León, Spin bath dynamics and dynamical renormalization group, Phys. Rev. B 100, 094308 (2019).