# Quantum memory effects in atomic ensembles coupled to photonic cavities

Adam Burgess<EMAIL_ADDRESS>Leverhulme Quantum Biology Doctoral
Training Centre, University of Surrey, Guildford, GU2 7XH, United Kingdom
Advanced Technology Institute, University of Surrey, Guildford, GU2 7XH,
United Kingdom  Department of Physics, University of Surrey, Guildford, GU2
7XH, United Kingdom  Marian Florescu<EMAIL_ADDRESS>Advanced
Technology Institute, University of Surrey, Guildford, GU2 7XH, United Kingdom
Department of Physics, University of Surrey, Guildford, GU2 7XH, United
Kingdom

###### Abstract

In this article we explore the dynamics of many-body atomic systems
symmetrically coupled to a single Lorentzian photonic cavity. Our study
reveals interesting dynamical characteristics including non-zero steady
states, superradiant decay, enhanced energy transfer and the ability to
modulate oscillations in the atomic system by tuning environmental degrees of
freedom. We also analyse a configuration consisting of a three-atom chain
embedded in a photonic cavity. Similarly, we find a strong enhancement of the
energy transfer rate between the two ends of the chain and identified specific
initial conditions that lead to significantly reduced dissipation between the
two atoms at the end of the chain. Another configuration of interest consists
of two symmetrical detuned reservoirs with respect to the atomic system. In
the single-atom case, we show that it is possible to enhance the decay rate of
the system by modulating its reservoir detuning, while in the many-atom case,
this results in dynamics akin to the on-resonant cavity. Finally, we examine
the validity of rotating wave approximation through a direct comparison
against the numerically exact hierarchical equations of motion approach. We
find good agreement in the weak coupling regime while in the intermediate
coupling regime, we identify qualitative similarities, but the rotating wave
approximation becomes less reliable. In the moderate coupling regime, we find
deviation of the steady states due to the formation of mixed photon atom
states.

††preprint: APS/123-QED

##  I Introduction

Recently, there has been a sharp increase in interest in using quantum
technologies and their applications to advance present computing techniques.
This has primarily been driven by significant improvements in the fabrication
of quantum technologies and the robustness of quantum computers with larger
numbers of logical qubits, as well as advancements in quantum algorithm
designs that have the potential to undermine current cryptographic techniques
[1, 2, 3, 4, 5]. However, the tendency of quantum systems to decohere, lose
quantum correlations, and dissipate energy into the environment due to
dephasing and relaxation poses a considerable challenge to their use. For the
field of quantum technologies to continue its unprecedented growth in the
future, it is essential to fully explore the ways in which the environment
affects the evolution of quantum systems and identify ways to control this
influence [6]. Additionally, designing new and improved artificial structures
that can process information on shorter spatial and temporal scales has been
greatly aided by our growing knowledge of the dynamics of atomic systems in
dissipative environments [7, 8, 9].

In the general theory of open quantum systems[10], spin systems coupled to
bosonic reservoirs are an archetypal model used to describe two-level atomic
systems interacting with an quantised electromagnetic field. However, it is
exceedingly difficult to unravel the dynamics of the two systems
simultaneously; typically, just the dynamics of the atomic system are of
interest. To eliminate the environmental degrees of freedom, various
approximations are deployed. A standard approximation is the Markov
approximation, in which induced memory effects of the environment are
neglected. This is appropriate for environments that recover rapidly after
interaction with the system [11]. Although this technique has been successful
in analyzing many systems, it fails in adequately capturing quantum-induced
memory effects. For instance, the local density of states for the
electromagnetic field in micro-structured photonic systems, such as photonic
crystals, fluctuates significantly for frequencies close to the photonic band
gap edge mode. This invalidates the Markov approximation [12, 13, 14], and
atomic systems coupled to such a photonic reservoir undergo considerable non-
Markovian dynamics. This leads to strong coupling between the atomic system
and the photonic reservoir generating dressed atomic states which result in
fractional steady state atomic populations, spectral splitting and sub-natural
line-widths of the atomic transitions [15].

Determining the full impact of structured electromagnetic reservoirs on the
dynamics of quantum emitters embedded in them is also crucial as we start to
scale up the technology of these artificial systems. It has recently been
demonstrated[16] that the coherence lifetime of multi-qubit systems may be
increased by properly characterising the non-Markovian noise. Understanding
the control mechanisms provided by structured environments may lead to novel
architectures for more robust quantum technologies with many applications,
including quantum networks [17] and clocks [18]. Another major difficulty
facing quantum technologies, beyond the preservation of entanglement, is the
generation of entanglement. An interesting recent development suggests that it
is possible to utilise the non-Markovianity of the environment to generate
entanglement between initially uncorrelated atomic systems [19, 20].
Understanding the bi-directional flow of information between system and
environment - the essence of non-Markovian systems - makes it possible to
transfer information and, therefore, quantum correlations between systems
using the environment as an intermediary. Photonic cavity systems are
naturally suited for generating such inter-atomic entanglement, as the atoms
couple strongly to the central cavity mode allowing for efficient energy
transfer between atoms, even without direct inter-atomic coupling. Here we
intend to study such systems.

This article is structured as follows: Section II covers the theoretical
formalism of for describing $N$ atomic systems symmetrically coupled to one
another inside a photonic reservoir associated with Lorentzian cavity and
explores the influence of the reservoir characteristics on the induced system
of interest dynamics. In Section III, we explore a configuration of a three-
atom chain coupled via nearest neighbour interactions inside a Lorentzian
cavity and identify target initial conditions that can lead to significantly
reduced dissipation. In Section IV, we consider a many-atom systems coupled to
two symmetrically detuned cavities and investigate the way in which the
structure of the reservoir impacts the decay rate of the atomic system
excitations. Finally, in Section V, we provide a direct comparison of the
results obtained under the RWA approximation against the predictions of the
numerically exact hierarchical equations of motion approach.

##  II Many Atoms in a Photonic Cavity

Figure 1: A schematic of many atoms embedded within a photonic cavity
reservoir.

(a)

(b)

(c)

(d)

Figure 2: (a) A plot of the time evolution of the excited state population of
a single initially excited atom coupled to varying numbers of atomic systems
with $\Lambda=\gamma=\lambda=0.1\omega_{0}$. (b) Enhancement of dissipation
rate with the number of coupled atomic systems. A plot of the time evolution
of the excited state population of a single initially excited atom coupled to
a varying number of atoms $N-1$.
$\Lambda=0,\lambda=0.1\omega_{0},\gamma=\omega_{0}$. (c) Enhancement of energy
transfer rate with the quality of the cavity. A plot of the time evolution of
the excited state population of a single initially excited atom coupled to 4
other atoms for varying spectral width $\gamma$.
$\Lambda=\lambda=0.1\omega_{0}$. (d) Tuneable frequency of oscillations by
modulation of the coupling strength $\lambda$ for five coupled atoms in the
Lorentzian cavity. $\Lambda=0,\gamma=0.1\omega_{0}$.

We begin our analysis by considering $N$-atoms embedded within a photonic
cavity. A schematic architecture of this system is shown in Fig. 1. The atoms
couple to the local electromagnetic field through the interaction of their
transition dipole with the electric field associated with the random
fluctuations in the environment and there is also a direct inter-atomic
coupling generated by dipole-dipole interaction. The Hamiltonian associated
with this model is given by

$\displaystyle H=$
$\displaystyle\sum_{i=1}^{N}\omega_{0}\sigma_{i}^{+}\sigma_{i}^{-}+\sum_{\lambda}\omega_{\lambda}a^{\dagger}_{\lambda}a_{\lambda}+$
$\displaystyle
i\sum_{i,\lambda}\sigma_{i}^{x}g_{\lambda}(a_{\lambda}-a^{\dagger}_{\lambda})+\sum_{i\neq
j}\Lambda\sigma^{+}_{i}\sigma_{j}^{-},$ (1)

where $\sigma_{i}^{+}$and $\sigma_{i}^{-}$, are the excitation and de-
excitation operators for the $i^{\textrm{th}}$ atomic system respectively,
$a_{\lambda}$ and $a^{\dagger}_{\lambda}$ are the bosonic field annihilation
and creation operators, $\omega_{0}$ and $\omega_{\lambda}$ are the atomic
transition and the $\lambda$-boson mode frequencies, $g_{\lambda}$ is the
coupling strength of the atomic system and the $\lambda$-boson mode [21, 22,
23] and $\Lambda$ is the dipole-dipole coupling strength between atoms. For
systems that interact with a local-electromagnetic field that is near resonant
with the systems’ transition energies, it is useful, and effective, to deploy
the rotating wave approximation (RWA). This allows for the neglecting of terms
that do not conserve the number of excitations in the system as they
contribute rapidly oscillating terms into the dynamics and have a vanishing
average. The RWA Hamiltonian is then given by

$\displaystyle H=$
$\displaystyle\sum_{i=1}^{N}\omega_{0}\sigma_{i}^{+}\sigma_{i}^{-}+\sum_{\lambda}\omega_{\lambda}a^{\dagger}_{\lambda}a_{\lambda}+$
$\displaystyle
i\sum_{i,\lambda}g_{\lambda}(a_{\lambda}\sigma^{+}_{i}-a^{\dagger}_{\lambda}\sigma^{-}_{i})+\sum_{i\neq
j}\Lambda\sigma^{+}_{i}\sigma_{j}^{-}$ (2)

This is equivalent to a Dicke model [24] for $N$-atoms with a coupling between
each of the atoms and is an archetypal model used in quantum optics. We note
here that on a practical side, the RWA allows for the generation of analytical
solutions to many quantum optical problems.

We have shown previously [25] that is possible to explore the single-
excitation dynamics of any system governed by the Hamiltonian above by using
the he environment correlation function and the associated characteristic
spectral density. The spectral density simply defines the distribution and
coupling strengths of our electromagnetic field modes,

$J(\omega)=\sum g_{\lambda}^{2}\delta(\omega-\omega_{\lambda}).$ (3)

For a lossy cavity, the spectral density is given by the Lorentzian
distribution.

$J(\omega)=\frac{\lambda\gamma}{\gamma^{2}+(\omega-\omega_{0})^{2}},$ (4)

where $\gamma$ represents the spectral width of the Lorentzian distribution
and, thus, provides a measure of how dissipative the cavity is. $\lambda$
defines the coupling strength of the cavity to the atomic system of interest.

If we consider the single excitation wavefunction given by

$\phi(0)=\sum_{i}^{N}c_{i}(0)\psi_{i}+\sum_{\lambda}c_{\lambda}(0)\psi_{\lambda},$
(5)

where $\psi_{i}=\ket{i}_{A}\ket{0}_{B}$ is the state wherein the
$i^{\textrm{th}}$ atom is in its excited state, and all other atoms and the
reservoir are in their ground states.
$\psi_{\lambda}=\ket{0}_{A}\ket{\lambda}_{B}$ represents the quantum states
with all atoms in their ground state, and the photon reservoir with its
$\lambda$ mode excited. The time evolution of this state is given by

$\phi(t)=c_{0}\psi_{0}+\sum_{i}^{N}c_{i}(t)\psi_{i}+\sum_{\lambda}c_{\lambda}(t)\psi_{\lambda},$
(6)

due to conservation of excitation number. It is convenient for our discussion
to consider the total polarisation of the system,
$c_{+}(t)=\sum_{i}^{N}c_{i}(t)$. By adiabatic cancelling of the reservoir
modes we generate the following equations of motion

$\displaystyle\dot{c}_{i}=$ $\displaystyle-
iJ(c_{+}-c_{i})-\int^{t}_{0}G(t-t_{1})c_{+}(t_{1})dt_{1},$
$\displaystyle\dot{c}_{+}=$ $\displaystyle-
iJ(N-1)c_{+}-N\int^{t}_{0}G(t-t_{1})c_{+}(t_{1})dt_{1}.$ (7)

We have introduced here the memory kernel

$G(t)=\sum_{\lambda}g_{\lambda}^{2}e^{i(\omega_{0}-\omega_{\lambda})t}.$ (8)

For the Lorentzian cavity system described by the spectral density 4, this
takes the form

$G(t)=\lambda e^{-\gamma t}.$ (9)

By utilising the Laplace transform, we can solve for the time evolution of the
$i^{\textrm{th}}$ atoms excited state amplitude $c_{i}$, as derived in
Appendix A, yielding

$\displaystyle c_{i}(t)=(c_{i}(0)-\frac{c_{+}(0)}{N})e^{i\Lambda t}+$ (10)
$\displaystyle\frac{c_{+}(0)}{N}e^{-\mu^{*}t}\left(\cosh(\Gamma
t)+\frac{\mu}{\Gamma}\sinh(\Gamma t)\right)$

where

$\Gamma=\sqrt{\mu^{2}-\lambda N},$ (11)

and

$\mu=\frac{1}{2}\big{(}\gamma-i\Lambda(N-1)\big{)}.$ (12)

We note that the effects of the environment are entirely determined by the
total polarisation of the system $c_{+}$ due to the symmetric nature of the
atomic couplings. In the long-time limit $t\rightarrow\infty$ the populations
become

$P_{i\infty}=\lim_{t\rightarrow\infty}|c_{i}(t)|^{2}=|c_{i}(0)-\frac{c_{+}(0)}{N}|^{2}.$
(13)

Our results shown Fig. 2(a) unveil the influence of size pf the atomic system
ensemble on the dynamics of the atomic populations. Clearly, the steady state
of $(1-1/N)^{2}$ is attained for initially excited atom 1
$c_{1}(0)=c_{+}(0)=1$. We note that due to the factor $\Gamma$ in the $\cosh$
and $\sinh$ functions, for any value of $N>1$, the populations of the atomic
systems will undergo oscillatory evolution, as long as the initial total
polarisation is non-zero ($c_{+}(0)\neq 0$). If we turn off the dipole-dipole
coupling ($\Lambda=0$), the system displays oscillatory behaviour only when
$\gamma<4\lambda N$. This shows that we can generate oscillatory behaviour by
increasing the number of atomic systems $N$ coupled to our reservoir. However,
by reducing the real component of $\Gamma$ we reduce the counter-action to the
exponential decay, increasing the effective decay rate in the system and
leading to superradiant decay for short times and some transient resurgence
due to the non-Markovian effects (see Fig. 2(b)) . We note that the
exponential decay rates scale with the number of atoms in the systems. To see
this more explicitly, we consider the Taylor expansion around $t=0$ up to the
fourth order. Assuming that $c_{i}(0)=\dfrac{1}{\sqrt{N}}$ then

$\displaystyle P_{i}(t)$ $\displaystyle=|c_{i}(t)|^{2}=$ (14)
$\displaystyle\frac{1}{N}\left(1-\lambda Nt^{2}+\frac{2}{3}\lambda\mu
Nt^{3}-\frac{1}{3}\lambda N(\mu^{2}-\lambda N)t^{4}\right).$

However, this superradiance is limited, as the dynamics are still enveloped by
the decay from $e^{-\mu^{*}t}$ as in Fig.2(b), the decay increases with atomic
system number, but this saturates for large values of $N$, it appears to
increase, however, this is due to an introduced oscillation in the system but
the enveloping decay is the same. Interestingly, we also note that as the
spectral width of the cavity $\gamma$ decreases, there is an increase in the
effective transfer rate between the atomic systems as

$\lim_{\gamma\rightarrow
0}\Gamma=i\sqrt{\frac{1}{4}\Lambda^{2}(N-1)^{2}+\lambda N}$ (15)

which implies that the coupling to the cavity mode increases the energy
transfer rate between the atomic systems. This has been shown in Fig. 2(c),
wherein by modulating the spectral width of the reservoir $\gamma$, where is
apparent that the excitation is transferred at a faster rate when $\gamma$
takes smaller values. This suggests that the environment supports energy
transfer and provides an additional channel. Furthermore, we can tune this
rate for small values of $\gamma$ by modulating the coupling strength to the
cavity mode $\lambda$. As shown in Fig. 2(d), we can now infer some control
over the dynamics of our system by modulating the environment.

##  III Atomic Chain in Lorentzian Cavity

(a)

(b)

(c)

(d)

(e)

(f)

Figure 3: Plots of the time evolution of the excited state population of the
three atoms in a chain in a photonic cavity (thick) and free (dashed) for the
first atom (blue), second (orange), third (green). (a) When the first atom is
initially excited. (b) The same as (a) but extended up to late times and
containing the total polarisation $|c_{+}|=|c_{1}+c_{2}+c_{3}|$ (dark orange).
(c) When all atoms start with the same polarisation
$c_{i}(0)=\frac{1}{\sqrt{3}}$. (d) When both the atoms at the ends of the
chain are initialised with an equal magnitude but opposite sign
$c_{1}(0)=-c_{3}(0)=\frac{1}{\sqrt{2}}$. (e) When both the first two atoms of
the chain are initialised with an equal magnitude but opposite sign
$c_{1}(0)=-c_{2}(0)=\frac{1}{\sqrt{2}}$. (f) as in (e) but extended to later
times, with total polarisation $|c_{+}|$ (dark orange).
$\Lambda=0.1\omega_{0}$, $\gamma=0.1\omega_{0}$ and $\lambda=0.1\omega_{0}$.

In the previous section, we considered the idealised system in which we
assumed symmetric all to all coupling between the atomic systems. A system of
this kind are realisable, but some creative engineering of the physical
platform is required, for example, employing coupled Josephson parametric
oscillators and cavity buses[26, 27]. A slightly more practical physical
instantiation would consist of atomic chains embedded in cavity structures,
such that the atomic systems are coupled to their nearest neighbours. Here, we
consider the dynamics of a three-atom chain in a photonic cavity.

For a chain of atomic systems of three atomic systems embedded within a single
photonic cavity, we system evolution is described by the following Hamiltonian

$\displaystyle H=$
$\displaystyle\sum_{i=1}^{3}\omega_{0}\sigma_{i}^{+}\sigma_{i}^{-}+\sum_{\lambda}\omega_{\lambda}a^{\dagger}_{\lambda}a_{\lambda}+$
(16) $\displaystyle
i\sum_{i,\lambda}g_{\lambda}(a_{\lambda}\sigma^{+}_{i}-a^{\dagger}_{\lambda}\sigma^{-}_{i})+\sum^{3}_{i}\Lambda\big{(}\sigma^{+}_{i}\sigma_{i+1}^{-}+\sigma^{+}_{i+1}\sigma_{i}^{-}\big{)}.$

where we are using the same notations as in Eqn. 2, but we note that now the
atoms only interact with their nearest neighbour.

Following the same procedure as before we obtain the following equations of
motion for the excited state amplitudes

$\displaystyle\dot{c}_{1}$ $\displaystyle=-i\Lambda
c_{2}-\int_{0}^{t}G(t-t_{1})c_{+}(t_{1})dt_{1}$ $\displaystyle\dot{c}_{2}$
$\displaystyle=-i\Lambda(c_{+}-c_{2})-\int_{0}^{t}G(t-t_{1})c_{+}(t_{1})dt_{1}$
$\displaystyle\dot{c}_{3}$ $\displaystyle=-i\Lambda
c_{2}-\int_{0}^{t}G(t-t_{1})c_{+}(t_{1})dt_{1}.$ (17)

These coupled differential equations in Eqn. 17 can be solved through Laplace
transformations. Before that, we can identify a few dynamical qualities by
looking at the form of these differential equations. If we turn off the inter-
atomic coupling $\Lambda=0$, we are led back to the previous system studied.
Let’s consider initial conditions where $c_{1}(0)=-c_{3}(0)$ and that
$c_{2}(0)=0$. In this case, there is no time evolution for the system as the
total polarisation of the atoms seen by the bath and by the central atom is
zero. We have demonstrated this for various initial conditions in Fig. 3. In
Fig. 3(a), we consider that the first atom in the chain is initially excited
and note that the rate of energy transfer between the first and third atoms
becomes higher when the cavity is introduced. This is because the cavity acts
as a secondary channel for the excitation to pass between atoms, supplementing
the conventional hopping along the chain. We also note that the system
achieves a lower occupation of the excited state; this is natural as the
cavity is dissipative, so we anticipate the loss in excitation. Furthermore,
as shown Fig. 3(b), at longer times, the effects of the cavity become weaker
and weaker, leading to free-like oscillations between the first and third
atoms in the system. Physically, this can be interpreted as follows: the total
polarisation $c_{+}$, which couples the atoms to the cavity, has decayed away
and the cavity can no longer interact these atomic systems. This a remarkable
effect as it provides direct means to effectively decouple these atomic
systems from their local photonic environment. When we prepare the atomic
systems with the same initial amplitudes ($c_{i}=\frac{1}{\sqrt{3}}$) as shown
in Fig. 3(c), the evolution of the atomic systems is akin to the single atom
case with slight modulations for the second atom as its coupling is different
to the two end atoms. By choosing the two atoms at the end of the chain to
have equal but opposite initial amplitudes
($c_{1}(0)=-c_{3}(0)=\frac{1}{\sqrt{2}}$), as shown in Fig. 3(d), we have no
evolution as predicted. Finally, in Fig. 3(e) and 3(f), the first two atoms in
the chain are prepared with initial amplitudes of equal magnitude but the
opposite sign ($c_{1}(0)=-c_{2}(0)=\frac{1}{\sqrt{2}}$). Here the total
polarisation $c_{+}$ is initially zero; as such, we do not anticipate much
dissipation, though this will grow as the excitation is passed towards the
third atom leading to some dissipation. Interestingly, in Fig. 3(f) we note
that as the system evolves towards later times, the amount of dissipation is,
in fact, incredibly low, and the two end atomic systems act as if they were
nearly free systems. Furthermore, their effective detuning is considerably
lower than in the case of Fig. 3(b), suggesting this would be a more useful
scheme to generating these free dynamics. Clearly, the results presented here
demonstrate that despite a lossy environment, it is still possible maintain
control over the system evolution, and to massively reduce the amount of
dissipation while retaining the energy transfer between atomic systems.

##  IV Double Dissipative cavity

(a)

(b)

(c)

(d)

(e)

(f)

Figure 4: Plots of the time evolution of the excited state population of a
single atom in a (thick) double (dashed) single detuned (black) on resonance
Lorentzian cavity. (a) For a single atomic system. (b) The magnitude of the
difference between the plots in (a) and the resonant system. (c) Coupled with
four atomic systems. (d) The magnitude of the difference between the plots in
(c) and the resonant system. For the two-cavity system
$\lambda=0.05\omega_{0}$, for the single cavities $\lambda=0.1\omega_{0}$,
$\gamma=0.1\omega_{0}$ and $\Lambda=0.1\omega_{0}$. (e),(f) Heatmaps of the
Cardano discriminant showing the region wherein the single atom system coupled
to the double symmetrically detuned cavities undergoes oscillatory dynamics
(light) no oscillations (dark) oscillations occur for
$\lambda=$(e)$0.01\omega$ (f) $0.1\omega_{0}$.

For a pair of cavities detuned from the transition energy of the atomic system
by a given value $\pm\Delta$, the spectral density associated is given by

$J(\omega)=\frac{\lambda\gamma}{(\omega-(\omega_{0}+\Delta))^{2}+\gamma^{2}}+\frac{\lambda\gamma}{(\omega-(\omega_{0}-\Delta))^{2}+\gamma^{2}},$
(18)

which leads to the associated memory kernel

$G(t)=2\lambda e^{-\gamma t}\cos{\Delta t}.$ (19)

In the following, we are analysing the excited state populations obtained by
solving the corresponding equations of motion given in Eqn. 7. In the Fig.
4(a)-(d), we are comparing the case of the two cavities with symmetric
detuning to the on-resonant case $\Delta=0$ and the single off-resonant cavity
with twice the coupling strength. In Fig. 4(a), we consider the time evolution
of a single atomic system coupled to these three different cavity
configurations. The double cavity system, although of the same coupling
strength, leads to more rapid dissipation of the atomic excitation when
compared to the on-resonant case. It further shows the inverse relation to the
single cavity system, wherein the greater dephasing induces more rapid
dissipation. Conversely, in Fig. 4(c), wherein we have five coupled atomic
systems, the double cavity system appears to be much more like the on-resonant
case. This is further exemplified in Fig. 4(d), where we have plotted the
relative difference between these configurations and the resonant cavity case.
We can see that the two-cavity system across time remains much closer to the
resonant case compared to the single on-resonant cavity. This suggests that if
detuned cavities are inevitable in a physical platform, we could engineer this
dephasing to make our system appear to act on resonance.

In Appendix.B we have shown the solution to the differential equation in the
Laplace transform frame. Via the inverse Laplace transform, we can infer that
the rates in our dynamics will be given by the roots of the polynomial given
by Eqn.51. From the Cardano formula, we can identify the regime under which
the system undergoes non-Markovian dynamics. This is done by determining when
these rates are real or complex, leading to Markovian and non-Markovian
dynamics respectively. Complex valued rates will lead to oscillations in the
population dynamics of the single atom. The discriminant is given by

$\Delta_{0}=\left(3\Delta^{2}-\gamma^{2}+6\lambda\right)^{3}+\left(9\Delta^{2}\gamma+\gamma^{3}-9\gamma\lambda\right)^{2},$
(20)

when $\Delta_{0}$ is negative, the system undergoes Markovian dynamics, and
for positive values, the dynamics are non-Markovian. Heatmaps of this result
are shown in Fig. 4(e) and 4(f). Remarkably, we find that as we expand the
non-Markovian regime for greater values of the detuning this enables the
atomic system to regain the excitation from the reservoir By increasing the
coupling strength, we can again increase this region of non-Markovianity.

##  V Hierarchical Equations of Motion

(a)

(b)

(c)

(d)

(e)

(f)

Figure 5: Plots of the time evolution of the excited state population of atoms
in an on-resonance Lorentzian cavity solved using HEOM (thick) and RWA
(dashed/dashed dot). (a) For the first atom initially excited for varying
numbers of atoms in the system, with all to all coupling in the weak coupling
regime $\lambda=0.01\omega_{0}$ with RWA predicted steady states (dashed). (b)
as in (a) in the intermediate coupling regime $\lambda=0.1\omega_{0}$. (c) Two
coupled atoms for varying spectral widths $\gamma$. (d) HEOM (thick) vs RWA
(dashed) predictions for a single atom inside of a detuned cavity with
detuning $\Delta=0.1\omega_{0}$(blue), $0.3\omega_{0}$ (red) and
$0.6\omega_{0}$ (green). (e) HEOM results for Three atom chain system, with
initial conditions $c_{1}(0)=-c_{2}(0)=\frac{1}{\sqrt{2}}$ for the population
of the first (blue) second (orange) third (green) atom in the chain. (f) HEOM
results for a single atom coupled strongly ($\lambda=\omega_{0}$) to the
cavity for a variety of spectral widths $\gamma$. For all plots unless
specified $\lambda=0.01\omega_{0}$, $J=0.1\omega_{0}$, $\gamma=0.1\omega_{0}$.

Throughout this paper, we have deployed the rotating wave approximation. This
is a common approximation within the field of quantum optics that greatly
reduces the complexity of physical problems by neglecting rapidly oscillating
terms in the Hamiltonian, allowing for the restriction of the Hilbert space to
a fixed excitation number. For the systems considered in this article (atomic
systems with equal transition energies and Lorentzian photonic cavities on
resonance), this approximation is appropriate as most frequencies are
approximately equal to the transition energy $\omega_{0}$. Therefore, most
counter-rotating terms will oscillate with frequencies $2\omega_{0}$ and they
can be effectively neglected.

In this section, we explore the accuracy of the RWA and validate our results
by deploying the hierarchical equations of motion (HEOM)[28]. These are a set
of differential equations capable of modelling the dynamics of a quantum
system linearly coupled to a bosonic environment. This is precisely the
configuration considered here. Furthermore, these equations of motion are
numerically exact, making them non-perturbative and non-Markovian, imposing no
limiting approximations on their solutions. The HEOM approach utilises a
discretisation of the environment to generate a numerically efficient approach
to capturing the effects of the collective degrees of freedom of the
environment. This is achieved by the construction of coupled differential
equations that form a hierarchy from repeated time-derivatives of the
influence functional, derived initially by Vernon and Feynman[29]. However,
with this approach come some shortcomings. Firstly, it is a numerical approach
and will not lead to analytical results making the search of the parameter
space less informed and less efficient. Secondly, HEOM requires that the bath
correlation function is of an exponential form. Generally, this is the case
for a limited number of spectral densities, but, fortunately, the Lorentzian
photonic cavity spectral density considered does take this form. The bath
correlation function at zero temperature is given by the Fourier transform of
the spectral density.

$C(t)=\int^{\infty}_{0}d\omega J(\omega)e^{-i\omega t},$ (21)

and the Lorentzian cavity system, this simply yields

$C(t)=\lambda e^{-(\gamma+i\omega_{0})t}.$ (22)

We can then rewrite this as an exponential sum over real and imaginary parts

$\displaystyle C(t)=C_{R}(t)+iC_{I}(t),$ (23) $\displaystyle
C_{R}(t)=\frac{\lambda}{2}(e^{-(\gamma+i\omega_{0})t}+e^{-(\gamma-i\omega_{0})t}=\sum_{k=1}^{2}c_{k}^{R}e^{-\gamma_{k}^{R}},$
(24) $\displaystyle
C_{I}(t)=\frac{\lambda}{2i}(e^{-(\gamma+i\omega_{0})t}-e^{-(\gamma-i\omega_{0})t}=\sum_{k=1}^{2}c_{k}^{I}e^{-\gamma_{k}^{I}}.$
(25)

In the non-rotating wave approximation, we can divide the total Hamiltonian
into the system Hamiltonian, containing only atomic degrees of freedom and the
interaction Hamiltonian, as follows

$\displaystyle H=H_{S}+H_{I},$ (26)

with the system Hamiltonian

$H_{S}=\sum_{i=1}^{N}\omega_{0}\sigma_{i}^{+}\sigma_{i}^{-}+\sum_{i<j}\Lambda\sigma^{x}_{i}\sigma_{j}^{x},$
(27)

and the interaction Hamiltonian

$\displaystyle
H_{I}=i\sum_{i,\lambda}\sigma_{i}^{x}g_{\lambda}(a_{\lambda}-a^{\dagger}_{\lambda})$
(28)
$\displaystyle=i\sigma_{T}^{x}\sum_{\lambda}g_{\lambda}(a_{\lambda}-a^{\dagger}_{\lambda}),$
(29)

and $\sigma_{T}^{x}=\sum_{i=1}^{N}\sigma_{i}^{x}$.

The hierarchical equations of motion are then[28]

$\displaystyle\dot{\rho}^{n}(t)$
$\displaystyle=\left(-iH^{\times}_{S}-\sum_{j=R,I}\sum_{k=1}^{2}n_{jk}\gamma^{j}_{k}\right)\rho^{n}(t)$
$\displaystyle-i\sum^{2}_{k=1}c_{k}^{R}n_{Rk}\sigma_{x}^{T\times}\rho^{n^{-}_{Rk}}(t)+\sum^{2}_{k=1}c_{k}^{I}n_{Ik}\sigma_{x}^{To}\rho^{n^{-}_{Ik}}(t)$
$\displaystyle-i\sum_{j=R,I}\sum_{k=1}^{2}\sigma_{x}^{T\times}\rho^{n^{+}_{jk}}(t),$
(30)

where we used the multi-index
$n=(n_{R1},n_{R2},n_{I1},n_{I2}),n_{R,Ii}\in\\{0,...,N_{c}\\}$ with $N_{c}$
the cut-off parameter defining the depth of the hierarchy needed for a
converged result. We have also introduced the commutator and anticommutator
notation $\sigma_{x}^{T\times}=[\sigma_{x}^{T},\cdot]$ and
$\sigma_{x}^{To}=\\{\sigma_{x}^{T},\cdot\\}$. In Eqn. 30, the only physical
density matrix is the $(0,0,0,0)$ indexed matrix which is the reduced atomic
systems’ density matrix. All other $\rho^{n}$ are auxiliary and encapsulate
the environmental attributes. The density operators with index $n^{\pm}_{jk}$
refer to auxiliary density operators with the $k$th index of the real or
imaginary term raised or lowered by one. For example, for index

$\displaystyle n$ $\displaystyle=(1,2,1,3)$ (31) $\displaystyle n^{+}_{R2}$
$\displaystyle=(1,3,1,3).$

In Fig. 5, we compare the RWA solutions to the numerically exact HEOM
solutions. We note that the RWA begins to break down as we increase the
coupling strength of the system, limiting the utility of this approximation to
weak coupling regime. In Fig. 5(a) we consider a coupling strength of
$\lambda=0.01\omega_{0}$ and we note a good agreement between both the HEOM
and RWA predictions. However, at later times, a phase difference between HEOM
and RWA emerges as expected since the counter-rotating terms become more
relevant. This has little effect on the steady state characteristics and the
RWA and HEOM predict nearly identical steady states for the system. This is
quite remarkable as in the RWA we restrict our dynamics to a greatly reduced
subspace of the entire Hilbert space. The total Hamiltonian does not conserve
excitation number and the dynamics can traverse larger domains of the Hilbert
space, yet we still have the RWA predicted steady state. Increasing the
coupling strength in Fig. 5(b), we note that the RWA approximation starts to
break down, and the agreement between the RWA and the HEOM approaches is only
qualitative. This is due to the discrepancy between the RWA Hamiltonians and
the total Hamiltonian, the system’s ground state is not equivalent. In RWA,
the ground state is associated with zero excitations in the entire system. In
the total Hamiltonian description, we have the ground state associated with
mixed states of the reservoir degrees of freedom and the atomic systems. The
degree of mixing is proportional to the coupling strength of the atomic
systems to the reservoir $\lambda$. Therefore, as we increase the coupling
strength, we have greater mixing in the total system steady state, and we see
a more significant divergence between the RWA and HEOM steady states.

Furthermore, as the coupling strength increases, we increase the number of
non-excitation number conserving interactions, leading to further deviation.
In Fig. 5(c), we have varied the spectral width of the cavity and compared the
HEOM and RWA solutions. Here we have two conflicting effects, as we decrease
spectral width, we increase the coupling strength to the central cavity mode,
making the RWA less reliable. However, by reducing the width, we also reduce
the effect of non-resonant modes in the cavity, which should improve the
accuracy of RWA. As such, we note that the intermediate value of the spectral
width $\gamma=0.1\omega_{0}$ provides the best agreement. In Fig. 5(d),
somewhat surprisingly, we find that in the weak coupling regime, the HEOM and
the RWA approximation results coincide very well for detuned cavities, even
across large values of the detuning, suggesting that we can learn about
detuned systems within the RWA approximation. Remarkably, we are able to
validate the phenomenon of reduced dissipation within the atomic chain which
leads to almost free dynamics for the two ends of the chain in Fig. 5(e),
confirming that is not an artefact of the RWA.

Due to the non-perturbative nature of HEOM we can also explore the strong
coupling regime. In Fig.5(f) we consider the dynamics of a single atom
strongly coupled ($\lambda=\omega_{0}$) to the cavity for a few parameters of
the spectral width ($\gamma$). Here we can see additional phenomena emerging.
Firstly, we can see that the for the case wherein $\gamma=\omega_{0}$ we have
exponential decay in HEOM, however, within the RWA, we expect to see
oscillatory dynamics. This shows that the strong coupling limit begins to
modulate the condition on oscillatory dynamics away from the standard
$\gamma^{2}<4\lambda$ relation. Furthermore, we note that for the narrow
spectral width $\gamma=0.1\omega_{0}$, we have an introduction of at least one
additional phase into the oscillatory dynamics moving away from a single phase
enveloped by a decay. Interestingly, we do see that the RWA predictions does
seem to coincide with the peaks of the HEOM calculations, showing that some of
the RWA dynamics remain, the RWA predictions also give good qualitative
agreement of the decay rates with HEOM. Finally, we find that the largest
spectral width $\gamma=5\omega_{0}$ has the greatest discrepancy from the
HEOM, this is particularly startling as the effective coupling strength here
is lower for higher spectral width, yet we see that this configuration has the
highest overall steady state.

As we contrast the RWA and HEOM approaches, we note that in the weak and
intermediate coupling regime, the RWA analytical formalism provides valuable
insights into the actual physical dynamics of the system and can guide our
search using numerically exact techniques. On the other hand, HEOM provides us
with a powerful tool to study parameter regimes inaccessible in RWA.

##  VI Conclusion

To conclude, we have explored the dynamics of many-body atomic systems
symmetrically coupled to a single Lorentzian photonic cavity and identified a
number of intriguing dynamical traits. This includes superradiant decay, non-
zero stable states, improved energy transfer, and the capacity to control
oscillations in the atomic system by adjusting environmental degrees of
freedom. We also explored multi-atomic systems coupled via nearest neighbour
interaction to each other and to the photonic reservoir of modes of the
cavity. In this case, we predict a similar increase in the rate of energy
transfer between the two ends of the chain. We have also identified the
initial conditions which result in a significant reduction of and exchange
excitations between the two atoms at the end of the chain. We further
demonstrated the ability of the structured photonic reservoir to influence the
atomic evolution in the case of two symmetrical detuned reservoirs. Here, for
the single-atom case, we show that is possible to increase the system’s decay
rate by modulating its detuning with respect to the reservoir spectral
features. In the many-atom case, we show that reservoir structuring brings the
system dynamics’ closer to the on-resonant cavity case. Finally, we
demonstrated the RWA predictions by deploying a numerically exact HEOM
approach and extended our results in the regime in which the RWA is no longer
valid. Our results demonstrate that despite the dissipative nature of the
photonic reservoirs, careful choices of initial conditions or engineering of
the environment can provide a robust control over the embedded atomic systems
dynamics.

###### Acknowledgements.

This work was supported by the Leverhulme Quantum Biology Doctoral Training
Centre at the University of Surrey were funded by a Leverhulme Trust training
centre grant number DS-2017-079, and the EPSRC (United Kingdom) Strategic
Equipment Grant No. EP/L02263X/1 (EP/M008576/1) and EPSRC (United Kingdom)
Grant EP/M027791/1 awards to M.F. We acknowledge helpful discussions with the
members of the Leverhulme Quantum Biology Doctoral Training Centre.

##  Appendix A Single Reservoir

###  A Fully Symmetric Coupling

We start with the (RWA) Hamiltonian

$\displaystyle H=$
$\displaystyle\sum_{i=1}^{N}\omega_{0}\sigma_{i}^{+}\sigma_{i}^{-}+\sum_{\lambda}\omega_{\lambda}a^{\dagger}_{\lambda}a_{\lambda}$
$\displaystyle+i\sum_{i,\lambda}g_{\lambda}(a_{\lambda}\sigma^{+}_{i}-a^{\dagger}_{\lambda}\sigma^{-}_{i})+\sum_{i\neq
j}\Lambda\sigma^{+}_{i}\sigma_{j}^{-},$ (32)

where $\sigma_{i}^{+}$and $\sigma_{i}^{-}$, are the excitation and de-
excitation operators for the $i^{\textrm{th}}$ atomic system respectively,
$a_{\lambda}$ and $a^{\dagger}_{\lambda}$ are the bosonic field annihilation
and creation operators, $\omega_{0}$ and $\omega_{\lambda}$ are the atomic
transition and the $\lambda$-boson mode frequencies, and $g_{\lambda}$ is the
coupling strength of the atomic system and the $\lambda$-boson mode.

By changing to the interacting picture, we can consider the interaction
Hamiltonian

$\displaystyle\tilde{H}_{I}=$ $\displaystyle
i\sum_{i,\lambda}g_{\lambda}(a_{\lambda}\sigma^{+}_{i}e^{i(\omega_{0}-\omega_{\lambda})t}-a^{\dagger}_{\lambda}\sigma^{-}_{i}e^{-i(\omega_{0}-\omega_{\lambda})t})$
$\displaystyle+\sum_{j\neq k}\Lambda\sigma^{+}_{j}\sigma_{k}^{-}.$ (33)

Such a Hamiltonian is convenient as it conserves the excitation number of any
wavefunction. As such, if we consider the single excitation wavefunction given
by

$\phi(0)=c_{0}\psi_{0}+\sum_{i}^{N}c_{i}(0)\psi_{i}+\sum_{\lambda}c_{\lambda}(0)\psi_{\lambda},$
(34)

its time evolution is given by

$\phi(t)=c_{0}\psi_{0}+\sum_{i}^{N}c_{i}(t)\psi_{i}+\sum_{\lambda}c_{\lambda}(t)\psi_{\lambda},$
(35)

where $\psi_{i}=\ket{i}_{A}\ket{0}_{B}$ the ith atom is in its excited state
and $\psi_{\lambda}=\ket{0}_{A}\ket{\lambda}_{B}$ the atoms are all in their
ground state, and the bosonic system has it is $\lambda$ mode excited. By
simply plugging this into the Schrödinger equation and determining the coupled
differential equations of motion for the state amplitudes, we get

$\tilde{H}_{I}\phi_{i}=-i\sum_{\lambda}g_{\lambda}\psi_{\lambda}e^{-i(\omega_{0}-\omega_{\lambda})t}+\sum_{j}^{N}\Lambda\psi_{j}(1-\delta_{ij}),$
(36)

$\tilde{H}_{I}\phi_{\lambda}=ig_{\lambda}e^{i(\omega_{0}-\omega_{\lambda})t}\sum_{i}\psi_{i}.$
(37)

By introducing the parameter

$c_{+}(t)=\sum_{i}c_{i}(t),$ (38)

we have

$\displaystyle\dot{c}_{i}=$
$\displaystyle-i\Lambda(c_{+}-c_{i})+\sum_{\lambda}c_{\lambda}g_{\lambda}e^{i(\omega_{0}-\omega_{\lambda})t},$
$\displaystyle\dot{c}_{+}=$
$\displaystyle-i\Lambda(N-1)c_{+}+N\sum_{\lambda}c_{\lambda}g_{\lambda}e^{i(\omega_{0}-\omega_{\lambda})t},$
$\displaystyle\dot{c}_{\lambda}=$ $\displaystyle-
g_{\lambda}e^{-i(\omega_{0}-\omega_{\lambda})t}c_{+}.$ (39)

Assuming the bosonic field is initially in its vacuum configuration
($c_{\lambda}(0)=0$), we can formally integrate up the last equation to get

$c_{\lambda}(t)=-\int^{t}_{0}g_{\lambda}e^{-(\omega_{0}-\omega_{\lambda})t_{1}}c_{+}(t_{1})dt_{1},$
(40)

and by introducing the memory kernel
$G(t)=\sum_{\lambda}g_{\lambda}^{2}e^{i(\omega_{0}-\omega_{\lambda})t}$ we can
rewrite the above equations as

$\displaystyle\dot{c}_{i}=$
$\displaystyle-i\Lambda(c_{+}-c_{i})-\int^{t}_{0}G(t-t_{1})c_{+}(t_{1})dt_{1},$
$\displaystyle\dot{c}_{+}=$
$\displaystyle-i\Lambda(N-1)c_{+}-N\int^{t}_{0}G(t-t_{1})c_{+}(t_{1})dt_{1}.$
(41)

It is notable that our memory kernel $G(t)$ is related to the spectral density
$S(\omega)=\sum_{\lambda}g_{\lambda}^{2}\delta(\omega-\omega_{\lambda})$
associated to the reservoir by the relation

$G(t)=\int d\omega S(\omega)e^{i(\omega_{0}-\omega)t}.$ (42)

In order to solve these equations, we utilise the Laplace transform.
Transforming the differential equations into algebraic ones.

$\displaystyle s\tilde{c}_{i}(s)-c_{i}(0)$
$\displaystyle=-i\Lambda(\tilde{c}_{+}(s)-\tilde{c}_{i}(s))-\tilde{c}_{+}(s)\tilde{G}_{1},$
$\displaystyle s\tilde{c}_{+}(s)-c_{+}(0)$
$\displaystyle=-i\Lambda(N-1)\tilde{c}_{+}(s)-N\tilde{c}_{+}(s)\tilde{G}_{1}.$
(43)

Where we have used that $\tilde{f}(s)=\mathcal{L}\\{f(t)\\}$. We can see that
we can solve for $\tilde{c}_{+}(s)$ getting

$\tilde{c}_{+}(s)=\frac{c_{+}(0)}{s+i\Lambda(N-1)+N\tilde{G}}.$ (44)

From this, we can solve for the single atom amplitude Laplace solution

$\tilde{c}_{i}(s)=\frac{c_{i}(0)}{s-i\Lambda}-\frac{c_{+}(0)(\tilde{G}+i\Lambda)}{(s-i\Lambda)(s+i\Lambda(N-1)+N\tilde{G})}.$
(45)

If we consider the form of $\tilde{G}(s)$ with respect to the spectral
density, we can see that

$\tilde{G}(s)=\int d\omega\frac{S(\omega)}{s-i(\omega_{0}-\omega)}.$ (46)

##  Appendix B Double Cavity System

The Laplace transform of the memory kernel is given by

$\tilde{G}(s)=\frac{2\lambda(\gamma+s)}{\Delta^{2}+(\gamma+s)^{2}}.$ (47)

In Laplace space the solutions to the equations of motion are given by

$\displaystyle c_{i}(s)=\frac{c_{i}(0)}{s-i\Lambda}-$ (48)
$\displaystyle\frac{c_{+}(0)\left(i\Lambda(\Delta^{2}+(\gamma+s)^{2})+2\lambda(\gamma+s)\right)}{(s-i\Lambda)\big{(}(i\Lambda(N-1)+s)(\Delta^{2}+(\gamma+s)^{2})+2\lambda
N(\gamma+s)\big{)}}.$ (49)

Considering just the single atom case $N=1$, $J=0$. We have that

$\displaystyle
c_{i}(s)=\frac{c_{i}(0)(\Delta^{2}+(\gamma+s)^{2})}{s(\Delta^{2}+(\gamma+s)^{2})+2\lambda(\gamma+s)}.$
(50)

The roots of the denominator will determine the rates of the equation by the
inverse Laplace transform. The expanded denominator will give

$\lambda\gamma+(\Delta^{2}+\gamma^{2}+2\lambda)s+2\gamma s^{2}+s^{3}.$ (51)

## References

  * Wang  _et al._ [2018] Y. Wang, Y. Li, Z.-q. Yin, and B. Zeng, npj Quantum Information 4, 46 (2018).
  * Van Meter and Itoh [2005] R. Van Meter and K. M. Itoh, Phys. Rev. A 71, 052320 (2005).
  * Markov and Saeedi [2013] I. L. Markov and M. Saeedi, Phys. Rev. A 87, 012310 (2013).
  * Saeedi and Pedram [2013] M. Saeedi and M. Pedram, Phys. Rev. A 87, 062318 (2013).
  * Shor [1994] P. Shor, in  _Proceedings 35th Annual Symposium on Foundations of Computer Science_ (1994) pp. 124–134.
  * Dong and Petersen [2010] D. Dong and I. Petersen, IET Control Theory & Applications 4, 2651 (2010).
  * Petta  _et al._ [2005] J. R. Petta, A. C. Johnson, J. M. Taylor, E. A. Laird, A. Yacoby, M. D. Lukin, C. M. Marcus, M. P. Hanson, and A. C. Gossard, Science 309, 2180 (2005).
  * Yamamoto  _et al._ [2003] T. Yamamoto, Y. A. Pashkin, O. Astafiev, Y. Nakamura, and J. S. Tsai, Nature 425, 941 (2003).
  * McDermott  _et al._ [2005] R. McDermott, R. W. Simmonds, M. Steffen, K. B. Cooper, K. Cicak, K. D. Osborn, S. Oh, D. P. Pappas, and J. M. Martinis, Science 307, 1299 (2005).
  * Breuer and Petruccione [2002] H. P. Breuer and F. Petruccione,  _The theory of open quantum systems_ (Oxford University Press, Great Clarendon Street, 2002).
  * Moy  _et al._ [1999] G. M. Moy, J. J. Hope, and C. M. Savage, Phys. Rev. A 59, 667 (1999).
  * Burgess and Florescu [2021] A. Burgess and M. Florescu, Opt. Mat. Express 11, 2037 (2021).
  * Florescu and John [2001] M. Florescu and S. John, Phys. Rev. A 64, 033801 (2001).
  * John and Florescu [2001] S. John and M. Florescu, Journal of Optics A: Pure and Applied Optics 3, 10.1088/1464-4258/3/6/361 (2001).
  * John and Quang [1994] S. John and T. Quang, Physical Review A 50, 1764 (1994).
  * White  _et al._ [2020] G. A. L. White, C. D. Hill, F. A. Pollock, L. C. L. Hollenberg, and K. Modi, Nature Communications 11, 6301 (2020).
  * Kimble [2008] H. J. Kimble, Nature 453, 1023 (2008).
  * Kómár  _et al._ [2014] P. Kómár, E. M. Kessler, M. Bishof, L. Jiang, A. S. Sørensen, J. Ye, and M. D. Lukin, Nature Physics 10, 582 (2014).
  * Mirkin  _et al._ [2019] N. Mirkin, P. Poggi, and D. Wisniacki, Phys. Rev. A 99, 062327 (2019).
  * Fleming  _et al._ [2012] C. H. Fleming, N. I. Cummings, C. Anastopoulos, and B. L. Hu, Journal of Physics A: Mathematical and Theoretical 45, 065301 (2012).
  * Pfeifer [1982] P. Pfeifer, Phys. Rev. A 26, 701 (1982).
  * Allen and Eberly [1987] L. Allen and J. Eberly,  _Optical Resonance and Two-level Atoms_, Dover books on physics and chemistry (Dover, 1987).
  * John and Quang [1995] S. John and T. Quang, Phys. Rev. Lett. 74, 3419 (1995).
  * Dicke [1954] R. H. Dicke, Phys. Rev. 93, 99 (1954).
  * Burgess and Florescu [2022] A. Burgess and M. Florescu, Phys. Rev. A 105, 062207 (2022).
  * Majer  _et al._ [2007] J. Majer, J. M. Chow, J. M. Gambetta, J. Koch, B. R. Johnson, J. A. Schreier, L. Frunzio, D. I. Schuster, A. A. Houck, A. Wallraff, A. Blais, M. H. Devoret, S. M. Girvin, and R. J. Schoelkopf, Nature 449, 443 (2007).
  * Onodera  _et al._ [2020] T. Onodera, E. Ng, and P. L. McMahon, npj Quantum Information 6, 48 (2020).
  * Tanimura [2020] Y. Tanimura, The Journal of Chemical Physics 153, 020901 (2020), https://doi.org/10.1063/5.0011599 .
  * Feynman and Vernon [1963] R. Feynman and F. Vernon, Annals of Physics 24, 118 (1963).


# Fixes That Fail: Self-Defeating Improvements in Machine-Learning Systems

Ruihan Wu
Cornell University
<EMAIL_ADDRESS>
&Chuan Guo Awni Hannun Laurens van der Maaten
Facebook AI Research
<EMAIL_ADDRESS>
Work performed during internship at Facebook AI Research.

###### Abstract

Machine-learning systems such as self-driving cars or virtual assistants are
composed of a large number of machine-learning models that recognize image
content, transcribe speech, analyze natural language, infer preferences, rank
options, _etc_. Models in these systems are often developed and trained
independently, which raises an obvious concern: _Can improving a machine-
learning model make the overall system worse?_ We answer this question
affirmatively by showing that improving a model can deteriorate the
performance of downstream models, even after those downstream models are
retrained. Such _self-defeating improvements_ are the result of entanglement
between the models in the system. We perform an error decomposition of systems
with multiple machine-learning models, which sheds light on the types of
errors that can lead to self-defeating improvements. We also present the
results of experiments which show that self-defeating improvements emerge in a
realistic stereo-based detection system for cars and pedestrians.

##  1 Introduction

Progress in machine learning has allowed us to develop increasingly
sophisticated artificially intelligent systems, including self-driving cars,
virtual assistants, and complex robots. These systems generally contain a
large number of machine-learning models that are trained to perform modular
tasks such as recognizing image or video content, transcribing speech,
analyzing natural language, eliciting user preferences, ranking options to be
presented to a user, _etc_. The models feed each other information: for
example, a pedestrian detector may use the output of a depth-estimation model
as input. Indeed, machine-learning systems can be interpreted as directed
acyclic graphs in which each vertex corresponds to a model, and models feed
each other information over the edges in the graph.

In practice, various models in this graph are often developed and trained
independently [21]. This modularization tends to lead to more usable APIs but
may also be motivated by other practical constraints. For example, some of the
models in the system may be developed by different teams (some models in the
system may be developed by a cloud provider); models may be retrained and
deployed at a very different cadence (personalization models may be updated
very regularly but image-recognition models may not); new downstream models
may be added continuously in the graph (user-specific models may need to be
created every time a user signs up for a service); or models may be non-
differentiable (gradient-boosted decision trees are commonly used for
selection of discrete features). This can make backpropagation through the
models in the graph infeasible or highly undesirable in many practical
settings.

The scenario described above leads to an obvious potential concern: _Can
improving a machine-learning model make the machine-learning system worse?_ We
answer this question affirmatively by showing how improvements in an upstream
model can deteriorate the performance of downstream models, _even if all the
downstream models are retrained_ after updating the upstream model. Such
_self-defeating improvements_ are caused by entanglement between the models in
the system. To better understand this phenomenon, we perform an error
decomposition of simple machine-learning systems. Our decomposition sheds
light on the different types of errors that can lead to self-defeating
improvements: we provide illustrations of each error type. We also show that
self-defeating improvements arise in a realistic system that detects cars and
pedestrians in stereo images.

Our study opens up a plethora of new research questions that, to the best of
our knowledge, have not yet been studied in-depth in the machine-learning
community. We hope that our study will encourage the community to move beyond
the study of machine-learning models in isolation, and to study machine-
learning systems more holistically instead.

##  2 Problem Setting

We model a machine-learning _system_ as a static directed acyclic graph (DAG)
of machine-learning _models_ , $G=(\mathcal{V},\mathcal{E})$, where each
vertex $v\in\mathcal{V}$ corresponds to a machine-learning model and each edge
$(v,w)\in\mathcal{E}$ represents the output of model $v$ being used as input
into model $w$. For example, vertex $v$ could be a depth-estimation model and
vertex $w$ a pedestrian-detection model that operates on the depth estimates
produced by vertex $v$. A model is _upstream_ of $v$ if it is an ancestor of
$v$. A _downstream_ model of $v$ is a descendant of $v$. An example of a
machine-learning system with three models is shown in Figure 1 (we study this
model in Section 3.2).

Figure 1: Example of a machine-learning system with three models,
$\mathcal{V}=\\{u,v,w\\}$, and two edges, $\mathcal{E}=\\{(u,w),(v,w)\\}$. The
system receives an example $\mathbf{x}$ as input, of which parts
$\mathbf{x}^{(u)}$, $\mathbf{x}^{(v)}$, and $\mathbf{x}^{(w)}$ are fed into
$u$, $v$, and $w$, respectively. Model $w$’s input is concatenated with the
outputs of its parents, $f_{u}(\mathbf{x})$ and $f_{v}(\mathbf{x})$.

Each vertex in the machine-learning system $G$ corresponds to a model
function, $v(\cdot)$, that operates on both data and outputs from its upstream
models. Denote the output from model $v$ by $f_{v}(\mathbf{x})$. The output of
a source model $v$ is $f_{v}(\mathbf{x})=v(\mathbf{x}^{(v)})$, where
$\mathbf{x}^{(v)}$ represents the part of input $\mathbf{x}$ that is relevant
to model $v$. The output of a non-source model $w$ is given by
$f_{w}(\mathbf{x})=w\left(\left[\mathbf{x}^{(w)},f_{v}(\mathbf{x}):v\in\mathcal{P}_{w}\right]\right)$,
where $\mathcal{P}_{w}=\\{v:(v,w)\in\mathcal{E}\\}$ are the parents of $w$.

Training. To train model $v$, we assume access to a training set
$\mathcal{D}_{v}=\\{(\mathbf{x}_{1},\mathbf{y}_{1}),\dots,(\mathbf{x}_{N_{v}},\mathbf{y}_{N_{v}})\\}$
with $N_{v}$ examples $\mathbf{x}_{n}\in\mathcal{X}$ and corresponding targets
$\mathbf{y}_{n}\in\mathcal{Y}_{v}$. For each model, we also assume a training
algorithm $A_{v}(\mathcal{D}_{v})$ that selects a model $v\in\mathcal{H}_{v}$
from hypothesis set $\mathcal{H}_{v}$ based on the training set. The learning
algorithm, $A_{v}$, the hypothesis set, $\mathcal{H}_{v}$, and the training
set, $\mathcal{D}_{v}$, are fixed during training but they may change each
time model $v$ is re-trained or updated. The data space $\mathcal{X}$, the
target space $\mathcal{Y}_{v}$, and the way in which inputs $\mathbf{x}^{(v)}$
are obtained do not change between model updates.

We assume that models in $G$ may be trained _separately_ rather than jointly,
that is, the learning signal obtained from training a downstream model, $v$,
may not be backpropagated into its upstream dependencies,
$w\in\mathcal{P}_{v}$. This assumption is in line with constraints that
commonly arise in real-world machine learning systems; see our motivation
below. Note that this assumption also implies that a downstream model cannot
be trained before all its upstream models are trained.

Evaluation. To be able to evaluate the performance of the task at $v$, we
assume access to a _fixed_ test set $\bar{\mathcal{D}}_{v}$ that is defined
analogously to the training set: it contains $M_{v}$ test examples
$\bar{\mathbf{x}}_{m}\in\mathcal{X}$ and corresponding targets
$\bar{\mathbf{y}}_{m}\in\mathcal{Y}_{v}$. In contrast to the training set, the
test set is fixed and does not change when the model is updated. However, we
note that the input distribution into a non-root model is not only governed by
the test set $\bar{\mathcal{D}}_{v}$ but also by upstream models: if an
upstream model changes, the inputs into a downstream model may change as well
even though the test set is fixed.

We also assume access to a _test loss_ function
$\ell_{v}(v,\bar{\mathcal{D}}_{v})$ for each model $v$, for example,
classification error (we assume lower is better). Akin to the test set
$\bar{\mathcal{D}}_{v}$, all test loss functions $\ell_{v}$ are fixed. Hence,
we always measure the generalization ability of a model, $v$, in the exact
same way. The test loss function may not have trainable parameters, that is,
it must be a “true” loss function.

Updating models. The machine-learning system we defined can be improved by
updating a model. In such an update, an individual model $v$ is replaced by an
alternative model $v^{\prime}$. We can determine if such a model update
constitutes an _improvement_ by evaluating whether or not its test loss
decreases:
$\ell_{v}(v^{\prime},\bar{\mathcal{D}}_{v})\leq\ell_{v}(v,\bar{\mathcal{D}}_{v})$.
When $v$ is not a leaf in $G$, the update to $v^{\prime}$ may affect the test
loss of downstream models as well. In practice, the improvement due to
$v^{\prime}$ often also improves downstream models (_e.g._ , [18]), but this
is not guaranteed. We define a _self-defeating improvement_ as a situation in
which replacing $v$ by $v^{\prime}$ improves that model but deteriorates at
least one of the downstream models $w$, even after all the downstream models
are updated (_i.e._ , re-trained using the same learning algorithm) to account
for the changes incurred by $v^{\prime}$.

######  Definition 2.1 (Self-defeating improvement).

Denoting the set of all downstream models (descendants) of $v$ in $G$ as
$\mathcal{C}_{v}$, a _self-defeating improvement_ arises when:

$\exists
w\in\mathcal{C}_{v}:\ell_{v}(v^{\prime},\bar{\mathcal{D}}_{v})\leq\ell_{v}(v,\bar{\mathcal{D}}_{v})\centernot\implies\ell_{w}(w^{\prime},\bar{\mathcal{D}}_{w})\leq\ell_{w}(w,\bar{\mathcal{D}}_{w})),$
(1)

where $w^{\prime}$ represents a new version of model $w$ that was trained
after model $v$ was updated to $v^{\prime}$ to account for the input
distribution change that the update of $v$ incurs on $w$.

Motivation of problem setting. In the problem setting described above, we have
made two assumptions: (1) models may be (re-)trained separately rather than
jointly and (2) the training algorithm, hypothesis set, and training data may
change between updates of the model111The training data is fixed during
training, _i.e._ , we do not consider online learning settings in our
setting.. Both assumptions are motivated by the practice of developing large
machine-learning systems that may comprise thousands of models [21, 27]. Joint
model training is often infeasible in such systems because:

  * •

Some of the models may have been developed by cloud providers and cannot be
changed.

  * •

Backpropagation may be technically infeasible or inefficient, _e.g._ , when
models are implemented in incompatible learning frameworks or when their
training data live in different data centers.

  * •

Some models may be re-trained or updated more often than others. For example,
personalization models may be updated every few hours to adapt to changes in
the data distribution, but re-training large language models [8] at that same
cadence is not very practical.

  * •

Upstream models may have hundreds or even thousands of downstream
dependencies, complicating appropriate weighting of the downstream learning
signals that multi-task learning [10] requires.

  * •

Some models may be non-differentiable, _e.g._ , gradient-boosted decision
trees are commonly used for feature selection in large sets of features.

We assume that changes in the training algorithm, hypothesis set, and training
data may occur because model owners (for example, cloud providers) constantly
seek to train and deploy improved versions of their models with the ultimate
goal of improving the system(s) in which the models are used.

##  3 Understanding Self-Defeating Improvements via Error Decompositions

To ensure that improvements to a part benefit the whole, traditional software
design relies heavily on _modules_ : independent and interchangeable
components that are combined to form complex software systems [7, §7.4]. Well-
designed modules can be modified independently without affecting the behavior
of other modules in the system. As a result, improvements to a module (_e.g._
, to its efficiency) are guaranteed to improve the entire system. Ideally, the
machine-learning models in the system from Section 2 would behave akin to
traditional software modules. If they did, improving any machine-learning
model would improve the entire system.

However, machine-learned models rarely behave like modules in practice because
they lack an explicit specification of their functionality [13]. As a result,
there is often a significant degree of _entanglement_ between different models
in the system [21]. This can lead to situations where improving a model does
not improve the entire system. To understand such self-defeating improvements,
we study Bayes error decompositions of simple machine-learning systems.

Figure 2: Illustration of a self-defeating improvement due to an increase in
the downstream approximation error. Upstream model $v$ predicts a 2D point for
each example. Downstream model $w$ is a linear classifier separating the
points into two classes (red and green). Left: Ground-truth target
distribution for the upstream model. Middle: Point predictions by upstream
model $v$, the optimal $v$-conditioned downstream classifier $w^{\dagger}$,
and the downstream classifier $w$ learned. Right: Point predictions by
improved upstream model $v^{\prime}$, the optimal $v^{\prime}$-conditioned
downstream classifier $w^{\dagger}$, and the downstream classifier
$w^{\prime}$ learned. Note how the predictions of $v^{\prime}$ better match
the ground-truth, but $w^{\prime}$ cannot exploit the improvement because it
is restricted to be linear.

###  3.1 Error Decomposition of System with Two Models

We adapt standard Bayes error decomposition [6, 32] to study a simple system
that has two vertices, $\mathcal{V}=\\{v,w\\}$, and a dependency between
upstream model $v$ and downstream model $w$, that is,
$\mathcal{E}=\\{(v,w)\\}$. We denote the Bayes-optimal downstream model by
$w^{*}$; the optimal downstream model conditioned on an upstream model by
$w^{\dagger}$; and the optimal upstream-conditional downstream model in the
hypothesis set $\mathcal{H}_{w}$ by $w^{\dagger}_{\mathcal{H}_{w}}$. Using
these definitions, we can decompose the _downstream risk_ of model $w$, given
upstream model $v$, as follows:

$\displaystyle\mathbb{E}[\ell_{w}(w\circ v)-\ell_{w}(w^{*})]=$
$\displaystyle\underbrace{\mathbb{E}[\ell_{w}(w^{\dagger}\circ
v)-\ell_{w}(w^{*})]}_{\text{upstream
error}}+\underbrace{\mathbb{E}[\ell_{w}(w^{\dagger}_{\mathcal{H}_{w}}\circ
v)-\ell_{w}(w^{\dagger}\circ v)]}_{\text{downstream approximation
error}}+\underbrace{\mathbb{E}[\ell_{w}(w\circ
v)-\ell_{w}(w^{\dagger}_{\mathcal{H}_{w}}\circ v)]}_{\text{downstream
estimation error}},$ (2)

where $\circ$ denotes function composition, and the expectations are under the
data distribution. The error decomposition is similar to standard
decompositions [6, 32] but contains three terms instead of two.222For brevity,
we do not include optimization error of [5] in our Bayes error decomposition.
Specifically, the _upstream error_ does not arise in prior error
decompositions, and the _downstream approximation_ and _downstream estimation_
errors differ from the standard decomposition.

Barring variance in the error estimate due to the finite size of test sets
$\bar{\mathcal{D}}_{w}$ (which is an issue that can arise in any model-
selection problem), a self-defeating improvement occurs when a model update
from $(v,w)$ to $(v^{\prime},w^{\prime})$ reduces the upstream risk but
increases the downstream risk. An increase in the downstream risk implies that
at least one of the three errors terms in the composition above must have
increased. We describe how each of these terms may increase due to a model
update:

  * •

Upstream error measures the error that is due to the upstream model $v$ not
being part of the Bayes-optimal solution $w^{*}$. It increases when an
improvement in the upstream loss does not translate in a reduction in Bayes
error of the downstream model, _i.e._ , when: $\ell_{w}(w^{\dagger}\circ
v^{\prime})>\ell_{w}(w^{\dagger}\circ v)$. An upstream error increase can
happen due the _loss mismatch_ : a setting in which the test loss functions,
$\ell$, of the upstream and downstream model optimize for different things.
For example, the upstream test loss function may not penalize errors that need
to be avoided in order for the downstream test loss to be minimized, or it may
penalize errors that are irrelevant to the downstream model.

Alternatively, the upstream error may increase due to _distribution mismatch_
: situations in which the upstream loss focuses on parts of the data-target
space $\mathcal{X}\times\mathcal{Y}_{v}$ that are unimportant for the
downstream model (and vice versa). For example, suppose the upstream model $v$
is an image-recognition model that aims to learn a feature representation that
separates cats from dogs, and $f_{v}(\mathbf{x})$ is a feature representation
obtained from that model. A downstream model, $w$, that distinguishes
different dog breeds based on $f_{v}$ may deteriorate when the improvement of
model $v$ collapses all representations of dog images in $f_{v}$. Examples of
this were observed in, _e.g._ , [18].

  * •

Downstream approximation error measures the error due to the optimal
$w^{\dagger}$ not being in the hypothesis set $\mathcal{H}_{w}$. The
approximation error increases when the downstream model, $w$, is unable to
exploit improvements in the upstream model, $v$, because exploiting those
improvements would require selection of a model that is not in that hypothesis
set.

Figure 2 illustrates how this situation can lead to a self-defeating
improvement. In the illustration, the upstream model, $v$, predicts a 2D
position for each data point. The downstream model, $w$, separates examples
into the positive (green color) and negative (red color) class based on the
prediction of $v$ using a linear classifier. The predictions of the improved
upstream model, $v^{\prime}$, better match the ground-truth targets: the
predictions of $v^{\prime}$ (right plot) match the ground-truth targets (left
plot) more closely than those of $v$ (middle plot). However, the re-trained
downstream model, $w^{\prime}$, deteriorates because the resulting
classification problem is more non-linear: the downstream model cannot
capitalize on the upstream improvement because it is linear. The resulting
increase in approximation error leads to the self-defeating improvement in
Figure 2.

Figure 3: Illustration of a self-defeating improvement due to an increase in
the downstream estimation error. Upstream model $v$ predicts a 2D point for
each example. Downstream model $w$ is a quadratic classifier separating the
points into two classes (red and green). Left: Ground-truth target samples and
distribution for the upstream model. Middle: Point predictions by upstream
model $v$, the corresponding optimal $v$-conditioned downstream model in
$\mathcal{H}_{w}$, $w^{\dagger}_{\mathcal{H}_{w}}$, and the downstream
classifier $w$ learned based on the training examples. Right: Point
predictions by improved upstream model $v^{\prime}$, the corresponding optimal
$v^{\prime}$-conditioned downstream model in $\mathcal{H}_{w}$,
$w^{\dagger}_{\mathcal{H}_{w}}$, and the downstream classifier $w^{\prime}$
learned based on the training examples. Note how the predictions of
$v^{\prime}$ better match the ground-truth, but $w^{\prime}$ deteriorates
because it only receives $N_{w}=6$ training examples.

  * •

Downstream estimation error measures the error due to the model training being
performed on a finite data sample with an imperfect optimizer, which makes
finding $w^{\dagger}_{\mathcal{H}_{w}}$ difficult in practice. The estimation
increases, for example, when the upstream model improves but the downstream
model requires more than $N_{w}$ training samples to capitalize on this
improvement.

Figure 3 shows an example of a self-defeating improvement caused by an
increase of the downstream estimation error. As before, the upstream model,
$v$, predicts a 2D position for each data point in the example. The downstream
model $w$ is selected from the set, $\mathcal{H}_{w}$, of quadratic
classifiers. It is tasked with performing binary classification into a
positive class (green color) and negative class (red color) based on the 2D
positions predicted by the upstream model. To train the binary classifier, the
downstream model $w$ (and $w^{\prime}$) is provided with $N_{w}=6$ labeled
training examples (the green and red markers in the plot). In the example, the
upstream model $v^{\prime}$ performs better than its original counterpart $v$:
the predictions of $v^{\prime}$ (right plot) match the ground-truth targets
(left plot) more closely than those of $v$ (middle plot). However, the
upstream improvement hampers the downstream model even though the optimal
downstream model $w^{\dagger}$ is in the hypothesis set $\mathcal{H}_{w}$ both
before and after the upstream improvement. After the upstream improvement, the
optimal downstream model that can be selected from the hypothesis set,
$w^{\dagger}_{\mathcal{H}_{w}}$, is more complex, which makes it harder to
find it based on the finite number of $N_{w}$ training examples. This results
in an increase in the estimation error, which causes the self-defeating
improvement in Figure 3.

Figure 4: Illustrations of self-defeating improvement due to an increase in
the upstream compatibility error. Green areas should be labeled as positive;
red areas as negative. Left (a and b): Original upstream model $u$ (top; blue
line and $+,-$ symbols) and its improved version $u^{\prime}$ (bottom; red).
Middle (a and b): Upstream model $v$. Right (a and b): Original downstream
model $w$ (top) and its improved version $w^{\prime}$ (bottom). In (a), the
errors of $u$ and $v$ are anti-correlated. As a result, the improved upstream
model $u^{\prime}$ negatively impacts the re-trained downstream model
$w^{\prime}$. In (b), the improved upstream model $u^{\prime}$ is identical to
the other upstream model $v$. As a result, $w^{\prime}$ looses the additional
information previously provided by $u$, negatively impacting its performance.

###  3.2 Error Decomposition of System with Two Upstream Models

More complex types of self-defeating improvements may arise when a downstream
model depends on multiple upstream models that are themselves entangled.
Consider the system of Figure 1 that has three models,
$\mathcal{V}=\\{u,v,w\\}$, and dependencies between upstream models $u$ and
$v$ and downstream model $w$, that is, $\mathcal{E}=\\{(u,w),(v,w)\\}$. The
error decomposition for this system is similar to Equation 2, but we can
further decompose the upstream error. Denoting the Bayes-optimal downstream
model by $w^{*}$, the optimal model given upstream models $u$ and $v$ by
$w^{\dagger}_{u,v}$, and the optimal upstream model $v$ given upstream model
$u$ by $v^{\dagger}_{u}$, we decompose the upstream error as follows:

$\displaystyle\underbrace{\mathbb{E}[\ell_{w}(w^{\dagger}_{u,v}\circ(u,v))-\ell_{w}(w^{*})]}_{\text{upstream
error}}=$
$\displaystyle\underbrace{\mathbb{E}[\ell_{w}(w^{\dagger}_{u,v}\circ(u,v))-\ell_{w}(w^{\dagger}_{u,v}\circ(u,v^{\dagger}_{u}))]}_{\text{upstream
compatibility
error}}+\underbrace{\mathbb{E}[\ell_{w}(w^{\dagger}_{u,v}\circ(u,v^{\dagger}_{u}))-\ell_{w}(w^{*})]}_{\text{excess
upstream error}}.$ (3)

The excess upstream error is similar to the upstream error in Equation 2:
upstream model $u$ may be suboptimal for the downstream task, for example,
because of loss mismatch or distribution mismatch. The key observation in the
error decomposition of a system with two upstream models is that the optimal
upstream model $v$ is a function of upstream model $u$ (and vice versa). The
upstream compatibility error captures the error due to upstream model $v$ not
being identical to the optimal $v^{\dagger}_{u}$. A self-defeating improvement
can occur when we update upstream model $u$ to $u^{\prime}$ because it may be
that $v^{\dagger}_{u}\neq v^{\dagger}_{u^{\prime}}$, which can cause the
upstream compatibility error to increase.

We provide two examples of this in Figure 4. The first example (left pane)
shows a self-defeating improvement due to upstream models $u$ and $v$ making
_anti-correlated errors_ that cancel each other out. The second example (right
pane) is a self-defeating improvement due to $u^{\prime}$ making more
_correlated predictions_ with $v$. We note that, in both examples, the optimal
$v$ depends on $u$ and the self-defeating improvement arises because of an
increase in the upstream compatibility error.

In the examples in Figure 4, all three models aim to distinguish two classes:
green (positive class) and red (negative class). Upstream models $u$ and $v$
do so based the $(x,y)$-location of points in the 2D plane. The downstream
model operates on the (hard) outputs of the two upstream models. In Figure
4(a), the original upstream models $u$ and $v$ make errors that are anti-
correlated. The original downstream model $w$ exploits this anti-correlation
to make perfect predictions. When upstream model $u$ is replaced by an
improved model $u^{\prime}$ that makes no errors, a self-defeating improvement
arises: downstream model $w$ no longer receives the information it needs to
separate both classes perfectly. In Figure 4(b), upstream model $u$ is
improved to $u^{\prime}$, which makes the exact same predictions as the other
upstream model $v$. As a result, downstream model $w^{\prime}$ no longer has
access to the complementary information that $u$ was providing in addition to
$v$, producing the self-defeating improvement.

While the examples in Figure 4 may seem contrived, anti-correlated errors and
correlated predictions can arise in practice when there are dependencies in
the targets used to train different models. For example, suppose the upstream
models are trained to predict ad clicks (_will a user click on an ad?_) and ad
conversions (_will a user purchase the product being advertised?_),
respectively. Because conversion can only happen after a click happens, there
are strong dependencies between the targets used to train both models, which
may introduce correlations between their predictions.

Figure 5: Overview of the car and pedestrian detection system used in our case
study in Section 4. The system consists of three models: (1) a depth-
estimation model $u$ that generates a pseudo-LiDAR representation of the
environment, (2) a car-detection model $v$ that operates on 3D point clouds,
and (3) a pedestrian-detection model $w$ that performs 3D pedestrian
detection.

##  4 Case Study: Pseudo-LiDAR for Detection

We perform a case study on self-defeating improvements in a system that may be
found in self-driving cars. In both cases, the system uses a pseudo-LiDAR [35]
as the basis for 3D detection of cars and pedestrians. Figure 5 gives an
overview of the models in the system:

  * •

A disparity-estimation model, $u$, predicts a map $\mathbf{d}$ with
disparities corresponding to every pixel in the left image of a stereo image
(a pair of rectified images from two horizontally aligned cameras).

  * •

A point cloud generator, $r$, uses disparity map $\mathbf{d}$ to create a
pseudo-LiDAR representation of the 3D environment. The point cloud generator,
$r$, is a simple non-learnable module [35].

  * •

The stereo image and disparity map are input into a model, $v$, that performs
3D detection of cars.

  * •

The same inputs are used in a separate model, $w$, that aims to perform 3D
pedestrian detection.

|  | Car detection | Pedestrian detection
---|---|---|---
|  | P-RCNN | F-PTNET | P-RCNN | F-PTNET
AP${}_{\text{3D}}$ | $\boldsymbol{u}$ | 39.91 | 33.64 | 34.42 | 43.02
$\boldsymbol{u^{\prime}}$ | 38.75 | 32.78 | 27.64 | 43.10
$\boldsymbol{y}$ | 68.06 | 55.46 | 55.19 | 62.16
AP${}_{\text{BEV}}$ | $\boldsymbol{u}$ | 50.34 | 43.79 | 38.51 | 50.49
$\boldsymbol{u^{\prime}}$ | 49.30 | 42.77 | 31.47 | 48.29
$\boldsymbol{y}$ | 78.43 | 67.35 | 63.00 | 65.92

Table 1: Average precision of car detection and pedestrian detection measured
for 3D box view (AP${}_{\text{3D}}$) and 2D birds-eye view
(AP${}_{\text{BEV}}$). Higher is better. Results are shown for baseline
disparity-estimation model $u$ and an improved disparity-estimation model
$u^{\prime}$ that is trained using a different loss. The AP of oracle
disparity-estimation model $\boldsymbol{y}$ is shown for reference. Improving
the disparity-estimation model leads to a self-defeating improvement in both
the car detection model and the pedestrian detection model.

###  4.1 System Design

The disparity-estimation and detection models used in our case study are
described below.

Disparity estimation. To perform disparity estimation, we adopt the PSMNet
model of Chang and Chen [11]. The model receives a stereo image,
$(\mathbf{x}_{\textrm{left}},\mathbf{x}_{\textrm{right}})$, as input and aims
to compute a disparity map,
$\mathbf{d}=u(\mathbf{x}_{\textrm{left}},\mathbf{x}_{\textrm{right}})$, that
contains a disparity estimate for each pixel in image
$\mathbf{x}_{\textrm{left}}$. We experiment with two training loss functions:
(1) the depth mean absolute error and (2) the disparity mean absolute error
(MAE). Given a ground-truth disparity map, $\mathbf{y}$, and a function that
maps disparities to depths, $g(\mathbf{d})=\frac{C}{\mathbf{d}}$, for some
camera constant $C$, the disparity MAE is proportional to
$\lVert\mathbf{d}-\mathbf{y}\rVert_{1}$ and the depth MAE is proportional to
$\lVert g(\mathbf{d})-g(\mathbf{y})\rVert_{1}$. To evaluate the quality of
model $u$, we measure disparity MAE:
$\ell_{u}(u,\bar{\mathcal{D}}_{u})=\frac{1}{M}\sum_{m=1}^{M}\lVert\mathbf{d}-\mathbf{y}\rVert_{1}$.
We expect that a model, $u^{\prime}$, trained to minimize disparity MAE will
provide better disparity estimates (per test loss $\ell_{u}$) than a model,
$u$, that is trained to minimize depth MAE. This has downstream effects on
performance of the detection models.

Car and pedestrian detection. We experiment with two different models for
performing detection of cars and pedestrians, _viz._ , the Point-RCNN model
(P-RCNN; Shi et al. [29]) and the Frustum PointNet detection model (F-PTNET;
Qi et al. [25]). Both models take stereo image
$(\mathbf{x}_{\textrm{left}},\mathbf{x}_{\textrm{right}})$ and point cloud
$r(\mathbf{d})$ as input. Before computing $r(\mathbf{d})$, we perform
winsorization on the prediction
$u(\mathbf{x}_{\textrm{left}},\mathbf{x}_{\textrm{right}})$: we remove all
points that are higher than $1$ meter in the point cloud (where the camera
position is the origin). Both 3D car detection model and 3D pedestrian
detection model are trained to detect their target objects at any distance
from the camera.

Following Geiger et al. [15], we evaluate the test loss of car-detection
model, $\ell_{v}$, using the negative average precision (AP) for an
intersection-over-union (IoU) of $0.7$. The test loss of the pedestrian-
detection model, $\ell_{w}$, is the negative AP for an IoU of $0.5$. The
pedestrian-detection test-loss is only evaluated on pedestrians whose distance
to the camera is $\leq\\!20$ meters. We measure both the car-detection and the
pedestrian-detection test losses for both the 3D box view
($-\text{AP}_{\text{3D}}$) and the 2D bird-eye view
($-\text{AP}_{\text{BEV}}$); see Geiger et al. [15] for details on the
definition of the test-loss functions.

Range | 0-max | 0-20 | 20-40 | 40-max
---|---|---|---|---
Model $\boldsymbol{u}$ | 1.28 | 1.32 | 1.12 | 1.11
Model $\boldsymbol{u^{\prime}}$ | 1.21 | 1.21 | 1.20 | 1.20

Table 2: Disparity MAE test loss of disparity-estimation models $u$ (trained
to minimize depth MAE) and $u^{\prime}$ (trained to minimize disparity MAE),
split out per range of ground-truth depth values (in meters). Model $u$ better
predicts the disparity of nearby points, but $u^{\prime}$ works better on
distant points.

###  4.2 Experiments

We evaluate our system on the KITTI dataset [15, CC BY-NC-SA 3.0], using the
training-validation split of [12]. Our baseline disparity-estimation model,
$u$, that is trained to minimize depth MAE obtains a test loss of 1.28 (see
Table 2). The improved version of that model, $u^{\prime}$, is trained to
minimize disparity MAE and obtains a test loss of 1.21, confirming the model
improvement.

|  | P-RCNN | F-PTNET
---|---|---|---
AP${}_{\text{3D}}$ | $\boldsymbol{u}$ | 30.53 | 44.10
$\boldsymbol{u^{\prime}}$ | 34.20 | 46.14
$\boldsymbol{y}$ | 55.19 | 62.16
AP${}_{\text{BEV}}$ | $\boldsymbol{u}$ | 34.09 | 51.84
$\boldsymbol{u^{\prime}}$ | 38.83 | 53.99
$\boldsymbol{y}$ | 64.25 | 65.23

Table 3: Average precision at IoU 0.5 of pedestrian detection measured for 3D
box view (AP${}_{\text{3D}}$) and 2D birds-eye view (AP${}_{\text{BEV}}$) for
two different detection models (P-RCNN and F-PTNET). Pedestrian detector is
trained by removing pedestrians farther than $20$ meters. Results are shown
for the baseline disparity estimation model, $u$, the improved model,
$u^{\prime}$, and an oracle model, $\boldsymbol{y}$. Higher is better.

Table 1 presents the test losses of the downstream models for car and
pedestrian detection, for both the baseline upstream model $u$ and the
improved model $u^{\prime}$. We observe that the improvement in disparity
estimation leads to self-defeating improvements on both the car-detection and
the pedestrian-detection tasks. The AP${}_{\text{3D}}$ of the car detector
drops from $39.91$ to $38.75$ (from $50.34$ to $49.30$ in
AP${}_{\text{BEV}}$). For pedestrian detection, AP${}_{\text{3D}}$ is
unchanged but AP${}_{\text{BEV}}$ drops from $50.49$ to $48.29$.
Interestingly, these self-defeating improvements are due to increases in
different error terms.

Self-defeating improvement in car detection. The self-defeating improvement
observed for car detection is likely due to an increase in the upstream error.
A disparity-estimation model that is trained to minimize depth MAE (_i.e._ ,
model $u$) focuses on making sure that _small_ disparity values (large depth
values), are predicted correctly. By contrast, a model trained to minimize
disparity MAE (model $u^{\prime}$) aims to predict _large_ disparity values
(small depth values) correctly, sacrificing accuracy in predictions of large
depth values (see Table 2). This negatively affects the car detector because
most cars tend to be relatively far away. In other words, the loss function
that improves model $u$ deteriorates the downstream model because it focuses
on errors that are less relevant to that model, which increases the upstream
error.

Self-defeating improvement in pedestrian detection. Different from car
detection, the self-defeating improvement in pedestrian detection is likely
due to an increase in the downstream approximation or estimation error. Table
2 shows that whereas the disparity MAE test loss of $u^{\prime}$ increases for
large depth values, it decreases in the 0-20 meter range. As the pedestrian-
detection evaluation focuses on pedestrians within 20 meters, we expect the
pedestrian detector $w^{\prime}$ to improve along with disparity-estimation
model $u^{\prime}$: that is, the upstream error likely decreases.

However, the point representations produced by model $u^{\prime}$ are likely
noisier for far-away pedestrians than those produced by $u$ (see Table 2).
Because the downstream model is trained to (also) detect these far-away
pedestrians, the noise in the pedestrian point clouds makes it harder for that
model to learn a good model of the shape of pedestrians. This negatively
affects the ability of $w^{\prime}$ to detect nearby pedestrians. Indeed, if
the pedestrian detector is only trained on nearby pedestrians, the self-
defeating improvement disappears; see Table 3. Hence, the observed self-
defeating improvement was due to an increase in a downstream error term: the
upstream model did improve for the downstream task, but the downstream model
was unable to capitalize on that improvement.

##  5 Related Work

This study is part of a larger body of work around _machine-learning
engineering_ [9, 20, 22]. Sculley et al. [27] provides a good overview of the
problems that machine-learning engineering studies. The self-defeating
improvements studied in this paper are due to the entanglement problem [1,
21], where “changing anything changes everything” [27]. Entanglement has also
been studied in work on “backward-compatible” learners [28, 31, 33]. Prior
work has also studied how the performance of humans operating a machine-
learning system may deteriorate when the system is improved [3].

Some aspects of self-defeating improvements have been studied in other
domains. In particular, the effect of an upstream model change on a downstream
model can be viewed as _domain shift_ [4, 24]. However, approaches to domain
shift such as importance weighting of training examples [30] are not directly
applicable because they require the input space of both models to be
identical. The study of self-defeating improvements is also related to the
study of _transfer learning_ , _e.g._ , [2, 18, 23]. Prior work on transfer
learning has reported examples of self-defeating improvements: _e.g._ , [18]
reports that some regularization techniques improve the accuracy of
convolutional networks on ImageNet, but negatively affect the transfer of
those networks to other recognition tasks.

Finally, work on _differentiable programming_ [17, 34] is relevant to the
study of self-defeating improvements. In some situations, a self-defeating
improvement may be resolved by backpropagating learning signals from a model
further upstream. However, such an approach is not a panacea as downstream
tasks may have conflicting objectives. Moreover, technical or practical
obstacles may prevent the use of differentiable programming in practice (see
Section 2).

##  6 Conclusion and Future Work

This study explored self-defeating improvements in modular machine-learning
systems. We presented a new error decomposition that sheds light on the error
sources that can give rise to self-defeating improvements. The findings
presented in this study suggest a plethora of directions for future work:

  * •

It may be possible to derive bounds on the error terms in Equation 2 that go
beyond the current art. For example, it is well-known that the trade-off
between approximation and estimation error leads to an excess risk that scales
between the inverse and the inverse square root of the number of training
examples [26, 37]. New theory may provide more insight into the effect of the
upstream hypothesis set and training set size on those errors, and may bound
the upstream error term in 2.

  * •

Our study focused on _first-order_ entanglement between two models (either an
upstream and a downstream model or two upstream models). Some settings may
give rise to higher-order entanglements, for example, between three upstream
models whose errors cancel each other out.

  * •

Our definition of different types of entanglement may pave the way for the
development of _diagnostic tools_ that identify the root cause of a self-
defeating improvement.

  * •

We need to develop _best practices_ that can help prevent self-defeating
improvements from occurring. One way to do so may be by tightening the “API”
of the machine-learning models, _e.g._ , by calibrating classifier outputs
[16] or whitening feature representations [19]. Such transforms constrain a
model’s output distribution, which may reduce the downstream negative effects
of model changes. Alternatively, we may train models that are inherently more
robust to changes in their input distribution, _e.g._ , using positive-
congruent training [36] or meta-learning [14].

  * •

We also need to study self-defeating improvements of other aspects. Indeed, an
important open question is: _Can making an upstream model fairer make its
downstream dependencies less fair?_

We hope this paper will encourage the study of such questions, and will
inspire the machine-learning community to rigorously study machine-learning
systems in addition to individual models.

## Acknowledgements

The authors thank Yan Wang, Yurong You, and Ari Morcos for helpful discussions
and advice.

## References

  * Amershi et al. [2019] S. Amershi, A. Begel, C. Bird, R. DeLine, H. Gall, E. Kamar, N. Nagappan, B. Nushi, and T. Zimmermann.  Software engineering for machine learning: A case study.  In _Proceedings of the International Conference on Software Engineering_ , 2019.
  * Azizpour et al. [2016] H. Azizpour, A. S. Razavian, J. Sullivan, A. Maki, and S. Carlsson.  Factors of transferability for a generic convnet representation.  _IEEE Transactions on Pattern Analysis and Machince Intelligence_ , 38(9):1790–1802, 2016.
  * Bansal et al. [2020] G. Bansal, B. Nushi, E. Kamar, D. S. Weld, W. S. Lasecki, and E. Horvitz.  Updates in human-AI teams: Understanding and addressing the performance/compatibility tradeoff.  In _Proceedings of the AAAI Conference on Artificial Intelligence_ , volume 33, pages 2429–2437, 2020.
  * Ben-David et al. [2010] S. Ben-David, J. Blitzer, K. Crammer, A. Kulesza, F. Pereira, and J. W. Vaughan.  A theory of learning from different domains.  _Machine Learning_ , 79:151–175, 2010.
  * Bottou and Bousquet [2008] L. Bottou and O. Bousquet.  The tradeoffs of large-scale learning.  In _Advances in Neural Information Processing Systems (NeurIPS)_ , pages 161–168, 2008.
  * Boucheron et al. [2005] S. Boucheron, O. Bousquet, and G. Lugosi.  Theory of classification: a survey of recent advances.  _ESAIM: Probability and Statistics_ , 9:323–375, 2005.
  * Brookshear [2009] J. Brookshear.  _Computer Science. An Overview ( $9^{th}$ edition)_.  Addison-Wesley, 2009.
  * Brown et al. [2020] T. B. Brown, B. Mann, N. Ryder, M. Subbiah, J. D. Kaplan, P. Dhariwal, A. Neelakantan, P. Shyam, G. Sastry, A. Askell, S. Agarwal, A. Herbert-Voss, G. M. Krueger, T. Henighan, R. Child, A. Ramesh, D. Ziegler, J. Wu, C. Winter, C. Hesse, M. Chen, E. Sigler, M. Litwin, S. Gray, B. Chess, J. Clark, C. Berner, S. McCandlish, A. Radford, I. Sutskever, and D. Amodei.  Language models are few-shot learners.  In _Advances in Neural Information Processing_ , 2020.
  * Burkov [2020] A. Burkov.  _Machine Learning Engineering_.  True Positive Inc., 2020.
  * Caruana [1997] R. Caruana.  Multi-task learning.  _Machine Learning_ , 28:41–75, 1997.
  * Chang and Chen [2018] J.-R. Chang and Y.-S. Chen.  Pyramid stereo matching network.  In _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition_ , pages 5410–5418, 2018.
  * Chen et al. [2015] X. Chen, K. Kundu, Y. Zhu, A. G. Berneshawi, H. Ma, S. Fidler, and R. Urtasun.  3D object proposals for accurate object class detection.  In _Advances in Neural Information Processing Systems_ , 2015.
  * D’Amour et al. [2020] A. D’Amour, K. Heller, D. Moldovan, B. Adlam, B. Alipanahi, A. Beutel, C. Chen, J. Deaton, J. Eisenstein, M. D. Hoffman, F. Hormozdiari, N. Houlsby, S. Hou, G. Jerfel, A. Karthikesalingam, M. Lucic, Y. Ma, C. McLean, D. Mincu, A. Mitani, A. Montanari, Z. Nado, V. Natarajan, C. Nielson, T. F. Osborne, R. Raman, K. Ramasamy, R. Sayres, J. Schrouff, M. Seneviratne, S. Sequeira, H. Suresh, V. Veitch, M. Vladymyrov, X. Wang, K. Webster, S. Yadlowsky, T. Yun, X. Zhai, and D. Sculley.  Underspecification presents challenges for credibility in modern machine learning.  In _arXiv:2011.03395_ , 2020.
  * Finn et al. [2017] C. Finn, P. Abbeel, and S. Levine.  Model-agnostic meta-learning for fast adaptation of deep networks.  In _Proceedings of the International Conference on Machine Learning_ , pages 1126–1135, 2017.
  * Geiger et al. [2012] A. Geiger, P. Lenz, and R. Urtasun.  Are we ready for autonomous driving? The KITTI vision benchmark suite.  In _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition_ , 2012.
  * Guo et al. [2017] C. Guo, G. Pleiss, Y. Sun, and K. Q. Weinberger.  On calibration of modern neural networks.  In _Proceedings of the International Conference on Machine Learning_ , pages 1321–1330, 2017.
  * Innes et al. [2019] M. Innes, A. Edelman, K. Fischer, C. Rackauckas, E. Saba, V. Shah, and W. Tebbutt.  A differentiable programming system to bridge machine learning and scientific computing.  In _arXiv:1907.07587_ , 2019.
  * Kornblith et al. [2019] S. Kornblith, J. Shlens, and Q. V. Le.  Do better ImageNet models transfer better?  In _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition_ , pages 2661–2671, 2019.
  * Krizhevsky and Hinton [2009] A. Krizhevsky and G. Hinton.  Learning multiple layers of features from tiny images.  Technical report, University of Toronto, 2009.
  * Lwakatarea et al. [2020] L. Lwakatarea, A. Raj, I. Crnkovica, J. Bosch, and H. H. Olsson.  Large-scale machine learning systems in real-world industrial settings: A review of challenges and solutions.  _Information and Software Technology_ , 127, 2020.
  * Nushi et al. [2017] B. Nushi, E. Kamar, E. Horvitz, and D. Kossmann.  On human intellect and machine failures: Troubleshooting integrative machine learning systems.  In _Proceedings of AAAI Conference on Artificial Intelligence_ , 2017\.
  * Ozkaya [2020] I. Ozkaya.  What is really different in engineering AI-enabled systems?  _IEEE Software_ , 37(4), 2020.
  * Pan and Yang [2010] S. J. Pan and Q. Yang.  A survey on transfer learning.  _IEEE Transactions on Knowledge and Data Engineering_ , 22:1345–1359, 2010.
  * Patel et al. [2015] V. M. Patel, R. Gopalan, R. Li, and R. Chellappa.  Visual domain adaptation: A survey of recent advances.  _IEEE Signal Processing Magazine_ , 32:53–69, 2015.
  * Qi et al. [2018] C. R. Qi, W. Liu, C. Wu, H. Su, and L. J. Guibas.  Frustum pointnets for 3D object detection from RGB-D data.  In _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition_ , pages 918–927, 2018.
  * Scovel and Steinwart [2005] C. Scovel and I. Steinwart.  Fast rates for support vector machines.  In _Proceedings of the Conference on Learning Theory (COLT 2005)_ , pages 279–294, 2005.
  * Sculley et al. [2015] D. Sculley, G. Holt, D. Golovin, E. Davydov, T. Phillips, D. Ebner, V. Chaudhary, M. Young, J.-F. Crespo, and D. Dennison.  Hidden technical debt in machine learning systems.  In _Advances in Neural Information Processing Systems_ , 2015.
  * Shen et al. [2020] Y. Shen, Y. Xiong, W. Xia, and S. Soatto.  Towards backward-compatible representation learning.  In _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition_ , pages 6368–6377, 2020.
  * Shi et al. [2019] S. Shi, X. Wang, and H. Li.  PointRCNN: 3D object proposal generation and detection from point cloud.  In _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition_ , pages 770–779, 2019.
  * Shimodaira [2000] H. Shimodaira.  Improving predictive inference under covariate shift by weighting the log-likelihood function.  _Journal of Statistical Planning and Inference_ , 90:227–244, 2000.
  * Srivastava et al. [2020] M. Srivastava, B. Nushi, E. Kamar, S. Shah, and E. Horvitz.  An empirical analysis of backward compatibility in machine learning systems.  In _Proceedings of the ACM SIGKDD Conference on Knowledge Discovery and Data Mining_ , 2020.
  * Vapnik [1982] V. N. Vapnik.  _Estimation of Dependences Based on Empirical Data_.  Springer-Verlag, Berlin, 1982.
  * Wang et al. [2020] C.-Y. Wang, Y.-L. Chang, S.-T. Yang, D. Chen, and S.-H. Lai.  Unified representation learning for cross model compatibility.  In _Proceedings of the British Machine Vision Conference_ , 2020.
  * Wang et al. [2018] F. Wang, J. Decker, X. Wu, G. Essertel, and T. Rompf.  Backpropagation with callbacks: Foundations for efficient and expressive differentiable programming.  In _Advances in Neural Information Processing Systems_ , volume 31, 2018.
  * Wang et al. [2019] Y. Wang, W.-L. Chao, D. Garg, B. Hariharan, M. Campbell, and K. Q. Weinberger.  Pseudo-LiDAR from visual depth estimation: Bridging the gap in 3D object detection for autonomous driving.  In _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition_ , pages 8445–8453, 2019.
  * Yan et al. [2021] S. Yan, Y. Xiong, K. Kundu, S. Yang, S. Deng, M. Wang, W. Xia, and S. Soatto.  Positive-congruent training: Towards regression-free model updates.  In _arXiv:2011.09161_ , 2021.
  * Zhang [2004] T. Zhang.  Statistical behavior and consistency of classification methods based on convex risk minimization.  _The Annals of Statistics_ , 32:56–85, 2004.