# Fate of the Mollow triplet in strongly-coupled atomic arrays

Orazio Scarlatella  Nigel R. Cooper  T.C.M. Group, Cavendish Laboratory,
University of Cambridge, J.J. Thomson Avenue, Cambridge CB3 0HE, UK

###### Abstract

Subwavelength arrays of quantum emitters have emerged as an interesting
platform displaying prominent collective effects. Here we study such arrays
under coherent driving, realizing an open quantum many-body problem in a
strongly non-linear regime. We show that the combination of dipolar
interactions and regular geometry have a dramatic effect on the spectrum of
emitted light: the famous Mollow triplet characterizing the emission of a
single atom develops a structured broadening with flat sidebands, with a
bandwidth determined by the dipolar interactions. This emission spectrum
characterizes atomic arrays and distinguishes them from disordered ensembles
as well as non-interacting emitters. Our predictions are based on a novel
dynamical mean-field theory (DMFT) approach to the problem, paving the way for
further studies of these systems.

The collective behaviour of ensembles of quantum emitters, resulting from the
coherent scattering of photons between them, has long attracted attention
[agarwalAgarwal1970, dickeDicke1954, lehmbergLehmberg1970]. Considering two-
level atoms, when the average inter-atomic distance becomes comparable to the
wavelength of the two-level transition, the resulting dipolar interactions
become strong and can lead to collective effects, such as superradiance and
subradiance [dickeDicke1954, grossHaroche1982a] and Lamb shifts
[friedbergManassah1973, keaveneyAdams2012].

Cold atoms systems are natural places where strong dipolar interactions can be
studied and exploited. There, the emitters can be arranged in ordered
configurations, reducing dephasing and enhancing collective effects
[bettlesAdams2015, glicensteinBrowaeys2020, srakaewZeiher2023]. Such arrays of
quantum emitters have recently been investigated for several applications in
quantum information and optics [pattiYelin2021, solntsevKivshar2021],
including realizing single atomic layer mirrors [bettlesAdams2016,
shahmoonYelin2017, ruiBloch2020] with coherent control [srakaewZeiher2023] for
photonic state preparation and photon-photon gates [bekensteinLukin2020], as
well as for photon storage [cechOlmos2023, asenjo-garciaChang2017,
facchinettiRuostekoski2016]. Several effects have been studied including the
shift of collective resonances [glicensteinBrowaeys2020, jenkinsBrowaeys2016],
electromagnetically induced transparency [fengScully2017, bettlesAdams2015,
srakaewZeiher2023] and superradiance and subradiance [ostermannYelin2023,
rubies-bigordaYelin2022a, massonAsenjo-Garcia2022, rubies-bigordaYelin2021,
asenjo-garciaChang2017]. Nevertheless, most works so far focused on linear
regimes, such as weakly-driven regimes in which the emitters behave as
classical scatterers [williamsonRuostekoski2020]. Less is known when non-
linearities become important, for which the system constitutes a far-from-
equilibrium quantum many-body problem with strong interactions.

Here we consider atomic arrays illuminated by a strong plane wave field and in
its driven-dissipative steady-state. Previous works mapped out a mean-field
phase diagram [parmeeCooper2018] and investigated the fate of steady-state
superradiance [ostermannYelin2023], while other works were mostly limited to
small system sizes [yuYu2016, jonesOlmos2017, olmosLesanovsky2014,
robicheauxSuresh2021].

Few theoretical methods are available in these nonlinear regimes, and the most
used approaches are generalized mean-field approaches based on a factorization
of spatial cumulants [kramerRitsch2015, plankensteinerRitsch2022]. These have
been successful in describing the collective decay of atomic arrays [rubies-
bigordaYelin2021, rubies-bigordaYelin2022a, rubies-bigordaYelin2023] and have
been recently applied to steady-state cases [ostermannYelin2023,
robicheauxSuresh2021], but become inaccurate for strong correlations such that
higher-order cumulants become important [robicheauxSuresh2021,
kramerRitsch2015]. Oh the other hand, a Dynamical Mean Field Theory (DMFT)
approach has recently been adapted for lattice models described by Markovian
master equations, as they commonly arise in atomic physics
[scarlatellaSchiro2021]. This technique, describing the limit of large lattice
coordination, has been in fact extremely successful in describing condensed
matter systems [georgesRozenberg1996], including in non-equilibrium conditions
[aokiWerner2014].

In this work we leverage the approach of [scarlatellaSchiro2021] and set up a
non-equilibrium DMFT technique for the steady states of driven-dissipative
arrays of two-level emitters, on the lines of a recent DMFT approach for
light-matter problems in equilibrium [lenkEckstein2022a], paving the way for
investigating strongly-correlated phenomena in these systems.

Using this approach, we study the spectrum of light emitted by arrays of close
two-level atoms strongly driven on resonance. For a single atom, this is
characterized by the famous Mollow triplet [mollowMollow1969], due to an atom-
photon dressing [cohen-tannoudjiReynaud1977]. How this is affected by dipolar
interactions in the case of a few atoms is well known [agarwalVetri1977,
senitzkySenitzky1978, carmichaelCarmichael1979, freedhoffFreedhoff1979,
kilinKilin1980]. Less is known instead for large ensembles, where recent works
investigated disordered clouds [ottKaiser2013, pucciBachelard2017,
jonesOlmos2017], while regular geometries were only investigated for
particular emission directions [jonesOlmos2017] or for cases of few atoms
[moniriDarsheshdar2023].

Here we find that for large arrays of close atoms in 1D and 2D square lattice
geometries, and regardless of the details of the drive and for generic
emission directions, the triplet acquires a structured broadening
characterized by flat sidebands. Those sidebands have a width determined by
the coherent dipolar interactions, in contrast with the single-atom case in
which this is set by the atom decay rate, and with the case of disordered
ensembles where some broadening was predicted but without any structure
[jonesOlmos2017]. This structured broadening of the triplet therefore
characterizes atomic arrays and distinguishes their emission from that of
disordered ensembles or of non-interacting emitters.

Figure 1: A 1D array of two-level atoms driven by an electric field. The atoms
interact via dipole-dipole interactions and collective dissipation. The
photons emitted are collected by a detector.

Model – We consider a large number $N$ of two-level atoms ordered in a
periodic array. The atoms are illuminated with a uniform plane wave and are
coupled to the free-space electromagnetic field, giving rise both to coherent
dipole-dipole interactions and to collective dissipation. These can be
described by means of a Markovian master equation [lehmbergLehmberg1970]:

$\displaystyle\dot{\rho}$ $\displaystyle=-\frac{{\rm i}}{\hbar}[H_{0}+H_{\rm
drive},\rho]+\mathcal{D}[\rho]$ (1) $\displaystyle H_{0}$
$\displaystyle=\frac{\hbar\Delta}{2}\sum_{i}^{N}{\sigma}_{i}^{z}+\sum_{ij;i\neq
j}^{N}\hbar V_{ij}\sigma_{i}^{+}\sigma_{j}^{-}$ (2) $\displaystyle H_{\rm
drive}$
$\displaystyle=\frac{\hbar\Omega}{2}\sum_{i}^{N}\left({\sigma}_{i}^{+}e^{{\rm
i}\mathbf{k}_{d}\cdot\mathbf{r}_{i}}+{\sigma}_{i}^{-}e^{-{\rm
i}\mathbf{k}_{d}\cdot\mathbf{r}_{i}}\right)$ (3)
$\displaystyle\mathcal{D}[\rho]$
$\displaystyle=\frac{1}{2}\sum_{i,j}\Gamma_{ij}\left(2\sigma_{i}^{-}\rho\sigma_{j}^{+}-\sigma_{i}^{+}\sigma_{j}^{-}\rho-\rho\sigma_{i}^{+}\sigma_{j}^{-}\right).$
(4)

Here $\sigma_{l}^{\alpha}$ with $\alpha=x,y$ or $z$ are the Pauli matrices on
site $l$ and ${\sigma}_{l}^{\pm}={\sigma}_{l}^{x}\pm{\rm i}{\sigma}_{l}^{y}$
the raising and lowering operators; $\Delta=\omega_{d}-\omega_{0}$ is the
detuning of the drive from the two-level transition energy $\omega_{0}$,
$\Omega=2\mathbf{d}\cdot\mathbf{E}/\hbar$ the Rabi coupling given by the
vector of transition dipole moments $\mathbf{d}$ and the electric field vector
$\mathbf{E}$, and $\mathbf{k}_{d}$ is the drive momentum.

The dipolar interactions and collective decay terms are given by
[kramerRitsch2015, asenjo-garciaChang2017]

$\displaystyle
V_{il}=\frac{3\Gamma}{4}\left\\{-\alpha\frac{\cos(k_{0}r_{ij})}{(k_{0}r_{ij})}+\beta\left[\frac{\sin(k_{0}r_{ij})}{(k_{0}r_{ij})^{2}}+\frac{\cos(k_{0}r_{ij})}{(k_{0}r_{ij})^{3}}\right]\right\\}$
(5)
$\displaystyle\Gamma_{il}=\frac{3\Gamma}{2}\left\\{\alpha\frac{\sin(k_{0}r_{ij})}{(k_{0}r_{ij})}+\beta\left[\frac{\cos(k_{0}r_{ij})}{(k_{0}r_{ij})^{2}}-\frac{\sin(k_{0}r_{ij})}{(k_{0}r_{ij})^{3}}\right]\right\\}$
(6)

where the single-atom decay rate is $\Gamma=$
$|\mathbf{d}|^{2}k_{0}^{3}/3\pi\epsilon_{0}\hbar$ with
$r_{il}=\left|\mathbf{r}_{i}-\mathbf{r}_{l}\right|=a|i-l|$, $\mathbf{r}_{i}$
is the position of atom $i$ on the lattice and $a$ the lattice spacing.
$k_{0}=2\pi/\lambda_{0}$ is the wavevector of two-level transition, with
wavelength $\lambda_{0}$. $\alpha=1-\cos^{2}\theta$ and
$\beta=1-3\cos^{2}\theta$, where $\theta$ represents the angle between the
line connecting atoms $i$ and $j$ and the common atomic dipole orientation.

Throughout the manuscript, we will assume units of $\Gamma=1$ and $\hbar=1$.
The interactions and losses are then completely determined by the parameter
$k_{0}a$: in the extreme limit $k_{0}a\rightarrow\infty$ interactions vanish
and dissipation becomes local with $\Gamma_{il}=\Gamma\delta_{il}$, recovering
the independent atoms case, while for $k_{0}a\rightarrow 0$ interactions
diverge while dissipation becomes negligible, going to a finite value while
also acquiring an all-to-all character $\Gamma_{il}=\Gamma$. Here we
concentrate on intermediate values (see [sm] for a plot).

DMFT for atomic arrays – To apply DMFT to the present case of a spin model,
which lacks a gaussian action describing a non-interacting limit, we decouple
the dipolar interactions and collective dissipation by introducing auxiliary
bosonic fields and apply a bosonic DMFT approach to those fields
[byczukVollhardt2008, huTong2009, andersWerner2010, andersWerner2011,
strandWerner2015, scarlatellaSchiro2021], on the lines of [lenkEckstein2022a]
and of an extended DMFT approach [sunKotliar2002].

Bosonic DMFT maps a lattice model onto an effective impurity model, describing
a single site of the lattice coupled to an effective field, as in a standard
Gutzwiller mean-field approach, and to an effective environment. Assuming a
local self-energy in space, both can be determined such as to satisfy a set of
self-consistent equations. The effective environment introduces corrections of
order $1/z$, where $z$ is the lattice coordination number, to Gutzwiller mean-
field corresponding to $1/z\rightarrow 0$ [strandWerner2015]. In the present
case in which interactions are not limited to nearest neighbors, the accuracy
of DMFT is expected to improve. We later confirm our results by comparing with
exact diagonalization (ED).

In the present case, the impurity model is a generalized spin-boson model,
that in the steady-state is described by the time-translation invariant
Keldysh action

$\begin{split}S^{\mathrm{imp}}&=S_{0}^{i}+\frac{1}{2}\int_{-\infty}^{\infty}dt\left[\varsigma_{i}^{\dagger}(t)\tau_{3}b+\text{hc}\right]\\\
&-\frac{1}{2}\int_{-\infty}^{\infty}dt\int_{-\infty}^{\infty}dt^{\prime}\varsigma_{i}^{\dagger}(t)\tau_{3}{\mathcal{W}}\left(t-t^{\prime}\right)\tau_{3}\varsigma_{i}\left(t^{\prime}\right)\end{split}$
(7)

Here $\varsigma_{i}=(d_{i+},\bar{d}_{i+},d_{i-},\bar{d}_{i-})^{T}$ is a vector
in Nambu and Keldysh formalisms, where $d$ and $\bar{d}$ are the scalar fields
corresponding to the spin operators $\sigma^{-}$ and $\sigma^{+}$ where the
bar indicates complex conjugation and their $+$ or $-$ indices indicate to
which branch of the Keldysh double contour they belong to. $b$ and
$\mathcal{W}$ represent the effective field and effective environment. These
are taken site-independent assuming a homogeneous steady-state, as verified by
a linear stability analysis. $\tau_{3}$ is a diagonal matrix with entries
$(1,1,-1,-1)$ arising in the Keldysh-Nambu formalism. Note that here a Nambu
formalism is needed to allow for finite values of the anomalous correlators
$\left\langle\sigma^{\pm}\sigma^{\pm}\right\rangle$, that arise due to the
drive term, breaking the $U(1)$ symmetry of the undriven problem.

The self-consistent conditions correspond to a set of matrix relations in
frequency and momentum space:
$\Pi_{\mathrm{loc}}=\left[\tau_{3}{\chi}^{-1}_{ii}\tau_{3}+\mathcal{W}\right]^{-1}$,
$U_{ii}=\frac{1}{N}\sum_{k}\left[W_{k}^{-1}-\Pi_{\mathrm{loc}}\right]^{-1}$
and $\mathcal{W}^{-1}=\Pi_{\mathrm{loc}}+U_{ii}^{-1}$ , and
$b_{\pm}=\left(\mathcal{W}^{R}(\omega={0})-W_{k=0}^{R}(\omega=0)\right)\langle{\varsigma_{i\pm}}\rangle$.
Here the superscript “$R$” indicates the retarded component and $b_{\pm}$ and
$\langle{\varsigma_{i\pm}}\rangle$ the $+$ or $-$ components of the
corresponding vectors. $W_{k}$ is a known matrix with entries defined by the
Fourier transforms of (5) and (6), reported in [sm]. Assuming one can compute
from the impurity action (7) the local connected Green’s function
$\chi_{ii}(\tau)=-{\rm
i}\lim_{t\rightarrow\infty}\left\langle\varsigma_{i}(t+\tau)\varsigma^{\dagger}_{i}(t)\right\rangle_{S_{\text{imp
}}^{\text{con }}}$ and expectation value
$\left\langle\varsigma_{i}\right\rangle$, then these equations determine
$\mathcal{W},b,\Pi_{\mathrm{loc}}$ and $U_{ii}$, such that the impurity
problem represents the original lattice problem.

The connected propagator of the lattice $\chi_{ij}(\tau)=-{\rm
i}\lim_{t\rightarrow\infty}\left\langle\varsigma_{i}(t+\tau)\varsigma^{\dagger}_{j}(t)\right\rangle^{\textrm{con}}$
can then be computed by $\chi^{-1}(k,\omega)=\tau_{3}(\Pi_{\rm
loc}^{-1}(\omega)-W_{k})\tau_{3}$.

Solving the impurity model is a challenge by itself and a viable numerically
exact approach is still lacking. Here we use a method based on a non-crossing
approximation (NCA) [schiroScarlatella2019, scarlatellaSchiro2021,
scarlatellaSchiro2024a]. More details on DMFT/NCA are found in [sm]. As a
technical remark, due to the Markovian component of the environment, one has
to carefully discretize the Green’s functions at $t=0$ such that the Keldsyh
causality structure [kamenevKamenev2023] is preserved.

Figure 2: (a) Sketch of the single-particle transitions determining the Mollow
triplet. In the case of many atoms (bottom) they take place between two bands
determined by the dipolar interactions. (b) The local atomic correlation
function for a 1D chain (solid) and 2D lattice (dashed) as a function of
lattice spacing and in DMFT. (c) The width of the sidebands, extracted by the
extrema of the derivatives in 1D (dashed lines in panel (b)), scales with the
nearest-neighbor interaction $V_{i,i+1}$ as a function of lattice spacing.
Here $\Omega=60\Gamma$, $\theta=0$ and
$\mathbf{k}_{d}~{}\cdot~{}\mathbf{r}_{i}=0$.

Results –  The spectrum of the emitted light from a single two-level atom
which is driven on resonance and strongly $\Omega\gg\Gamma$ shows a central
peak at the two-level transition energy $\omega=\omega_{0}$, accompanied two
side peaks at $\omega=\omega_{0}\pm\Omega$, known as the Mollow triplet
[mollowMollow1969]. This is understood by a picture in which the atom is
dressed by the laser field [cohen-tannoudjiReynaud1977, scullyZubairy1997].
The three transition energies correspond to incoherent transitions induced by
the jump operators between the eigenstates $|+\rangle$,$|-\rangle$ of the
Hamiltonian (2) (3) in the case of a single atom, where it reduces to
$H=\Omega\sigma^{x}/2$. These are represented graphically in Fig. 2 (a), where
as (2) is defined in the frame rotating at $\omega_{0}$, the resonances are
shifted at $\omega=0$ and $\omega=\pm\Omega$.

In the case of several atoms and for $V_{ij}\ll\Omega$, single-particle
transitions can be understood as taking place between two bands of width set
by $V_{ij}$, separated by the drive intensity $\Omega$, as depicted in Fig. 2.
The width of the bands corresponds to a broadening of the triplet lines, where
interband transitions correspond to the central peak and intraband ones to the
sidebands. The spectrum of emitted light in the steady-state is related to the
atomic correlation functions
$g_{ij}(\tau)=\lim_{t\rightarrow\infty}\langle{\sigma_{i}^{+}(t+\tau)\sigma_{j}^{-}(t)}\rangle$
[lehmbergLehmberg1970].

Predicting how triplet broadens is nevertheless non-trivial. No broadening is
predicted by a Gutzwiller mean-field approach, as well as by approximating the
lattice self-energy with that of a single site, that allows to analytically
estimate the Hubbard bands of a Mott insulator [georgesRozenberg1996,
strandWerner2015a]. This motivates a DMFT approach.

In DMFT, a natural quantity to compute is the local Green’s function
$g_{ii}(\tau)$. This could be measured by near-field microscopy
[courjonBainier1994, courjonCourjon2003]. It is shown in Fig. 2 panel (b) for
different values of $k_{0}a$, for a 1D chain and 2D square lattice. We assume
there a perpendicular drive $\mathbf{k}_{d}~{}\cdot~{}\mathbf{r}_{i}=0$ and a
drive polarization corresponding to $\theta=0$, though similar results are
found for different choices. For large $k_{0}a$ one recovers the single-atom
triplet, composed by three Lorentzian peaks, while these transform into more
complicated structures at smaller $k_{0}a$, eventually merging when
$\Omega\sim V_{ij}$. The broadening has a remarkable structure with the
sidebands turning into flat and broad features, while the central feature
remaining peaked at its center. While for a single atom the peak widths are
determined by the decay rate $\Gamma$ [scullyZubairy1997], for the arrays the
sideband widths are set by the dipolar interactions: Fig. 2 (c) shows that
this scales as the nearest-neighbor interaction $V_{i,i+1}$, which goes like
$1/(k_{0}a)^{3}$ in the corresponding range of values. In this regime of small
$k_{0}a$ in fact the dipolar interaction dominates over dissipation, as we
further confirm in [sm] by modulating the dissipation strength. As we will
demonstrate throughout the paper, a broadening with this structure occurs in a
variety of conditions and constitutes a hallmark of the emission properties of
these atomic arrays.

The spectrum of emitted light in the far field is in fact proportional to the
momentum-resolved propagator $g(\mathbf{k},\omega)$, computed at
$\mathbf{k}=k_{0}{\mathbf{n}}$ where ${\mathbf{n}}$ is a unit vector
indicating the emission direction [lehmbergLehmberg1970]. In the regime of
large drive considered, we find that this propagator becomes momentum
independent in DMFT and coincides with the local propagator
$g(\mathbf{k},\omega)\approx g_{ii}(\omega)$. This indicates that the
structured broadening of the latter also characterizes the far-field emission.

Since DMFT assumes a momentum-independent self-energy, we confirm our DMFT
predictions with exact diagonalization (ED) of a small 1D chain with open
boundary conditions. To access larger systems sizes, we disregard the
dissipator as the dipolar interaction is dominant for small $k_{0}a$ as
discussed. Also, in this large-drive regime the stationary-state is
featureless, approaching infinite temperature. The local Green’s function can
then be approximated by weighting the density of states (DoS)
$\sum_{n,m}\delta({\omega-E_{n}-E_{m}})$ with the matrix elements
$|\left[\sigma_{i}^{-}\right]_{nm}|^{2}=|\langle\psi_{n}|\sigma_{i}^{-}|\psi_{m}\rangle|^{2}$,
see e.g. [scarlatellaSchiro2019a], namely $g_{ii}(\omega)\approx
2\pi/2^{N}\sum_{n,m}|\left[\sigma_{i}^{-}\right]_{nm}|^{2}\delta({\omega-
E_{n}-E_{m}})$, where $|\psi_{n}\rangle$ and $E_{n}$ are the eigenstates and
eigenvalues of the Hamiltonian. The latter is compared with the DMFT local
Green’s function in panel (a) of Fig. 3 in the top row, displaying very good
agreement. We note that such an agreement in 1D is remarkable, since DMFT best
describes large dimensional systems. A possible understanding is that the
physics at play is very short-ranged, as we find that replacing the long range
interactions and dissipation (5)(6) with short range ones our results are
essentially unchanged [sm]. The DoS, in the 3rd row, is similar in
correspondence of the three peaks, corresponding to an equal width of the
bands sketched in Fig. 2. The structured broadening of the Green’s function
therefore reveals some structure of the many-body eigenstates beyond the DoS,
that is captured by the matrix element in the 2nd row.

Figure 3: A comparison of the DMFT Green’s functions with an approximation
from exact diagonalization (ED) of the Hamlitonian for a chain of $N=10$
sites. In panel (a) the DMFT local Green’s function (top row) is compared with
a combination of the relevant matrix element and density of states, calculated
as a histogram; these are also shown in the 2nd and 3rd row. In panel (b) the
DMFT momentum-resolved Green’s function $g(k_{z},\omega)$ is compared to ED in
an analogous approximation for $k_{0}a=0.7$. Here $\Omega=60\Gamma$,
$\theta=0$ and $\mathbf{k}_{d}~{}\cdot~{}\mathbf{r}_{i}=0$.

Leveraging this understanding, we approximate also the momentum-resolved
Green’s function with a weighted sum of the DoS with matrix elements
$|\left[\sigma_{k}^{-}\right]_{nm}|^{2}/N=|\langle\psi_{n}|\sigma_{k}^{-}|\psi_{m}\rangle|^{2}/N$
in panel (b) of Fig. 3. This compares well with the momentum-independent DMFT
prediction at momenta away from $k_{z}=0$. Approaching $k_{z}=0$ instead the
resonances narrow down, such that neglecting the dissipator is not justified
there, and comparison with DMFT deteriorates. This corresponds to a momentum
dependence that is missed in DMFT, rather than a finite-size effect [sm]. The
perfectly narrow matrix element for the central peak at $k_{z}=0$ can be
explained by a total spin conservation argument [sm]. The narrow peaks in this
case also explain why a structured broadening was not predicted in
[jonesOlmos2017], where emission from a 1D chain has been studied, but only in
an emission direction orthogonal to the chain corresponding to $k_{z}=0$.

We remark that the far-field emission can only probe the Brillouin zone for
$\left\lvert k_{z}\right\rvert\leq k_{0}$, where $|k_{z}|=k_{0}$ corresponds
to a direction of emission parallel to the chain [lehmbergLehmberg1970].
Nevertheless, there is a finite window of momenta in which such a broadening
can be observed, as shown in Fig. 3. Changing the drive direction provides a
natural way around this constraint.

Similar spectra are in fact found for different directions, as Fig. 4 (a)
shows in 1D. A finite drive momentum along the chain $k_{d,z}$ can be
accounted for within DMFT by a gauge transformation
$\sigma_{j}^{+}e^{ik_{d,z}aj}\rightarrow\sigma_{j}^{+}$ that trades the drive
phases with a shift in momentum of dipolar interactions and dissipation
$V_{k_{z}}\rightarrow V_{k_{z}-k_{d,z}}$ and
$\Gamma_{k_{z}}\rightarrow\Gamma_{k_{z}-k_{d,z}}$. We find that in DMFT the
atoms Green’s functions are independent of the drive direction. Similarly, in
ED they change very little as shown in Fig. 4 (a), apart that they are shifted
in momentum by $k_{d,z}$. This allows to scan the Brillouin zone by varying
the drive direction at a fixed emission direction.

Similar spectra are also found across different drive polarizations $\theta$,
as shown in Fig. 4 panel (b) for the local propagator of a 1D chain.

Figure 4: (a) Momentum-resolved atomic Green’s function for different drive
directions labelled by $k_{z,d}$, approximated as in Fig. 3 in ED for a 1D
chain of $N=10$ sites and for $k_{z}a=k_{0}a=0.5$ and a polarization
$\theta=\pi/2$. The black lines corresponds instead to the DMFT prediction.
(b) Local atomic Green’s function for different drive polarizations $\theta$
for a 1D chain in DMFT and $\mathbf{k}_{d}~{}\cdot~{}\mathbf{r}_{i}=0$. Here
$\Omega=60\Gamma$.

Conclusions –  We showed that the spectrum of light emitted by atomic arrays
strongly driven on resonance is characterized by a structured broadening of
the famous Mollow triplet of a single atom, that acquires flat sidebands with
a width controlled by the dipolar interactions. This structure is predicted in
several conditions and therefore constitutes a hallmark of the emission of
atomic arrays and distinguishes it from that of disordered ensembles or non-
interacting emitters. In particular, we discussed how it can be observed in
near-field microscopy and in the far-field emission.

We have also introduced a Dynamical Mean-Field Theory approach for the steady-
states of atomic arrays, paving the way for investigating a variety of
strongly-correlated phenomena in these systems.

Acknowledgements –  We thank Mike Gunn for insightful discussions. This work
was supported by the Engineering and Physical Sciences Research Council [grant
number EP/W005484], and a Simons Investigator Award [Grant No. 511029]. For
the purpose of open access, the authors has applied a creative commons
attribution (CC BY) licence to any author accepted manuscript version arising.

The data (code) to reproduce the results of the manuscript will provided by
the author under request.

## References

## Supplemental Material to “Fate of the Mollow triplet in strongly-coupled
atomic arrays”

This supplemental is organized as follows. In App. A we give additional
details on the model, in App. B we discuss further the triplet dipendence on
parameters. In App. C we report a further discussion of exact diagonalization
results. The DMFT formalism is discussed in App. D, its numerical
implementation in App. E, the NCA impurity solver in App. F and the definition
of the Green’s functions at $t=0$ in App. G.

##  Appendix A Details of the model

The dipolar interaction and collective dissipation are plotted in Fig. S.1 for
several values of $k_{0}a$ and in 1D. Note that for small values of $k_{0}a$
the interaction becomes larger than the dissipation rate.

Figure S.1: Dipolar interactions $V_{il}$ and collective dissipation rate
$\Gamma_{il}$ for several values of $k_{0}a$ for a 1D chain and for
$\theta=0$.

##  Appendix B Other results in DMFT

In Fig S.2 we show that by artificially changing the dissipation strength
$\Gamma_{ij}\rightarrow\kappa\Gamma_{ij}$ (dotted line) while keeping the
dipole interactions constant essentially leaves the sidebands widths unchanged
and mainly affects the amplitudes, though the lesser the smaller $k_{0}a$.

We also find that switching to a local dissipation
$\Gamma_{ij}\rightarrow\delta_{ij}\Gamma_{ii}$ (dashed line) or nearest-
neighbor interactions $V_{ij}=\delta_{j,i\pm 1}V_{i,i\pm 1}$ leaves the curves
unchanged in DMFT. Similar results are also found in ED.

Figure S.2: Local atomic correlation function changing the dissipation
strength $\Gamma_{ij}\rightarrow\kappa\Gamma_{ij}$ but keeping the coherent
dipolar interaction fixed (dotted line), or by making the dissipation local
$\Gamma_{ij}\rightarrow\delta_{ij}\Gamma_{ii}$ (dashed line). A similar result
to the latter is found making the interactions short ranged
$V_{ij}\rightarrow\delta_{j,i\pm 1}V_{i,i\pm 1}$. Here $\Omega=20\Gamma$.

##  Appendix C Other results in exact diagonalization

Fig. S.3 (a) shows the momentum-resolved propagator in ED, approximated as in
main text, throughout the entire Brillouin zone. In panel (b) the widths of
the triplet resonances are extracted from the momentum-resolved matrix element
in ED by a threshold value on the amplitude. These panel show that the
momentum dependence is most pronounced around the points $k_{z}a=0$ and
$k_{z}a=\pi$ where the widths of the resonances decreases significantly and
comparison with DMFT is the worst. Panel (c) show that the same behaviour is
found changing the system sizes, suggesting that those features are not
finite-size effects, but rather a true momentum dependence that cannot be
resolved by DMFT.

Figure S.3: (a) Momentum-resolved propagator in DMFT and in ED, in the same
approximation of neglecting the dissipator as in the main text, scanning the
whole Brillouin zone. (b) and (c): the widths of the sidebands and central
peak vs momentum, extracted from the momentum-resolved matrix element, defined
by a threshold value on the amplitude of $0.05$. (a),(b) are for N=10 atoms.
(c) shows the dependence on system size $N$, respectively close to $k=0$ and
$k=\pi$.

####  C.0.1 Narrow central peak at k=0

The fact that the matrix element of $\sigma_{k=0}^{-}=\sum_{i}\sigma_{i}^{-}$
has a vanishing width for the central peak in Figs. 3 and S.3 arises due to
total spin conservation. This can be understood in perturbation theory. For
this we need to calculate the matrix element of $\sum_{i}\sigma_{i}^{-}$ on
the eigenstates of the Hamiltonian. It is useful to make a rotation such that
$\sigma^{x}\rightarrow\sigma^{z}$, for which
$\sigma^{-}=(\sigma^{x}-i\sigma^{y})/2\rightarrow(\sigma^{z}-i\sigma^{x})/2$
and the Hamiltonian reads
$H\rightarrow\frac{\Omega}{2}\sum_{i}\sigma_{i}^{z}+\sum_{ij}V_{ij}/4(\sigma_{i}^{z}+i\sigma_{i}^{x})(\sigma_{i}^{z}-i\sigma_{i}^{x})$.
Assuming $\Omega\gg V$ and treating $V$ as a perturbation, the perturbed
eigenstates at zeroth order are labelled by $|\sigma_{\rm tot}^{z},n\rangle$,
obtained by diagonalizing the perturbation in the degenerate subspaces of the
non-interacting Hamiltonian at fixed $\sigma_{\rm tot}^{z}$. In the new basis,
we want to compute $\langle\sigma_{\rm
tot1}^{z},n_{1}|\sum_{i}(\sigma_{i}^{z}-i\sigma_{i}^{x})/2|\sigma_{\rm
tot2}^{z},n_{2}\rangle$, where the transitions contributing to the central
peak are given by the case $\sigma_{\rm tot1}^{z}=\sigma_{\rm tot2}^{z}$. In
this case $\langle\sigma_{\rm tot}^{z},n_{1}|\sigma_{i}^{x}|\sigma_{\rm
tot}^{z},n_{2}\rangle=0$ because
$\sigma_{i}^{x}=\sigma_{i}^{+}+\sigma_{i}^{-}$ changes total $\sigma^{z}$ by 1
(expanding in the basis of tensor products of $\sigma^{z}_{i}$ eigenstates the
overlap vanishes term by term), while $\sum_{i}\sigma_{i}^{z}$ gets replaced
by its eigenvalue and the matrix element vanishes for $n_{1}\neq n_{2}$ by
orthogonality of different eigenstates

$\langle\sigma_{\rm
tot}^{z},n_{1}|\sum_{i}(\sigma_{i}^{z}-i\sigma_{i}^{x})/2|\sigma_{\rm
tot}^{z},n_{2}\rangle=\sigma_{\rm tot}^{z}\delta_{n_{1},n_{2}}/2$ (8)

This corresponds to the narrow central peak at $k=0$. The numerics show that
this holds beyond this simple perturbative argument. For $k\ll 1$ the pointed
shape of the central peak can be understood as those $n_{1}\neq n_{2}$ matrix
elements, corresponding to a broadening of the central peak, acquiring a
finite yet small weight compared to the $n_{1}=n_{2}$ terms.

####  C.0.2 Few atoms case

In Fig. S.4 we show the local propagator computed from an exact
diagonalization of the full Liouvillian (including the dissipator). This shows
how the broadening in the thermodynamic limit emerges from the limit of a few
atoms. We note that in the case of very few atoms that could be solved
analytically, is difficult to infer the structure in the thermodynamic limit.

Figure S.4: Green’s function computed in ED by diagonalizing the full
Liouvillian (including the dissipator) for increasing number of atoms and for
$k_{0}a=0.7$, $\Omega=60$ and for $i=N/2$ ($i=(N+1)/2$) for even (odd) $N$.

##  Appendix D Nonequilibrium steady-state DMFT

Here we present the main steps in the derivation of the DMFT approach for the
steady-state of the master equation (1).

###  D.1 Keldysh action

First we switch from the master equation (1) to a Keldysh formalism
[kamenevKamenev2023]. Markovian equations can be readily translated in Keldysh
actions that are time local, where some care has to be taken about operator
ordering [siebererDiehl2016, siebererDiehl2014]. The model is then described
in terms of the action $S_{\text{eff }}=S_{0}+S_{\text{int }}$ where the non-
local terms are given by

$\begin{split}S_{\mathrm{int}}&=-\int_{-\infty}^{\infty}dt\left[\sum_{ij,i\neq
j}V_{ij}\left(\bar{d}_{i_{+}}d_{j+}-\bar{d}_{i-}d_{j-}\right)\right.\\\
&+\left.\frac{i}{2}\sum_{i,j}\Gamma_{ij}\left(2d_{i+}\bar{d}_{j-}-\bar{d}_{i+}d_{j+}-\bar{d}_{i-}d_{j-}\right)\right]\end{split}$
(9)

Here the fields $d_{+},\bar{d}_{+},d_{-},\bar{d}_{-}$ defined in the main text
are evaluated at the integration time $t$. The explicit expression of the
local part of the action $S_{0}$, that is non-gaussian for spin models
[altlandSimons2012], is not important for our purposes.

For DMFT, it is convenient to use the Keldysh-Nambu vectors notation
introduced in the main text. To work with unrestricted sums, we define
$V_{ii}=0$, and to use the Nambu formalism we note that $(2,2)$ Nambu
component is related to the $(1,1)$ by transposing space-time and Keldysh
indices. Note also that $\Gamma,V$ are symmetric $\Gamma_{ji}=\Gamma_{ij}$,
$V_{ji}=V_{ij}$. The action then can be written as

$\begin{split}S_{\mathrm{int}}=-\frac{1}{2}\int_{-\infty}^{\infty}dt\sum_{ij}\begin{pmatrix}\vdots\end{pmatrix}_{i}^{\dagger}\begin{pmatrix}{V}_{ij}-\frac{i}{2}{\Gamma_{ij}}&0&0&0\\\
0&{V}_{ij}-\frac{i}{2}{\Gamma_{ij}}&0&i\Gamma_{ij}\\\
i{\Gamma_{ij}}&0&-{V}_{ij}-\frac{i}{2}{\Gamma_{ij}}&0\\\
0&0&0&-{V}_{ij}-\frac{i}{2}{\Gamma_{ij}}\end{pmatrix}\begin{pmatrix}d_{+}\\\
\bar{d}_{+}\\\ d_{-}\\\ \bar{d}_{-}\end{pmatrix}_{j}\end{split}$ (10)

In frequency and momentum space it becomes

$\begin{split}S_{\mathrm{int}}&=-\frac{1}{2}\int_{-\infty}^{\infty}\frac{d\omega}{2\pi}\frac{1}{N}\sum_{k}\begin{pmatrix}\vdots\end{pmatrix}^{\dagger}\begin{pmatrix}{V}_{k}-\frac{i}{2}{\Gamma_{k}}&0&0&0\\\
0&{V}_{k}-\frac{i}{2}{\Gamma_{k}}&0&i\Gamma_{k}\\\
i{\Gamma_{k}}&0&-{V}_{k}-\frac{i}{2}{\Gamma_{k}}&0\\\
0&0&0&-{V}_{k}-\frac{i}{2}{\Gamma_{k}}\end{pmatrix}\begin{pmatrix}d_{+}(\omega,k)\\\
\bar{d}_{+}(-\omega,-k)\\\ d_{-}(\omega,k)\\\
\bar{d}_{-}(-\omega,-k)\end{pmatrix}\\\
&=-\frac{1}{2}\int_{-\infty}^{\infty}\frac{d\omega}{2\pi}\frac{1}{N}\sum_{k}\begin{pmatrix}\vdots\end{pmatrix}^{\dagger}\tau_{3}W_{k}\tau_{3}\begin{pmatrix}d_{+}(\omega,k)\\\
\bar{d}_{+}(-\omega,-k)\\\ d_{-}(\omega,k)\\\
\bar{d}_{-}(-\omega,-k)\end{pmatrix}\end{split}$ (11)

which defines the matrix $W_{k}$ appearing in the main text.

To proceed, we write the real space and real time action (10) in the compact
form

$S_{\mathrm{int}}=-\frac{1}{2}\int_{-\infty}^{\infty}dtdt^{\prime}\sum_{rr^{\prime}}\varsigma_{r}^{\dagger}(t)\tau_{3}W_{rr^{\prime}}\left(t-t^{\prime}\right)\tau_{3}\varsigma_{r^{\prime}}\left(t^{\prime}\right)$
(12)

where we defined the site positions $r=ia$ and where in the present case of a
Markovian master equation
$W_{rr^{\prime}}\left(t-t^{\prime}\right)=W_{rr^{\prime}}\delta(t-t^{\prime})$
is time local.

In the following, we generalize the equilibrium DMFT approach of
[lenkEckstein2022a] to Keldysh field theory in the stationary state. Our
derivation mostly follows this reference, therefore we report the main
equations here and refer to [lenkEckstein2022a] for more details on the
derivation.

###  D.2 Bosonic representation of dipolar interaction and collective
dissipation

We introduce auxiliary bosonic fields $\phi_{r}(t)$ at each site using a
Hubbard-Stratonovich transformation. Namely we define a new action
$S_{\text{HS }}$ such that the original action is reproduced if the auxiliary
gaussian fields are integrated out: $e^{iS_{\text{eff
}}}=\int\mathcal{D}[\phi]e^{iS_{\mathrm{HS}}}.$ Up to an irrelevant constant
that we will omit throughout the manuscript, the Hubbard-Stratonovich action
reads

$S_{\mathrm{HS}}=S_{0}+S_{\phi\phi}+S_{\phi\sigma},$ (13)

with a quadratic term

$S_{\phi\phi}=\frac{1}{2}\int_{-\infty}^{\infty}dtdt^{\prime}\sum_{rr^{\prime}}\phi_{r}^{\dagger}(t)\tau_{3}\left[W^{-1}\right]_{rr^{\prime}}\left(t-t^{\prime}\right)\tau_{3}\phi_{r^{\prime}}\left(t^{\prime}\right),$

and a local linear coupling term

$S_{\phi\sigma}=\frac{1}{2}\sum_{r}\int_{-\infty}^{\infty}dt\phi^{\dagger}_{r}(t)\tau_{3}\varsigma_{r}(t)+\text{hc}$

where now $\phi$ is a Nambu vector of complex bosonic fields.

###  D.3 Impurity model and self-consistent equations

The Hubbard-Stratonovic action defines a bosonic field theory coupled to spin
fields, that introduce an anharmonicity for the bosons. Since the coupling is
local, we can apply a bosonic DMFT approach. This maps the Hubbard-Stratonovic
action to the local impurity problem

$\displaystyle
S_{\mathrm{HS}}^{\mathrm{imp}}=S_{0}^{i}+S_{\phi\phi}^{\mathrm{imp}}+S_{\phi\sigma}^{\mathrm{imp}},$
(14) $\displaystyle
S_{\phi\sigma}^{\mathrm{imp}}={+}\frac{1}{2}\int_{-\infty}^{\infty}dt\phi^{\dagger}_{i}(t)\tau_{3}\left[\varsigma_{i}(t)-h(t)\right]+\text{hc},$
$\displaystyle
S_{\phi\phi}^{\mathrm{imp}}=\frac{1}{2}\int_{-\infty}^{\infty}dtdt^{\prime}\phi^{\dagger}_{i}(t)\tau_{3}\mathcal{W}^{-1}\left(t-t^{\prime}\right)\tau_{3}\phi_{i}\left(t^{\prime}\right).$

where $S_{0}^{i}$ is the action of a single isolated atom on site $i$. The
interaction of the impurity site with the rest of the lattice is mimicked by
the effective field $h(t)$, analogously an in a Gutzwiller mean-field
approach, and by the effective Weiss field $\mathcal{W}$, assumed homogeneous
in space.

From the impurity model, one can compute the local bosonic expectation value
$\left\langle\phi_{i}(t)\right\rangle_{S_{\mathrm{HS}}}$ and the connected
correlation function
$U_{ii}(t-t^{\prime})=-i\left\langle\phi_{i}(t)\phi^{\dagger}_{i}(t^{\prime})\right\rangle_{S_{\mathrm{HS}}^{\mathrm{imp}}}^{\mathrm{con}}$.
This defines the impurity self-energy via the Dyson equation

$\Pi_{\mathrm{loc}}=\mathcal{W}^{-1}-U_{ii}^{-1}.$ (15)

At the same time, assuming translation invariance, the lattice Green’s
function for the bosons
$U_{rr^{\prime}}(t-t^{\prime})=-i\left\langle\phi_{r}(t)\phi^{\dagger}_{r^{\prime}}(t^{\prime})\right\rangle_{S_{\mathrm{HS}}}^{\text{con
}}$ obeys in momentum space the Dyson equation
$U_{k}^{-1}=W_{k}^{-1}-\Pi_{k}$. In DMFT the self energy $\Pi_{k}$ is assumed
to be local in space and is identified with the impurity self-energy
$\Pi_{k}=\Pi_{\text{loc }}$. Thus, summing over momentum, the local
correlation function must also satisfy

$\begin{split}U_{ii}&=\frac{1}{N}\sum_{k}\left[W_{k}^{-1}-\Pi_{\mathrm{loc}}\right]^{-1}\end{split}$
(16)

which from (15) determines the Weiss field $\mathcal{W}$. In addition,
assuming $\left\langle\phi_{r}(t)\right\rangle_{S_{\mathrm{HS}}}=\phi$ for all
$r$ and $t$, the field $h$ is given by

$h_{+(-)}=-\left([W_{k=0}^{-1}(\omega=0)]^{R}-[\mathcal{W}(\omega=0)^{-1}]^{R}\right)\phi_{+(-)},$
(17)

where $h$ and $\phi$ take the same value on the $+$ and $-$ contours, as a
classical field and an average value, and $R$ indicates the retarded component
of the Fourier transform of the correlation function (this is a $2\times 2$
matrix in Nambu space).

Given an effective field $h$ and a Weiss field $\mathcal{W}$, the values of
the local correlation function $U_{ii}$ and average bosonic field $\phi$ are
completely determined by solving the impurity problem, with an appropriate
method. But these values must also satisfy (15) (16) and (17), for the
impurity problem to be representative of the lattice model. These conditions
therefore provide a set of closed equations for $h$ and $\mathcal{W}$, that
must be satisfied by a DMFT solution.

###  D.4 Alternative impurity model

The impurity action (14) defines a spin-boson model, where a two-level system
is coupled to a continuum of bosonic modes, with propagator $\mathcal{W}$. To
solve the impurity model, it is convenient to integrate out the bosonic
fields, leading to an equivalent spin model with action

$S^{\mathrm{imp}}=S_{0}^{i}+S_{\mathrm{int},1}^{\mathrm{imp}}+S_{\mathrm{int},2}^{\mathrm{imp}}$
(18)

with an effective field term

$S_{\mathrm{int},1}^{\mathrm{imp}}=\frac{1}{2}\int_{-\infty}^{\infty}dt\varsigma_{i}^{\dagger}(t)\tau_{3}b+\text{hc},$

and retarded interaction

$S_{\mathrm{int},2}^{\mathrm{imp}}=-\frac{1}{2}\int_{-\infty}^{\infty}dtdt^{\prime}\varsigma_{i}^{\dagger}(t)\tau_{3}{\mathcal{W}}\left(t-t^{\prime}\right)\tau_{3}\varsigma_{i}\left(t^{\prime}\right).$

Defining the spin connected Green’s function as
$\chi_{ii}(t-t^{\prime})=-i\left\langle\varsigma_{i}(t)\varsigma^{\dagger}_{i}(t^{\prime})\right\rangle_{S_{\text{imp
}}^{\text{con }}}$, this is related to the bosonic Green’s function by $\quad
U_{ii}=\mathcal{W}+\mathcal{W}\circ\tau_{3}{\chi}_{ii}\tau_{3}\circ\mathcal{W}$,
where $\circ$ indicates real-time convolutions. Using Eq. (15), this yields
the expression for the self-energy in terms of the spin Green function

$\Pi_{\mathrm{loc}}=\left[\tau_{3}{\chi}^{-1}_{ii}\tau_{3}+\mathcal{W}\right]^{-1}$
(19)

Also, from the gaussian integration, the effective field has expression
$b_{+(-)}=\mathcal{W}^{R}(\omega=0)h_{+(-)}$ and demanding that the action is
stationary around the mean-field order parameter $\phi_{+(-)}$ one obtains
$\left\langle\varsigma_{+(-)}\right\rangle$ as a function of $h_{+(-)}$.
Combining those expressions yields

$b_{+(-)}=\left(\mathcal{W}^{R}(\omega={0})-W_{k=0}^{R}(\omega=0)\right)\left\langle\varsigma_{i,+(-)}\right\rangle,$
(20)

These two equations together with (15) and (16) defines an equivalent set of
DMFT equations to those given in the previous paragraph. These are in terms of
spin quantities which are most natural to compute from our impurity solver.

##  Appendix E Numerical implementation of DMFT

Starting from a guess for $b$ and $\mathcal{W}$, such as the mean field
solution for the first and $\mathcal{W}=1/N\sum_{k}W_{k}$ for the second, the
following steps are iterated until a fixed point is reached:

  * •

Solving the impurity model: given $\mathcal{W}$, $b$, the spin model is solved
using the impurity solver, computing the steady-state density matrix
$\rho_{s}$ and atomic correlation function
$\chi=\begin{pmatrix}\chi^{{++}}&\chi^{+-}\\\
\chi^{-+}&\chi^{--}\end{pmatrix}$ for $t>0$, where the entries are $2\times 2$
matrices in Nambu space. Then $\chi$ at negative times $t<0$ is obtained
assuming the steady-state relation

$\chi_{\alpha\beta}^{ab}(-t)=-[\chi^{ba}_{\bar{\beta}\bar{\alpha}}]^{*}(t)$
(21)

where the conjugate-transpose of both Nambu $a,b$ and Keldysh $\alpha\beta$
indices is taken and the Keldysh indices $\alpha,\beta\in[+,-]$ are negated,
such that $\bar{\alpha}=-\alpha$.

  * •

updating the Weiss field $\mathcal{W}$: in frequency and momentum space, where
inverses of Green’s functions are simply given by a matrix inverse of their
$4\times 4$ Nambu-Keldysh stucture, we evaluate the self-consistent equations
reported in the main text:
$\Pi_{\mathrm{loc}}=\left[\tau_{3}{\chi}^{-1}_{ii}\tau_{3}+\mathcal{W}\right]^{-1}$,
$U_{ii}=\frac{1}{N}\sum_{k}\left[W_{k}^{-1}-\Pi_{\mathrm{loc}}\right]^{-1}$,
$\mathcal{W}^{-1}=\Pi_{\mathrm{loc}}+U_{ii}^{-1}$ , and
$b_{\pm}=\left(\mathcal{W}^{R}(\omega={0})-W_{k=0}^{R}(\omega=0)\right)\langle{\varsigma_{i}}\rangle_{\pm}$,
respectively (19),(16),(15) and (20) in this Supplemental. These give new
values for $b$ and $\mathcal{W}$, that are transformed back into real times to
interate the procedure.

An important point of the update procedure concerns how to define the Green’s
functions at time $t=0$ such as to interface the NCA impurity solver discussed
below and Keldsyh field theory. This will be discussed in Sec. G.

##  Appendix F NCA impurity solver in the Nambu formalism

The real-time dynamics of the impurity problem is computed using an impurity
solver based on a non-crossing approximation for the dynamical map
$\hat{\mathcal{V}}(t)$ that evolves the density matrix
$\rho(t)=\hat{\mathcal{V}}(t)\rho(0)$ [schiroScarlatella2019,
scarlatellaSchiro2021, scarlatellaSchiro2024a]. This approach corresponds to a
first-order self-consistent scheme (equivalent to a non-crossing
approximation, NCA), obtained by expanding in the system-bath coupling
(hybridization expansion). The dynamics is obtained by solving the Dyson
equation for the superoperator $\hat{\mathcal{V}}(t)$

$\partial_{t}\hat{\mathcal{V}}(t)=\hat{\mathcal{H}}_{\rm
eff}\hat{\mathcal{V}}(t)+\int_{0}^{t}dt_{1}\hat{\mathcal{D}}\left(t-t_{1}\right)\hat{\mathcal{V}}\left(t_{1}\right)$
(22)

where $\hat{\mathcal{H}}_{\rm eff}=-i\left[H,\bullet\right]+i\sum_{\alpha\beta
ab}\varsigma_{\alpha a}^{\dagger}[\tau_{3}]_{\alpha a,\beta b}b_{\beta b}$,
with the NCA self-energy given by

$\displaystyle\hat{\mathcal{D}}\left(t\right)=$
$\displaystyle-\frac{i}{2}\sum_{\alpha\beta
ab}\alpha\beta\left[\mathcal{W}_{ba}^{\beta\alpha}\left(t\right)\hat{\varsigma}_{\beta
b}^{\dagger}\hat{\mathcal{V}}\left(t\right)\hat{\varsigma}_{\alpha
a}+\mathcal{W}_{ab}^{\alpha\beta}\left(-t\right)\hat{\varsigma}_{\beta
b}\hat{\mathcal{V}}\left(t\right)\hat{\varsigma}_{\alpha a}^{\dagger}\right]=$
(23) $\displaystyle=\frac{1}{2}\sum_{\alpha\beta
ab}\alpha\beta(-i)\mathcal{W}_{ba}^{\beta\alpha}\left(t\right)\hat{\varsigma}_{\beta
b}^{\dagger}\hat{\mathcal{V}}\left(t\right)\hat{\varsigma}_{\alpha
a}+\text{hc}$

Here we introduced the vectors of spin super-operators
$\hat{\varsigma}=(\sigma^{-}\bullet,\sigma^{+}\bullet,\bullet\sigma^{-},\bullet\sigma^{+})^{T}$
and
$\hat{\varsigma}^{\dagger}=(\sigma^{+}\bullet,\sigma^{-}\bullet,\bullet\sigma^{+},\bullet\sigma^{-})$,
where the bullets indicate whether the corresponding operators apply from the
left or from the right, that are analogous to the Nambu-Keldysh vectors of
fields in the Keldysh formalism. We refer to their components, i.e.
$\hat{\varsigma}_{\alpha a}$, using a Keldysh index $\alpha=+,(-)$ addressing
the first(second) half of the vector and a Nambu index $a=1(2)$ taking the
corresponding first (second) component. An analogous notation is used for the
other symbols. The second equality is obtained by using that the baths
correlation functions at time $-t$ are contained in the Nambu components at
$t$, and are found using the steady-state relation (21).

Note that here all quantities are written as functions of a single time
assuming a time-independent stationary-state. Note also that the $(2,2)$ Nambu
component gives a contribution to the dissipator equal to the $(1,1)$
component, since
$\mathcal{W}_{22}^{\beta\alpha}(t)=\mathcal{W}_{11}^{\alpha\beta}(-t)$.
Analogously, the two contributions given by the $(1,2)$ and $(2,1)$ components
are equal to each other by (21). This justifies the overall $1/2$ in (23).

Correlation functions are computed using the generalized quantum regression
theorem as in [scarlatellaSchiro2021], that in our steady-state case and for
$t\geq 0$ reads

$\displaystyle\chi_{ab}^{\alpha\beta}\left(t\right)=-i\operatorname{tr}\left[\hat{\zeta}_{\alpha
a}\hat{\mathcal{V}}\left(t\right)\hat{\zeta}_{\beta
b}^{\dagger}\rho_{s}\right]+i\text{tr}\left[\hat{\zeta}_{\alpha
a}\hat{\mathcal{V}}\left(t\right)\rho_{s}\right]\text{tr}\left[\hat{\zeta}_{\beta
b}^{\dagger}\rho_{s}\right]$ (24)

while negative times are computed by (21).

Integrating (22), (23) needs using an implicit discretization scheme for
numerical stability. We use here a second-order Runge-Kutta scheme
[scarlatellaSchiro2024a]. Note that this needs to be adapted as to capture
correctly a $~{}\delta(t)$ contribution in the Weiss field $\mathcal{W}$
[scarlatellaSchiro2024a].

We remark that in the limit of independent atoms the impurity problem becomes
equivalent to a single-site of the original lattice problem, described by a
Lindblad master equation. In this limit our NCA approach exactly reduces to a
Lindblad master equation [scarlatellaSchiro2024a], introducing no further
approximations.

##  Appendix G Definition of Green’s functions at t=0

The Green’s functions defined on the Keldysh contour don’t have a uniquely
defined $t\rightarrow 0$ limit: from the definition of time-ordering they
satisfy the relations

$\displaystyle\mathcal{W}^{++}$ $\displaystyle=\mathcal{W}^{-+}$
$\displaystyle\mathcal{W}^{--}$ $\displaystyle=\mathcal{W}^{+-}$
$\displaystyle\rm{for\quad}t>0$ (25) $\displaystyle\mathcal{W}^{++}$
$\displaystyle=\mathcal{W}^{+-}$ $\displaystyle\mathcal{W}^{--}$
$\displaystyle=\mathcal{W}^{-+}$ $\displaystyle\rm{for\quad}t<0$ (26)

from which one sees that the $t\rightarrow 0^{+}$ and $t\rightarrow 0^{-}$
limits of the $++$ and $--$ components are different.

When discretizing time for numerics though, there’s the problem of how to
define the Green functions at $t=0$. This affects the evaluation of the DMFT
equations, that consist of a series of convolutions in time or matrix products
in frequency, where the value strictly at $t=0$ would not be important, unless
there’s $\delta(t)$ contributions (frequency-independent contributions) to the
Green functions. Since we study a Markovian problem, Markovian contributions
are always generated: for example a natural initial condition for the
effective bath $\mathcal{W}$ is taking $\mathcal{W}(\omega)=W_{r,r^{\prime}}$
frequency independent, that encodes the local physical Markovian dissipator
plus the instantaneous dipole-dipole interaction term. In this case, it is
necessary to define appropriately the Green’s functions at $t=0$, such that
the “causality” structure of Keldysh Green’s functions, i.e. the triangular
structure of Green’s functions in the basis of classical and quantum fields
[kamenevKamenev2023], is preserved by the DMFT equations.

In addition, we have to be able to resolve the $t\rightarrow 0^{\pm}$ limits
of the Green functions, which are still a necessary input for the NCA self-
energy evaluated at $t=0$ (23). These cannot be extrapolated from finite
times, but must be deduced from the $t=0$ value of the Green’s functions. For
this reason we resort to a parametrization that allow to reconstruct the
$t\rightarrow 0^{\pm}$ limits from the $t=0$ value and vice-versa.

###  G.1 t=0 contribution to the NCA dissipator

The Markovian component of NCA self-energy (23) is

$\displaystyle\hat{\mathcal{D}}\left(t=0\right)=$
$\displaystyle-\frac{i}{2}\sum_{\alpha\beta
ab}\alpha\beta\left[\mathcal{W}_{ba}^{\beta\alpha}\left(t\rightarrow
0^{+}\right)\hat{\varsigma}_{\beta b}^{\dagger}\hat{\varsigma}_{\alpha
a}+\mathcal{W}_{ab}^{\alpha\beta}\left(t\rightarrow
0^{-}\right)\hat{\varsigma}_{\beta b}\hat{\varsigma}_{\alpha
a}^{\dagger}\right]$ (27)

using $\hat{\mathcal{V}}\left(0\right)=\mathbb{1}$, where the two separate
limits of $\mathcal{W}$ from positive and negative times enter. For the Nambu
components $\mathcal{W}_{12}$ and $\mathcal{W}_{21}$ those limits are equal
(from their definition), while they’re different for $\mathcal{W}_{11}$ and
$\mathcal{W}_{22}$. Since $\mathcal{W}_{22}$ and $\mathcal{W}_{11}$ are
related, we focus on the latter, whose limits can be parametrized as (25)

$\displaystyle\mathcal{W}_{11}(t\rightarrow 0^{+})$
$\displaystyle=\begin{pmatrix}E_{\kappa}-i\frac{\kappa}{2}&E_{f}-i\frac{f}{2}\\\
E_{\kappa}-i\frac{\kappa}{2}&E_{f}-i\frac{f}{2}\end{pmatrix}$
$\displaystyle\mathcal{W}_{11}(t\rightarrow 0^{-})$
$\displaystyle=\begin{pmatrix}-E_{f}-i\frac{f}{2}&-E_{f}-i\frac{f}{2}\\\
-E_{\kappa}-i\frac{\kappa}{2}&-E_{k}-i\frac{\kappa}{2}\end{pmatrix}$ (28)

where the two are related by (21).

We will read out $\kappa,f,E_{\kappa},E_{f}$ from the Green’s functions at
$t=0$, defined in the following, but there only the difference
$E_{\kappa}-E_{f}$ enters there rather than $E_{\kappa}$ and $E_{f}$
separately. Therefore, we rather use the alternative parametrization:

$\displaystyle\mathcal{W}_{11}(t\rightarrow 0^{+})$
$\displaystyle=\begin{pmatrix}E_{\kappa}-E_{f}-i\frac{\kappa}{2}&-i\frac{f}{2}\\\
-i\frac{\kappa}{2}&-i\frac{f}{2}\end{pmatrix}$
$\displaystyle\mathcal{W}_{11}(t\rightarrow 0^{-})$
$\displaystyle=\begin{pmatrix}-i\frac{f}{2}&-i\frac{f}{2}\\\
-i\frac{\kappa}{2}&-E_{k}+E_{f}-i\frac{\kappa}{2}\end{pmatrix}$ (29)

Using the spin anti-commutation relations, one can check that both
parametrizations yields the same $t=0$ value of the $11$ Nambu component of
the self-energy, defined by
$\hat{\mathcal{D}}=\sum_{ab}\hat{\mathcal{D}}_{ab}$, that reads

$\begin{split}\hat{\mathcal{D}}_{11}(t=0)\rho&=\frac{1}{2}\left[\left(\frac{\kappa}{2}+iE_{\kappa}\right)\left(d\rho
d^{\dagger}-d^{\dagger}d\rho\right)+\left(\frac{f}{2}+iE_{f}\right)\left(d^{\dagger}\rho
d-dd^{\dagger}\rho\right)+\text{hc}\right]=\\\
&=\frac{1}{2}\left\\{-iE_{\kappa}\left[d^{\dagger}d,\rho\right]-iE_{f}\left[dd^{\dagger},\rho\right]+\left[\frac{\kappa}{2}\left(d\rho
d^{\dagger}-d^{\dagger}d\rho+\text{hc}\right)+\frac{f}{2}\left(d^{\dagger}\rho
d-dd^{\dagger}\rho+\text{hc}\right)\right]\right\\}\end{split}$ (30)

where we called $d=\sigma^{-}$ to avoid confusion with Keldysh indices. Then
$\mathcal{W}_{22}(t\rightarrow 0^{+})=\mathcal{W}^{T}_{11}(t\rightarrow
0^{-})$ is obtained from the relation
$\mathcal{W}(t)_{22}^{\alpha\beta}=\mathcal{W}(-t)_{11}^{\beta\alpha}$, while
note that $\mathcal{W}_{12(21)}(t\rightarrow 0^{\pm})$ have all entries equal
by definition (by a similar argument as (25)), thus those are trivially given
by the $t=0$ value.

###  G.2 $t=0$ definition of Keldysh Green functions

The proper $t=0$ structure of Keldysh Green’s functions arises naturally when
describing the effect of a Markovian environment [siebererDiehl2016,
siebererDiehl2014, kamenevKamenev2023]. In this case the $11$ Nambu component
of the Keldysh action reads

$S_{11}=-\frac{1}{2}\int_{-\infty}^{\infty}dt\begin{pmatrix}\vdots\end{pmatrix}^{\dagger}\tau_{3}{\mathcal{W}}_{11}(t=0)\tau_{3}\begin{pmatrix}d_{+}\\\
d_{-}\end{pmatrix}$ (31)

with the Green function of the Markovian bath

${\mathcal{W}}_{11}(t=0)=\begin{pmatrix}E_{k}-i\frac{\kappa}{2}{-E_{f}-i\frac{f}{2}}&-if\\\
-i\kappa&{-E_{k}-i\frac{\kappa}{2}}+E_{f}-i\frac{f}{2}\end{pmatrix}$ (32)

Note that this simply corresponds to
$\mathcal{W}_{11}(t=0)=\mathcal{W}_{11}(t\rightarrow
0^{+})+\mathcal{W}_{11}(t\rightarrow 0^{-})$, with the definition of (29).
Note that (32) enforces the “causal” structure of the Keldsyh action,
corresponding to the triangular structure when a rotation to the classical-
quantum basis is performed.

Finally, the $22$ Nambu component of the action is simply related to the $11$
component by ${\mathcal{W}}_{22}(t=0)={\mathcal{W}}_{11}(t=0)^{T}$. Instead
the $12$ and $21$ terms in the action are
${\mathcal{W}}_{12(21)}(t=0)=2{\mathcal{W}}_{12(21)}(t\rightarrow 0^{\pm})$
since the 2 limits are equal and, when translating a Markovian dissipator to
the corresponding Keldysh action, they add up.

###  G.3 t=0 component in DMFT iterations

In practice, in the numerical iterative scheme to solve the DMFT E equations,
the following steps must be added for the $t=0$ component:

$\displaystyle\rightarrow{\mathcal{W}}(0)\quad{\longrightarrow}\quad\mathcal{W}(0^{+}),\mathcal{W}(0^{-})\underset{\rm
NCA\,impurity\,solver}{\longrightarrow}\chi(0^{+})\chi(0^{-})\longrightarrow\quad{\chi}(0)\quad\underset{\textrm{DMFT\;\;equations}}{\rightarrow}$
(33)

Namely, before passing the effective environment ${\mathcal{W}}$ to the NCA
impurity solver, the $t\rightarrow 0^{\pm}$ limits are resolved, extracting
the parameters $E_{\kappa}-E_{f}$, $\kappa$ and $f$ from ${\mathcal{W}}(t=0)$
and using them to determine $\mathcal{W}(0^{-})$ and $\mathcal{W}(0^{+})$ with
(29). Then, the atom correlation function $\chi(0^{\pm})$ is computed by the
NCA solver (24) and $\chi(0)$ is defined by (32), before being inputted in the
DMFT equations.

## References