# Modeling and frequency domain analysis of nonlinear compliant joints for a
passive dynamic swimmer

Juan Pablo Carbajal, Rafael Bayma, Marc Ziegler and Zi-Qiang Lang

(August 0 d , 2024)

###### Abstract

In this paper we present the study of the mathematical model of a real life
joint used in an underwater robotic fish. Fluid-structure interaction is
utterly simplified and the motion of the joint is approximated by Düffing’s
equation. We compare the quality of analytical harmonic solutions previously
reported, with the input-output relation obtained via truncated Volterra
series expansion. Comparisons show a trade-off between accuracy and
flexibility of the methods. The methods are discussed in detail in order to
facilitate reproduction of our results. The approach presented herein can be
used to verify results in nonlinear resonance applications and in the design
of bio-inspired compliant robots that exploit passive properties of their
dynamics. We focus on the potential use of this type of joint for energy
extraction from environmental sources, in this case a Kármán vortex street
shed by an obstacle in a flow. Open challenges and questions are mentioned
throughout the document.

##  1 Introduction

How much of the diverse behavior we observe in animals is a direct expression
of the dynamics of the individual’s body? In the last two decades, many
characteristics of animal locomotion on land were successfully linked to the
mechanical properties of the legs. The springy behavior observed in the
trajectories of the center of mass during running, walking and jumping can be
explained by the mechanical properties of the limbs and its tunning (Farley et
al., (1993); Farley, (1996); Ferris et al., (1998); Dickinson et al., (2000);
Kerdok et al., (2002); Moritz and Farley, (2003); Roberts and Azizi, (2011)).
The hypothesis that running, walking and jumping is tuned to the resonance of
the underlying mechanical system is strongly supported by experimental
evidence and by machines constructed based on this idea, the passive dynamic
walkers (Alexander, (1990); Ahlborn and Blake, (2002); McGeer, (1990);
Thompson and Raiber, (1989); Collins et al., (2005)). The consequence of such
setting is energy efficient performance and alleviation of the controller.

Due to similarities between running and swimming (Bejan and Marden, (2006);
Kokshenev, (2010)), it is not surprising (but not less exiting) that efficient
locomotion in fluids was reported to relay on the dynamics of the body of the
animal. Living trouts have been observed to exploit the energy in the flow
they inhabit to reduce their swimming efforts (Liao et al., (2003)). Later,
euthanized trouts performed passive self-propulsion when placed in the von
Kármán vortex street shed by an obstacle in a flow (Liao, (2007)). The
experimental results are supported by mathematical models, analytical and
numerical (Eldredge and Pisani, (2008); Kanso and Newton, (2009); Alben,
(2009)). These models do not fully agree with each other, however this is
expected due to the mathematical complexity of the interaction between
structures and fluids, of which the passive case is the worst scenario:
unprescribed motion of the interface boundary. This alone represents an open
challenge for the mathematical modeling community. Lest the challenge remains
unsolved for too long, robotic researchers design their machines with less
fluid-dynamic rigor, pursuing the first passive dynamic swimmer. The final
objective is to build a swimming machine that can perform at least as well as
fish, and at comparable power ratings (Harper et al., (1997); Lauder et al.,
(2007)).

For a robot to extract the energy in the surroundings, its mechanical
properties have to be tuned to the environmental energy storage. Stated this
way, the problem is one of energy harvesting, were nonlinear properties are
believed to be beneficial (Cottone et al., (2009)). Therefore, we need to pin
down the resonance characteristic of the actuators and joints to be used in
the robot that, as their biological counterparts, are generally nonlinear. As
if difficulty was lacking, the study of resonance of nonlinear systems has
suffer a very slow developmental process that started early in the 1960’s and
has yet not overcome its infancy. In a technical report from 1958 by
Brilliant, (1958) we read:

> Sometimes nonlinearity is avoided, not because it would have an undesired
> effect in practice, but simply because its effect cannot be computed.

Nowadays the situation is not completely different. However, we have more
powerful computers, new simulation methods and some novel uses of classical
tools promise to open the path ahead Peng et al., (2007); Vakakis et al.,
(2009).

In this paper we present the study of the joint of a robot fish from a
mathematical point of view. The amplitude of the oscillations of the joint in
response to periodic forcing is studied. In section 2 we briefly introduce the
real mechanical device and we move to its mathematical model in section 3. In
the same section the model for fluid-structure interaction is briefly
described. Approximated solutions methods are introduced in section 4. Results
are presented in section 5. Finally, we close the paper in section 6 with a
discussion on the implications of the results and the relevance of the
approach to the design of robots and the test of controllers based on
resonances.

##  2  A simple compliant joint

To extract energy from the environment, the robotic platform has to be
optimally driven by interaction forces. In this case, by the interaction
between the rigid body of the robot and the surrounding turbulent flow. The
feat can not be performed if the angle trajectories of the joints are fully
prescribed by the controller. Therefore, the joints of our robot are
compliant, i.e. the motion of the joint is not only defined by the controller,
but also by external actions. In Figure 1a we see the details of the joints of
the robot fish used in Ziegler et al., (2011). Each joint behaves as a
rotational spring. The restoring torque is generated when the relative angle
between the two connected bodies is not zero. The force producing the torque
is given by the extension of a linear spring fixed to the first body. The
spring is connected via an inelastic thread to an appendage of the second
body. When the deflection angle is zero, the extension of the spring is
minimum as well as the force it exerts. We call tension to this minimum force
value and it is referred with the letter $F$. Measurements of the torque for
$F=$0.73\text{\,}\mathrm{N}$$, are given in Figure 1b. The values of the
parameters used throughout this paper are given in Table 1. The two parameters
$r,d$ are distances that can be seen in the figures. The elastic constant of
the linear spring is $K$ and $I$ denotes the moment of inertia of the joint
around the axis of rotation. The linear specific damping coefficient of the
joint is denoted with $\zeta$. The parameters $\mathcal{Q}_{0}$ and $\Omega$
correspond to the amplitude and frequency of the external forcing,
respectively.

Name | Value
---|---
$r$ |  $20.24\pm 0.02\text{\,}\mathrm{m}$ m
$d$ |  $27.68\pm 0.02\text{\,}\mathrm{m}$ m
$K$ |  $81\pm 1\text{\,}\mathrm{N}$/ m
$I$ |  $(3.1\pm 0.1)\times 10^{-5}\text{\,}\mathrm{kg}$ m ^ 2
$\zeta\cdot I$ |  $(2.2\pm 0.1)\times 10^{-4}\text{\,}\mathrm{N}$ m s
$\mathcal{Q}_{0}\cdot I$ |  $1\times 10^{-4}\text{\,}\mathrm{N}$ m
$\nicefrac{{\Omega}}{{2\pi}}$ | $\left(0,3\right]\text{\,}\mathrm{Hz}$
Table 1: Value of the parameters used here and in the model studied in Ziegler
et al., (2011) . Figure 1: Description of the joints in WandaX. a) Details of
the joints used in Ziegler et al., (2011), note the axis of rotation and the
appendage. b) Measured torques applied to the joint under controlled
deflection. Replacing the parameters given in Table 1 into Eq. (1) the fit
provides $F=$0.73\pm 0.05\text{\,}\mathrm{N}$$.

##  3 Mathematical model

A geometrical representation of the joint is given in Figure 2. The parameters
$r$,$d$ and $K$ are fixed at construction time and always $d>r$ (Table 1). The
tension in the spring at its shortest length ($F$), is the controlled
parameter and a servomotor can change it dynamically.

Figure 2: Schematic of the rotational spring used to derive Eq. (1).

The torque applied to the bodies connected to the joint is thus,

$\begin{split}\tau=&\\\
&\frac{K\left(\sqrt{\epsilon^{2}+4S\sin^{2}\frac{\theta}{2}}-\epsilon\right)+F}{\sqrt{\epsilon^{2}+4S\sin^{2}\frac{\theta}{2}}}S\sin\theta.\end{split}$
(1)

Where we have defined the parameters $\epsilon=d-r$ and $S=rd$. The formula is
obtained by calculating the deformation of the spring as a function of the
deflection angle of the appendage ($\theta$). The deformation is inside
parenthesis in the numerator of equation (1) and when multiplied by the
stiffness $K$, it gives the force due to the deformation of the linear spring.
The outer factors come from the product of the force and the moment arm. Note
that the the torque $\tau$ is linear in the controlled input $F$.

For $\theta\ll 1$ the third order Taylor expansion gives

$\tau(\theta,F)=\kappa(F)\theta+\alpha(F)\theta^{3}+\mathcal{O}(\theta^{5}).$
(2)

Where

$\displaystyle\kappa(F)$ $\displaystyle=$ $\displaystyle\frac{SF}{\epsilon},$
(3) $\displaystyle\alpha(F)$ $\displaystyle=$
$\displaystyle\frac{SF}{\epsilon}\left[\frac{S}{2\epsilon^{2}}\left(\frac{\epsilon
K}{F}-1\right)-\frac{1}{6}\right].$ (4)

And the equation of motion of the deflection angle is

$\begin{split}\ddot{\theta}+\zeta\dot{\theta}+\frac{\tau(\theta,F)}{I}&\approx\\\
\ddot{\theta}+\zeta\dot{\theta}+k(F)\theta&+a(F)\theta^{3}\;=\;\Gamma.\end{split}$
(5)

Where $I$ denotes the moment of inertia around the axis of rotation. The
specific damping is given by $\zeta$, and we have defined
$k=\nicefrac{{\kappa}}{{I}}$, $a=\nicefrac{{\alpha}}{{I}}$. $\Gamma$ is the
specific net effect of all other external torques acting on the joint. The
approximating equation of motion is the well studied Düffing’s equation.

To quantify the error introduced by the approximation, we calculated the angle
at which the difference between the torque produced by equation (1) and
equation (2) is equal to a reference error given by $\Delta\tau=r\Delta F$,
where $\Delta F=$0.05\text{\,}\mathrm{N}$$ is a reasonable resolution for a
force sensor working in a $10\text{\,}\mathrm{N}$ range. These angles are
plotted in Figure 3 for different values of the tension.

Figure 3: Approximation error. Values of the deflection of the joint for which
the error reaches the reference value given in the text. The maximum
deflection has a peak close to $F_{3}=$(5.90\pm 0.01)\times
10^{-1}\text{\,}\mathrm{N}$$ meaning that terms of degree $>3$ almost cancel
each other.

###  3.1 Hardening, linear and softening spring

As can be seen from equation (2), $\alpha(F)$ modulates the intensity of the
cubic non-linear term. This term vanishes when

$F=F_{*}=\frac{\epsilon K}{\cfrac{\epsilon^{2}}{3S}+1}.$ (6)

rendering a linear spring for the angles where the third order approximation
is valid. For bigger values of $F$ the spring will be softening ($\alpha<0$)
and for smaller values it will be a hardening spring ($\alpha>0$) (Fig. 4).
However, the full expression of the torque (Eq. (1)) contains higher order
terms. Hence, the linear behavior will be even more evident if the higher
order terms cancel each other. In Figure 4 curves of torque versus angle for
several values of $F$ are shown. In particular we show the curve for which up
to seventh order nonlinearities give a minimum contribution (found via
optimization $F_{0}=0.9F_{*}$) together with the curve at $F=F_{*}$. These
illustrates the power of the actuation chose, since we can control the
dynamical properties of a virtual rotational spring at the joint (more details
in Ziegler et al., (2011)).

Figure 4: Plots of the torque function for values of tension
$F=[0.5,\,1,\,0.9,\,1.5]\cdot F_{*}$. Note the hardening behavior for small
tension and softening behavior for higher tension. Almost linear curves are
found for $F=F_{*}$ and $F=F_{0}$, as a reference the linear curve is shown.

###  3.2 Forcing model

As explained before, the external torques acting on the joint are due to
fluid-structure interaction. This kind of interaction still poses a great
challenge for the mathematical modeling community. In the search for
simplified models of the forces generated in the interaction between flexible
bodies and turbulent flows, we found the work of Kanso and Newton, (2009) and
Alben, (2009) instructive. Using equation (3.8) given by Alben, (2009), we can
calculate the pressure difference on the boundary of a slender body in a
vortex street. If the width of the body is bigger than the separation among
vortices and it is placed in the middle of the vortex street, the forces on
the body can be approximated by a sine function $f(t)=f_{0}\sin(\Omega
t+\phi)$, where $t$ is time and $\phi$ is an initial offset that sets the
balance between the sine and cosine behavior of the forcing. The frequency
$\Omega$ is proportional to the speed of the flow plus the vorticity of the
vortices (assumed to be equal for each vortex) scaled by a factor that depends
on the geometrical properties of the wake. The amplitude $f_{0}$ is
proportional to the density of the fluid and the square of the vorticity of
each vortex. With a few additional assumptions, the torque acting on the joint
can be made proportional to this force. Here we adopt this over-simplified
forcing model to avoid diverting the attention of the reader from the core
ideas of our work. In this manner we postpone a detailed study with a more
elaborated forcing model.

##  4  Solution methods

In this section we present two independent methods to estimate the amplitude
of oscillations of the joint under periodic forcing. The first method turns
out to be excellent for estimating the amplitude of the first harmonic. The
second method is based on a Volterra series expansion of the Duffing equation.
It allows to study the response of the joint to more general forcing
conditions, but the tension range for which it is valid is smaller.

###  4.1 Harmonic solutions of Düffing’s equation

Under periodic forcing $\Gamma=\mathcal{Q}_{0}\sin(\Omega t)$, equation (5)
has been extensively studied (see Holmes and Rand, (1976); Luo and Han, (1997)
and references therein). Following these analyses, we show here how to
maximize the amplitude of the periodic response of the joint by tunning the
tension parameter.

The key point of the analysis in Luo and Han, (1997), is that we search for
the amplitude of solutions of (5) that are periodic (this rules out sub-
harmonics, supra-harmonics and chaotic motion), $\theta(t)=A\sin(\Omega
t+\psi)$. Under this assumption it can be shown that the amplitude $A$ of such
solutions is given by the roots of the polynomial,

$\begin{split}A^{3}\frac{9a^{2}}{16}+A^{2}\frac{3\left(k-\Omega^{2}\right)a}{2}+&\\\
A\left[\left(\Omega^{2}+\zeta^{2}-2k\right)\Omega^{2}+k^{2}\right]&-\mathcal{Q}_{0}^{2}=0.\end{split}$
(7)

The roots of this polynomial can be obtained analytically using a computer
algebra system as Maxima, (2009) and they establish the relation between the
amplitude of the oscillations and the parameters of the equation,
$A\left(k,a,\Omega,\mathcal{Q}_{0}\right)$. In the case at hand, we have
$k(F)$ and $a(F)$, therefore $A\left(F,\Omega,\mathcal{Q}_{0}\right)$. For a
given value of the parameters, only one root corresponds to the observed
amplitude (there are unstable amplitudes). This implies, that is not enough to
look at the roots, but we must also check their stability. This adds some
complexity to the evaluation of the results that adds to the limitation of
pure harmonic inputs.

###  4.2  Volterra series expansion

When the amplitude of the forcing is sufficiently small, the behavior of
Duffing’s oscillator in the neighborhood of the origin can be described by
polynomial Volterra functionals(Theorem 3.1 of Rugh, (1981)). This means that
the angle of the oscillations can be described by an expansion of the form

$\theta(t)\approx\sum_{i=1}^{n}y_{i}(t),$ (8)

where $n$ is the order of the expansion and each term is the multi-dimensional
integral

$y_{i}(t)=\int_{-\infty}^{+\infty}h_{i}(t-\sigma_{1},\ldots,t-\sigma_{i})\\\
\Gamma(\sigma_{1})\ldots\Gamma(\sigma_{i})\mathrm{d}\sigma_{1}\ldots\mathrm{d}\sigma_{i}.$
(9)

where $\Gamma(t)$ is the forcing signal and we wish to determine the kernel
functions $\left\\{h_{i}\right\\}_{i=1}^{n}$. However, for the objective at
hand, it is more useful to calculate the Fourier transform $Y(\omega)$ of
these functionals. We have not found this process described in the literature,
hence we describe it succinctly. We used the inverse Fourier transform
definition on the input signal and substitute it into the functional
definition Eq. (9). By inverting the order of integration and splitting
variables, we make the kernel transform $H(\omega)$ appear. Then, we compare
this expression to the formal definition of the inverse Fourier transform of
$Y(\omega)$ and isolate the desired result. This process can be applied for
the first three orders, but it gets cumbersome for higher ones. Therefore,
before showing the results, we will briefly introduce some notation used here
to simplify the presentation.

An ordered set of arguments $(s_{1},s_{2},\ldots,s_{n})$ will be denoted
$s_{1:n}$. In general we have,

$(s_{k},s_{k+1},\ldots,s_{n})=s_{k:n}\quad k\leq n.$

For example, the fifth order kernel $H_{5}(s_{1},s_{2},s_{3},s_{4},s_{5})$
will be written

$\displaystyle H_{5}(s_{1:5})=$ $\displaystyle-3a\,H_{1}\left(\Sigma
s_{1:5}\right)H_{1}(s_{1})H_{1}(s_{2})H_{3}(s_{3:5}),$

where $\Sigma s_{1:5}=\sum_{i=1}^{5}s_{i}$. For the integration variable in
multiple integrals we will write
$\mathrm{d}s_{1:n}=\mathrm{d}s_{1}\mathrm{d}s_{2}\cdots\mathrm{d}s_{n}$. For
multiple summations over $n$ indexes we will write
$\sum_{k_{1:n}}=\sum_{k_{1}}\ldots\sum_{k_{n}}$. Having defined the notation,
we continue with our presentation.

For zero initial conditions, $(\theta,\dot{\theta})=0$, and working in the
frequency domain, each kernel can be obtained recursively based on lower order
kernels (details of the calculation are given in Peyton-Jones and Billings,
(1989)). The first order kernel derived from the model in Eq. (5), is
equivalent to the frequency response function (also known as transfer
function) of a linear second order system,

$H_{1}(s)=\frac{1}{s^{2}+\zeta s+k}.$ (10)

The following equations present the Volterra kernels up to seventh order using
recurrent relations,

$\displaystyle\begin{gathered}\begin{split}H_{3}(s_{1:3})=-a&H_{1}\left(\Sigma
s_{1:3}\right)\cdot\\\ &H_{1}(s_{1})H_{1}(s_{2})H_{1}(s_{3}),\end{split}\\\
\begin{split}H_{5}(s_{1:5})=-3a&H_{1}\left(\Sigma s_{1:5}\right)\cdot\\\
&H_{1}(s_{1})H_{1}(s_{2})H_{3}(s_{3:5}),\end{split}\\\
\begin{split}H_{7}(s_{1:7})=&-3aH_{1}\left(\Sigma s_{1:7}\right)\cdot\\\
&\Big{[}\;H_{1}(s_{1})H_{1}(s_{2})H_{5}(s_{3:7})+\\\
&H_{1}(s_{1})H_{3}(s_{2:4})H_{3}(s_{5:7})\;\Big{]}.\end{split}\end{gathered}$
(14)

Due to the symmetry of the eq. (5), all even order kernels are zero. Note that
these kernels are valid for any second order cubic oscillator. The recursive
formulas were obtained using the probing method of Peyton-Jones and Billings,
(1989), which proceeds as follows. We start from a polynomial nonlinear
ordinary differential equation (as Eq. (5)) and assume the response can be
represented by Volterra series, i.e. Eq. (8). We substitute this into the
system’s differential equation, and equate similar terms. This last step
produces equations relating different order kernels and the expressions are
always recursive (we obtain Eqs. (14) in our case). The factorial in the
recursive relations appear from the multinomial expansion of nonlinear
polynomial terms.

We proceed to calculate the Fourier transform of the output given by Eq. (9)
subject to harmonic forcing. We will use Eqs. (14) in order to obtain an
explicit relation between the parameters of our model
$k,a,\zeta,\mathcal{Q}_{0},\Omega,F$ and the amplitude of the first harmonic
of the output.

(a)

(b)

Figure 5: Amplitude of periodic response in the plane $(F,\Omega)$.
LABEL:sub@fig:JointORFLuo Amplitude of the oscillations according to Eq. (7).
The line of maxima is marked with dashes. The vertical dotted lines indicate
selected values of $\Omega$. LABEL:sub@fig:ORFError Relative differences
between the predicted amplitudes by the method described in section 4.1 and
the Volterra expansion described in section 4.2.

We recall that the Fourier transform of a general bandwith limited periodic
function is a scaled Dirac Comb (Shah function, Impulse train, Dirac train,
etc.) (see Schetzen, (1980)),

$\displaystyle\gamma\left(\omega\right)=\sum_{l=-L}^{+L}X_{i}\delta\left(\omega-\Omega
l\right)$ (15)

where $L$ is a positive integer, $X_{i}$ is the complex amplitude of the
$i$-th harmonic and $\delta(\cdot)$ is the impulse function (Dirac delta
function). The Fourier transform of each of the output terms in (8) is

$\begin{gathered}Y_{i}(\omega)=\left(2\pi\right)^{1-i}\int_{-\infty}^{\infty}H_{i}\left(\Delta\omega_{i}\right)\cdot\\\
\gamma(w_{1})\cdots\gamma(w_{i-1})\gamma\left(\omega-\Sigma
w_{1:i-1}\right)\mathrm{d}w_{1:i-1}=\\\
(2\pi)^{1-i}\sum_{l_{1:i}=-L}^{+L}B_{l_{1}:l_{i}}\delta\left(\omega-\Omega\Sigma
l_{1:i}\right),\end{gathered}$ (16)

where we have used Eq. (15) and
$\Delta\omega_{i}=(w_{1},\ldots,w_{i-1},\omega-\Sigma w_{1:i-1})$ and,

$\displaystyle B_{l_{1}:l_{i}}=H_{i}\left(\Omega l_{1:i}\right)X_{l_{1}}\cdots
X_{l_{i}}.$ (17)

The integral in (16), has an interesting geometric interpretation in terms of
convolutions in hyperplanes, we refer the interested reader to Lang and
Billings, (1996). Therein, the output frequency range of nonlinear systems
that are representable by Volterra series is analytically calculated.
Additionally, in Eq.(16) of that paper the frequency spectrum of the output
signal is represented as the superposition of contributions from the
nonlinearities. In the case studied herein, the input has only one frequency,
therefore $L=1$, $X_{0}=0$, $X_{\pm 1}=\pm j\mathcal{Q}_{0}\pi$ with $j$ the
imaginary unit. Replacing these values in all the equations and using the
relations in Eq. (14) we obtain the desired result,

$Y(\Omega)=\frac{-j\pi}{64}\Big{\\{}3a^{3}\mathcal{Q}_{0}^{7}\big{[}\phantom{aaaaaaaaaaaaaaaaaaaaaaaaa}\\\
2H_{1}(3\Omega)H_{1}(-3\Omega)H_{1}^{3}(-\Omega)H_{1}^{5}(\Omega)\\\
+6H_{1}(3\Omega)H_{1}^{4}(-\Omega)H_{1}^{5}(\Omega)\phantom{aaaaaaa}\\\
+6H_{1}^{2}(3\Omega)H_{1}^{3}(-\Omega)H_{1}^{5}(\Omega)\phantom{aaaaaaa}\\\
+3H_{1}(-3\Omega)H_{1}^{4}(-\Omega)H_{1}^{5}(\Omega)\phantom{aaaaaa}\\\
-15H_{1}(3\Omega)H_{1}^{3}(-\Omega)H_{1}^{6}(\Omega)\phantom{aaaaaaa}\\\
+45H_{1}^{3}(-\Omega)H_{1}^{7}(\Omega)\phantom{aaaaaaaaaaaa}\\\
+45H_{1}^{4}(-\Omega)H_{1}^{6}(\Omega)+18H_{1}^{5}(-\Omega)H_{1}^{5}(\Omega)\;\big{]}\phantom{aaaaaaaaaaa}\\\
-12a^{2}\mathcal{Q}_{0}^{5}\big{[}H_{1}(3\Omega)H_{1}^{2}(-\Omega)H_{1}^{4}(\Omega)\phantom{aa}\\\
+6H_{1}^{2}(-\Omega)H_{1}^{5}(\Omega)\phantom{aaaaaa}\\\
+3H_{1}^{3}(-\Omega)H_{1}^{4}(\Omega)\;\big{]}\phantom{aaaaa}\\\
+48a\mathcal{Q}_{0}^{3}H_{1}(-\Omega)H_{1}^{3}(\Omega)-64\mathcal{Q}_{0}H_{1}(\Omega)\Big{\\}}.$
(18)

Where $H_{1}(x)$ is given by Eq. (10) with $s=jx$. The amplitude of the
oscillation is obtanied by taking the double of the modulus of the complex
number $Y(\Omega)/\pi$.

##  5  Results

Figure 5a shows the amplitude of periodic oscillations in the $(F,\Omega)$
plane, according to Eq. (7). The line of maxima is shown with dashes. This
line describes the value of the tension that produces maximum amplitude for a
given forcing frequency. In Fig. 5b we plot the logarithm of the relative
difference between the amplitudes given by (7) and (18). The two approximation
differ for regions of low frequency and low tension where the system is most
nonlinear (Luo and Han, (1997)). To compare with the observed amplitude
(obtained by simulating (5) without any approximation), we extracted the
curves of amplitude against tension for
$\nicefrac{{\Omega}}{{2\pi}}=$0.5,1.5\text{\,}\mathrm{Hz}$$ (vertical dotted
lines in Fig. 5). These curves are shown in Figure 6 together with the
amplitude calculated using Eq. (18) corresponding to the Volterra kernel
expansion. For the $1.5\text{\,}\mathrm{Hz}$ frequency, both models predict
the simulated amplitude accurately. For the $0.5\text{\,}\mathrm{Hz}$
frequency, the amplitude predicted by Eq. (7) drifts away from the simulated
value for lower tensions. The amplitude calculated using the Volterra
expansion diverges for low tensions at this frequency.

Figure 6:  Amplitude of periodic response of the joint for forcing with
frequencies $\Omega=[$0.5,1.5\text{\,}\mathrm{]}$\mathrm{Hz}$. Amplitudes
according to Eq. (7) in solid line, amplitudes obtained from the Volterra
expansion Eq. (18) in dashes and amplitudes from simulations without
approximation in circles.

In the same figure we show the amplitudes measured in simulations using the
exact expression for the torque (1) (Octave (Eaton, (2002)) function ode45,
absolute and relative tolerances, $10^{-6}$). The amplitude of the
oscillations is obtained using the absolute value of the analytic signal of
the angle (Octave function hilbert). The agreement between the approximated
results and the simulated ones is noteworthy.

These results show that the approximation from Luo and Han, (1997) can be used
for a tension controller designed for a joint of this kind, to obtain periodic
responses with the same frequency as the forcing. The Volterra approximation
fails when the system becomes strongly nonlinear, however these expansion can
be used for more general responses or when the inputs have a more complicated
frequency spectrum (as in a realistic scenario).

##  6 Discussions and conclusion

Though the quote from Brilliant, (1958) remains valid, we used the knowledge
about Düffing’s equation to understand analytically a compliant joint. This
gives a corner stone to study more complicated setups (e.g. chain of joints,
as in the robotic fish) and define the robust engineering design of robots
with the desired resonance properties. Additionally, the results showed here
provide a reference solution to the problem of finding the right pretension
for a given external forcing, related to the idea of adaptive controllers. For
a forcing with a sufficiently slow varying frequency, the tension could be
adjusted to maintain the amount of energy transferred into the system to its
maximum possible value, i.e. keep the system close to the line of maxima in
Fig. 5a. With the results herein the adjustment could be done by direct
calculations. However, other methods like the frequency oscillators(AFO)
(Buchli et al., (2006)) were used to achieve a similar objective. Applications
of the latter has been reported in Buchli and Ijspeert, (2008), however, due
to lack of theoretical results on the resonance frequency of the nonlinear
platform studied therein, the theoretical solution to the problem was unknown
and a grounded evaluation of the performance was not possible. Approximations
as the one presented here could be used to verify the results generated by the
adaptive oscillators.

In our results we exclude any information about oscillations with frequencies
different from that of the forcing since that is the case where Eq. (7) is
valid. The higher order harmonics in the system’s output could be determined
and compared with the simulated results. Traditional harmonic balance method
can not be used to find these results. Additionally, a more complete panorama
of the frequency response of the system could be obtained with a nonlinear
analysis method as the NOFRF (proposed in Lang and Billings, (2005) and
showcased in Peng et al., (2007)), which is based on Volterra series
expansions as the one presented in this work. NOFRF give the response of the
system to a given set of inputs, therefore the natural ensemble of forcing
signals must be known to acquire useful information. This ensemble has to be
compiled from (or modeled based on) fluid-structure interactions data.
Comparing the present results with the information provided by NOFRF will be
the objective of our next work. Proved that such analysis is helpful for the
design of compliant robots, we will extend the analysis to biped and quadruped
robots. In that situation the forcing comes from the ground reaction forces,
that are generally not smooth (impacts) and contain a rich frequency spectrum.
We consider those scenarios as highly complex problems to be addressed with a
mature toolbox of methods.

Concluding, we have provided a method that allows the construction of a
controller that would maximize the harvesting of environmental energy, under
the circumstances defined herein. The method will be extended to cover more
realistic situations than the simplified periodic forcing, including forcing
with wider frequency spectrum. Additionally, the same methodology presented
here can be used to verify results obtained using heuristic methods as AFO or
neural networks. Our results can be easily validated and we invite our
colleagues to do it.

## Acknowledgements

The authors would like to thank Cecilia Tapia Siles for sharing insights and
information about passive swimmers.

#### Funding.

The research leading to these results has received funding from the European
Community’s Seventh Framework Programme FP7/2007-2013-Challenge 2 –Cognitive
Systems, Interaction, Robotics– under [grant agreement No 248311-AMARSi]. RB
is founded by the CAPES Foundation, Ministry of Education of Brazil. MZ
received funding for this work from the SNSF [project no. 122279] (From
locomotion to cognition).

#### Author contributions.

JPC, RB and ZQL worked on the mathematics of the models and in the
simulations. MZ designed and built the joint studied in this work. All authors
collaborated for the production of the manuscript.

## References

  * Ahlborn and Blake, (2002) Ahlborn, B. K. and Blake, R. W. (2002).  Walking and running at resonance.  Zoology _105_ , 165–74.
  * Alben, (2009) Alben, S. (2009).  On the swimming of a flexible body in a vortex street.  J. Fluid Mech. _635_ , 27.
  * Alexander, (1990) Alexander, R. M. (1990).  Three Uses for Springs in Legged Locomotion.  Int. J. Robot. Res. _9_ , 53–61.
  * Bejan and Marden, (2006) Bejan, A. and Marden, J. H. (2006).  Unifying constructal theory for scale effects in running, swimming and flying.  J. Exp. Biol. _209_ , 238–48.
  * Brilliant, (1958) Brilliant, M. B. (1958).  Theory of the analysis of nonlinear systems.  Technical Report 345 MIT - Research laboratory of electronics.
  * Buchli and Ijspeert, (2008) Buchli, J. and Ijspeert, A. J. (2008).  Self-organized adaptive legged locomotion in a compliant quadruped robot.  Autonomous Robots _25_ , 331.
  * Buchli et al., (2006) Buchli, J., Righetti, L. and Ijspeert, A. J. (2006).  Engineering entrainment and adaptation in limit cycle systems : From biological inspiration to applications in robotics.  Biol. Cybern. _95_ , 645–64.
  * Collins et al., (2005) Collins, S., Ruina, A., Tedrake, R. and Wisse, M. (2005).  Efficient bipedal robots based on passive-dynamic walkers.  Science _307_ , 1082–5.
  * Cottone et al., (2009) Cottone, F., Vocca, H. and Gammaitoni, L. (2009).  Nonlinear Energy Harvesting.  Phys. Rev. Lett. _102_ , 80601.
  * Dickinson et al., (2000) Dickinson, M. H., Farley, C. T., Full, R. J., Koehl, M. A. R., Kram, R. and Lehman, S. (2000).  How Animals Move: An Integrative View.  Science _288_ , 100–106.
  * Eaton, (2002) Eaton, J. W. (2002).  GNU Octave Manual.  Network Theory Limited, http://www.octave.org.
  * Eldredge and Pisani, (2008) Eldredge, J. D. and Pisani, D. (2008).  Passive locomotion of a simple articulated fish-like system in the wake of an obstacle.  J. Fluid Mech. _607_.
  * Farley, (1996) Farley, C. T. (1996).  Leg stiffness and stride frequency in human running.  J. Biomech. _29_ , 181–186.
  * Farley et al., (1993) Farley, C. T., Glasheen, J. and McMahon, T. A. (1993).  Running springs: speed and animal size.  J. Exp. Biol. _185_ , 71–86.
  * Ferris et al., (1998) Ferris, D. P., Louie, M. and Farley, C. T. (1998).  Running in the real world: adjusting leg stiffness for different surfaces.  Proc. R. Soc. Lond. [Biol] _265_ , 989–94.
  * Harper et al., (1997) Harper, K., Berkemeier, M. and Grace, S. (1997).  Decreasing the energy costs of swimming robots through passive elastic elements.  In Proc. ICRA 1997 vol. 3, pp. 1839 –1844,.
  * Holmes and Rand, (1976) Holmes, P. and Rand, D. (1976).  The bifurcations of duffing’s equation: An application of catastrophe theory.  J. Sound Vib. _44_ , 237–253.
  * Kanso and Newton, (2009) Kanso, E. and Newton, P. K. (2009).  Passive locomotion via normal-mode coupling in a submerged spring-mass system.  J. Fluid Mech. _641_ , 205.
  * Kerdok et al., (2002) Kerdok, A. E., Biewener, A. A., McMahon, T. A., Weyand, P. G. and Herr, H. M. (2002).  Energetics and mechanics of human running on surfaces of different stiffnesses.  J. Appl. Physiol. _92_ , 469–78.
  * Kokshenev, (2010) Kokshenev, V. B. (2010).  Key principle of the efficient running, swimming, and flying.  EPL (Europhysics Letters) _90_ , 48005.
  * Lang and Billings, (1996) Lang, Z. Q. and Billings, S. A. (1996).  Output frequency characteristics of nonlinear systems.  International Journal of Control _64_ , 1049–1067.
  * Lang and Billings, (2005) Lang, Z.-Q. and Billings, S. a. (2005).  Energy transfer properties of non-linear systems in the frequency domain.  International Journal of Control _78_ , 345–362.
  * Lauder et al., (2007) Lauder, G. V., Anderson, E. J., Tangorra, J. and Madden, P. G. A. (2007).  Fish biorobotics: kinematics and hydrodynamics of self-propulsion.  J. Exp. Biol. _210_ , 2767–80.
  * Liao, (2007) Liao, J. C. (2007).  A review of fish swimming mechanics and behaviour in altered flows.  Philos. T. Roy. Soc. B _362_ , 1973–93.
  * Liao et al., (2003) Liao, J. C., Beal, D. N., Lauder, G. V. and Triantafyllou, M. S. (2003).  Fish exploiting vortices decrease muscle activity.  Science _302_ , 1566–9.
  * Luo and Han, (1997) Luo, A. C. J. and Han, R. P. S. (1997).  A quantitative stability and bifurcation analyses of the generalized duffing oscillator with strong nonlinearity.  J. Franklin I. _334_ , 447–459.
  * Maxima, (2009) Maxima (2009).  Computer Algebra System. Version 5.18.1.  http://maxima.sourceforge.net/.
  * McGeer, (1990) McGeer, T. (1990).  Passive Dynamic Walking.  Int. J. Robot. Res. _9_ , 62–82.
  * Moritz and Farley, (2003) Moritz, C. T. and Farley, C. T. (2003).  Human hopping on damped surfaces: strategies for adjusting leg mechanics.  Proc R Soc Lond [Biol] _270_ , 1741–6.
  * Peng et al., (2007) Peng, Z., Lang, Z.-Q. and Billings, S. (2007).  Resonances and resonant frequencies for a class of nonlinear systems.  J. Sound Vib. _300_ , 993–1014.
  * Peyton-Jones and Billings, (1989) Peyton-Jones, J. C. and Billings, S. A. (1989).  Recursive algorithm for computing the frequency response of a class of non-linear difference equation models.  International Journal of Control _50_ , 1925–1940.
  * Roberts and Azizi, (2011) Roberts, T. J. and Azizi, E. (2011).  Flexible mechanisms: the diverse roles of biological springs in vertebrate movement.  J. Exp. Biol. _214_ , 353–361.
  * Rugh, (1981) Rugh, W. J. (1981).  Nonlinear system theory : the Volterra/Wiener approach.  Baltimore, London, Johns Hopkins University Press.
  * Schetzen, (1980) Schetzen, M. (1980).  The Volterra and Wiener theories of nonlinear systems.  Wiley-Interscience.
  * Thompson and Raiber, (1989) Thompson, C. M. and Raiber, M. H. (1989).  Passive Dynamic Running pp. 74–83.  Berlin/Heidelberg: Springer-Verlag.
  * Vakakis et al., (2009) Vakakis, A. F., Gendelman, O. V., Bergman, L. A., McFarland, D. M., Kerschen, G. and Lee, Y. S. (2009).  Nonlinear Targeted Energy Transfer in Mechanical and Structural Systems, vol. 156, of Solid mechanics and its applications.  Springer Netherlands, Dordrecht.
  * Ziegler et al., (2011) Ziegler, M., Hoffmann, M., Carbajal, J. P. and Pfeifer, R. (2011).  Varying Body Stiffness for Aquatic Locomotion.  In Proc. ICRA 2011, Shanghai, China.


# Altering level shifts and spontaneous decay rates of distant atoms
using partially-transparent asymmetric mirror interfaces

Nicholas Furtak-Wells  The School of Physics and Astronomy, University of
Leeds, Leeds LS2 9JT, United Kingdom  Benjamin Dawson  The School of Physics
and Astronomy, University of Leeds, Leeds LS2 9JT, United Kingdom  School of
Chemical and Process Engineering, University of Leeds, Leeds LS2 9JT, United
Kingdom  Thomas Mann  School of Chemical and Process Engineering, University
of Leeds, Leeds LS2 9JT, United Kingdom  Gin Jose  School of Chemical and
Process Engineering, University of Leeds, Leeds LS2 9JT, United Kingdom  Almut
Beige  The School of Physics and Astronomy, University of Leeds, Leeds LS2
9JT, United Kingdom

###### Abstract

In three dimensions, dipole-dipole interactions which alter atomic level
shifts and spontaneous decay rates only persist over distances comparable to
the wavelength of the emitted light. To provide novel tools for quantum
technology applications, like quantum sensing, many attempts have been made to
extend the range of these interactions. In this paper we show that this can be
achieved with the help of partially-transparent asymmetric interfaces without
involving negative index metamaterials. Suppose two atoms are placed on
opposite sides of the interface, each at the position of the mirror image of
the other. In this case, their emitted light interferes exactly as it would
when the atoms are right next to each other. Hence their dipole-dipole
interaction assumes an additional maximum, even when the actual distance of
the atoms is several orders of magnitude larger than their transition
wavelength.

##  I Introduction

In 1982, Scully and Drühl [1] proposed a double-slit experiment in which the
slits are two two-level atoms. As illustrated in Fig. 1(a), the atoms are kept
at a constant distance, are continuously driven by laser light and emit
photons at a constant rate. When their distance is comparable to the
wavelength of the emitted light, an interference pattern forms on a far-away
screen. Averaged over many photons, this pattern very closely resembles the
interference pattern of classical double-slit experiments. It only disappears
when information about the origin of each photon becomes available [2]. As in
classical two-slit interference experiments, the distance between the
intensity minima and maxima depends on the distance between the atoms [3].

Figure 1: (a) Schematic view of the double-slit experiment of Eichmann et al.
[2] in which the slits are two laser-driven atoms. When the actual distance
$|x_{a}-x_{b}|$ of the atoms is of the same order of magnitude as the
wavelength $\lambda_{0}$ of the emitted light, the radiation from the atoms
interferes on a far-away screen, thereby causing a short-range dipole-dipole
interaction. (b) Schematic view of two atoms on opposite sides of a partially-
transparent mirror. In case of laser driving, the radiation from the atoms
interferes in an analogous way, although with reduced visibility, when the
distance of atom $a$ from the mirror image of atom $b$, i.e.
$|x_{a}|-|x_{b}|$, is comparable to $\lambda_{0}$. Now the interference is an
indication for the possible presence of an ultra-long range dipole-dipole
interaction. (c) Upon reflection, a partially-transparent interface adds phase
factors to the radiation coming from the atoms. As we shall see below, in case
of a symmetric mirror, these phase factors unfortunately cancel the predicted
atomic interactions. These require the presence of a partially-transparent
asymmetric interface which can be realised for example with a help of a thin
metallic layer with a different surface roughness on both sides.

When this two-atom double-slit experiment was first performed by Eichmann et
al. in 1993 [2], it raised many questions, like, how can spontaneously emitted
photons interfere [4]. A closer look at the experiment shows that it is best
not to think of the atoms as sources of spontaneously emitted photons but as
continuously radiating dipole antennae [5]. Both atoms constantly transfer
energy into the surrounding free radiation field which only manifests itself
as “individual photons” upon detection [6, 7]. When an individual photon is
registered on a photographic plate, it contains in general energy from both
atoms. Depending on its direction of emission, radiation either interferes
constructively or destructively, thereby resulting in a spatial dependence of
the intensity of the emitted light. Moreover, interference effects result in a
spatial dependence of first and second order photon correlations [8, 9, 10].
By now, the interference of light from distant atoms is relatively well
understood and has already found applications in distributed quantum computing
[12, 11, 13, 14, 15].

Different from the classical case, interference in the two-atom double-slit
experiment depends in addition on the internal state of the slits. Different
entangled atomic states radiate light in different preferred directions [16].
Suppose two atoms are right next to each other and share a single energy
quantum. If the atoms in Fig. 1(a) are in their symmetric state, all of their
emitted light interferes constructively, their coupling to the free radiation
field is collectively enhanced and a photon is emitted at twice the usual
rate. However, if the atoms are in their anti-symmetric state, their efforts
to transfer their energy into the free radiation field cancel each other out
and the spontaneous decay rate tends to zero. At larger atomic distances, the
alterations of the spontaneous decay rates of certain atomic states are
synonymous with dipole-dipole interactions [17, 18, 19, 20]. The key to atomic
ultralong-range interactions is therefore the recreation of the interference
effects of the original two-atom double-slit experiment, even when the atoms
are relatively far away from each other.

Taking this into account, in this paper, we predict a mirror-mediated,
targeted remote dipole-dipole interaction which can persist over distances
that are many orders of magnitude longer than the wavelength of the emitted
light. Suppose two atoms, $a$ and $b$, are separated by a two-sided partially-
transparent mirror and the position of each atom coincides with the position
of the mirror image of the atom on the opposite side. Moreover, a far-away
screen collects their spontaneously emitted photons, as illustrated in Fig.
1(b). Comparing Figs. 1(a) and (b) and viewing the atoms again as radiating
dipole antennae, we see that—for half of the emitted light—the two paths from
a source to a certain point on the far-away screen are always of the same
length. The resulting interference pattern is therefore the same as in the
above two-atom double-slit experiment [1], apart from a reduction in
visibility.

Since atomic dipole-dipole interactions are the result of interference effects
and the interference of the spontaneously emitted photons are the same for two
nearby atoms and two atoms on opposite sides of a partially-reflecting mirror,
the above discussion suggests an additional maximum of the dipole-dipole
interaction between two atoms on opposite sides of a partially-transparent
mirror. As we shall see below, the strength of this ultralong-range
interaction does not depend on the actual distance of the atoms but on the
distance between atom $a$ and the mirror image of atom $b$. Unfortunately, for
symmetric mirrors, constructive interference of light on one-side implies
destructive interference on the other [21]. Hence when interference effects
due to the reflection of light on one side of the mirror increase atomic
interactions, reflection on the other side reduces them such that there is no
net enhancement. Remote dipole-dipole interactions therefore require the
presence of an asymmetric mirror interface. This can be realised with the help
of mirror coatings with different surface roughness on both sides. For
simplicity, we assume in the following that the surface of the mirror
interface is smooth on one side but rough on the other, as illustrated in Fig.
1(c). Such an interface can be realised, for example, by placing a very thin
layer of tiny metallic droplets, e.g. gold droplets, onto a glass surface with
some space between them. Although the atomic interactions which we predict
here are weaker than the standard dipole-dipole interactions of nearby atoms,
they are expected to find applications, for example, in quantum sensing.

Atomic dipole-dipole interactions have already been studied in different
environments but so far, they have always been relatively short-range [22,
23]. Theoretical and experimental studies usually consider atom-mirror
interactions [24, 25, 26, 27], interactions between atoms on the same side of
an interface [28, 29, 30], atoms separated by negative index metamaterials and
other thin films [31, 32, 33, 34, 35, 36, 37, 38, 39] and atoms near one-
dimensional nanofibers and wave guides [40, 41, 42, 43]. In this paper we take
an alternative theoretical approach and emphasise that, if a physical system
can be modelled by classical physics, then the expectation values of the
corresponding quantum model must evolve in exactly the same way. Otherwise, it
would be possible to contradict quantum physics with classical experiments. In
this paper, we show that this observation can be used to simplify the
derivations of quantum optical master equations, thereby enabling the
modelling of the experimental setup shown in Fig. 1(c). All we need to know
is, how classical optics models light propagation in the absence of the atoms
but in the presence of the asymmetric mirror interface [44, 45].

##  II Results

In free space, the complex electric field observable
$\boldsymbol{E}(\boldsymbol{r})$ at position $\boldsymbol{r}$ can be written
as a superposition of local complex electric field contributions
$\boldsymbol{E}_{\boldsymbol{k}\lambda}(\boldsymbol{r})$ of travelling waves
with polarisations $\lambda$, wave vectors $\boldsymbol{k}$ and directions of
propagation $\boldsymbol{s}=\boldsymbol{k}/\|\boldsymbol{k}\|$,

$\displaystyle\boldsymbol{E}(\boldsymbol{r})$ $\displaystyle=$
$\displaystyle\sum_{\lambda=1,2}\int{\rm
d}^{3}\boldsymbol{k}\,\boldsymbol{E}_{\boldsymbol{k}\lambda}(\boldsymbol{r})\,.$
(1)

Since all $(\boldsymbol{k},\lambda)$ components travel at the speed of light
$c$ in straight lines, we therefore know that [44, 45]

$\displaystyle U_{\rm
F}(t,0)\,\boldsymbol{E}_{\boldsymbol{k}\lambda}(\boldsymbol{r})\,|0_{\rm
F}\rangle=\boldsymbol{E}_{\boldsymbol{k}\lambda}(\boldsymbol{r}+\boldsymbol{s}ct)\,|0_{\rm
F}\rangle\,,$ (2)

where $|0_{\rm F}\rangle$ and $U_{\rm F}(t,0)$ denote the vacuum state and the
time evolution operator of the free field Hamiltonian $H_{\rm F}$. Creating a
local field excitation with wave vector $\boldsymbol{k}$ and polarisation
$\lambda$ at an initial time $t=0$ at a position $\boldsymbol{r}$ and evolving
the resulting state vector for a time $t$ results in the same final state as
creating the same field excitation at time $t$ at position
$\boldsymbol{r}+\boldsymbol{s}ct$.

Now suppose a partially-transparent mirror with (real) transmission and
reflection rates $t_{i}$ and $r_{i}$ (cf. Fig. 1(c)) is placed in the $x=0$
plane, reduces electric field amplitudes and changes directions of propagation
but does not alter the general shape of any incoming wave packets. In this
case, light still travels along straight lines most of the time. However, for
example, for a field excitation which has been created at the position of atom
$b$ at $t=0$ and which has been transmitted or reflected within the time
interval $(0,t)$, Eq. (1) needs to be replaced by

$\displaystyle U_{\rm
F}(t,0)\,\boldsymbol{E}_{\boldsymbol{k}\lambda}(\boldsymbol{r}_{b})\,|0_{\rm
F}\rangle$ $\displaystyle=$ $\displaystyle
r_{b}\,\boldsymbol{E}_{\tilde{\boldsymbol{k}}\lambda}(\tilde{\boldsymbol{r}}_{b}+\tilde{\boldsymbol{s}}ct)\,|0_{\rm
F}\rangle$ (3)
$\displaystyle+t_{b}\,\boldsymbol{E}_{\boldsymbol{k}\lambda}(\boldsymbol{r}_{b}+\boldsymbol{s}ct)\,|0_{\rm
F}\rangle\,.~{}~{}~{}~{}$

The tilde indicates that a minus sign has been added to the $x$ component of a
vector, thereby ensuring that electric field vectors are always orthogonal to
their direction of propagation. For simplicity, we assume here that
transmission and reflection rates do not depend on the angle of incidence and
the frequency of the incoming light. Such a dependence would alter the
strength of the predicted interactions but we expect that our results remain
valid, at least qualitatively.

In the following, we denote the electron charge and the position and the
dipole moment operator of atom $i$ by $e$, $\boldsymbol{r}_{i}$ and
$\boldsymbol{D}^{(i)}$, respectively. Then, within the dipole and the rotating
wave approximation, the interaction Hamiltonian between the atoms and the
surrounding free radiation field can be written as [18, 19, 20]

$\displaystyle H_{\rm AF}$ $\displaystyle=$ $\displaystyle
e\sum_{i=a,b}\boldsymbol{D}^{(i)\dagger}\cdot\boldsymbol{E}(\boldsymbol{r}_{i})^{\dagger}+{\rm
H.c.}$ (4)

In Methods, we analyse the dynamics generated by this Hamiltonian with the
help of a Dyson series expansion. As we shall see below, as long as we know
how $\boldsymbol{D}^{(i)}$ and $\boldsymbol{E}(\boldsymbol{r}_{i})$ evolve in
the absence of atom-field interactions, the dynamics of the two atoms in Fig.
1(c) can be analysed in a relatively straightforward way.

Quantum optical master equations describe the dynamics of atomic density
matrices $\rho_{\rm A}(t)$ on a coarse grained time scale $\Delta t$ which is
much larger than the inverse transition frequency $1/\omega_{0}$ but also much
smaller than the atomic lifetime $1/\Gamma_{\rm free}$ [6, 7, 18, 19, 20]. To
obtain them, we assume that the free radiation field is initially in its
vacuum state $|0_{\rm F}\rangle$, evolve the atoms and the field for a time
$\Delta t$ with the total Hamiltonian $H$ and follow this dynamics with a
measurement as to whether or not a photon has been emitted. Proceeding as
described in Methods, one can then show that the time derivative
$\dot{\rho}_{\rm A}$ equals

$\displaystyle\dot{\rho}_{\rm A}=-{{\rm i}\over\hbar}\left(H_{\rm
cond}\rho_{\rm A}-\rho_{\rm A}H_{\rm cond}^{\dagger}\right)+{\cal L}(\rho_{\rm
A})~{}~{}$ (5)

to a very good approximation. The reset operator ${\cal L}(\rho_{\rm A})$ and
the non-Hermitian Hamiltonian $H_{\rm cond}$ in this equation can be used to
analyse the dynamics of the two atoms in a time interval $(0,\Delta t)$ under
the condition of an emission and of no emission, respectively.

Figure 2: Dependence of the real part of $\Gamma_{\rm mir}^{(ab)}$ in Eq. (7)
as a function of the effective distance $\xi$ between atom $a$ and the mirror
image of atom $b$ for different orientations of the atomic dipole moment
vectors and for $t_{a}r_{b}=0.5$. This constant represents the changes to the
spontaneous decay rates of certain collective atomic states in the presence of
a mirror-mediated atomic dipole-dipole interaction. For simplicity, we assume
here that
${\boldsymbol{D}}^{(a)}_{12}={\boldsymbol{D}}^{(b)}_{12}={\boldsymbol{D}}_{12}$
and that ${\boldsymbol{D}}_{12}$ is a real vector.

Taking into account that light that has been emitted by one atom may travel to
the opposite side of the mirror interface where it interferes with the
reflected light originating from the other atom, one can show that the
operators $\mathcal{L}(\rho_{\rm A})$ and $H_{\rm cond}$ are of the form

$\displaystyle\mathcal{L}(\rho_{\rm A})$ $\displaystyle=$
$\displaystyle\sum_{i,j=a,b}{\rm Re}\left(\Gamma_{\rm
mir}^{(ij)}\right)\,\sigma_{i}^{-}\rho_{\rm A}\sigma_{j}^{+}\,,$
$\displaystyle H_{\rm cond}$ $\displaystyle=$ $\displaystyle H_{\rm
A}-\frac{{\rm i}\hbar}{2}\sum_{i,j=a,b}\Gamma_{\rm
mir}^{(ij)}\,\sigma_{j}^{+}\sigma_{i}^{-}\,.$ (6)

As we shall see below, the complex constants $\Gamma_{\rm mir}^{(ij)}$ in this
equation depend on the properties of the atoms and on the reflection and
transmission rates of the mirror interface. Here we are especially interested
in the case where the position $\boldsymbol{r}_{a}=(x_{a},y_{a},z_{a})$ of
atom $a$ and the position $\tilde{\boldsymbol{r}}_{b}=(-x_{b},y_{b},z_{b})$ of
the mirror image of atom $b$ are almost the same. For simplicity, let us
assume that $y_{a}=y_{b}$ and $z_{a}=z_{b}$ such that the relative effective
distance between $\boldsymbol{r}_{a}$ and $\tilde{\boldsymbol{r}}_{b}$ becomes
$\xi=k_{0}(x_{a}+x_{b})$. Using this notation and considering real
transmission and reflection rates [46], one can show that $\Gamma^{(aa)}_{\rm
mir}=\Gamma^{(bb)}_{\rm mir}=\Gamma_{\rm free}$ and $\Gamma_{\rm
mir}^{(ba)}=\Gamma_{\rm mir}^{(ab)*}$ with

$\displaystyle{\rm Re}\left(\Gamma_{\rm mir}^{(ab)}\right)$ $\displaystyle=$
$\displaystyle-{3\over 8}t_{a}r_{b}\Gamma_{\rm
free}\left[\left(\hat{\boldsymbol{D}}^{(a)}_{12}\cdot\hat{\boldsymbol{x}}\right)\left(\hat{\boldsymbol{D}}^{(b)}_{12}\cdot\hat{\boldsymbol{x}}\right)\right.$
(7)
$\displaystyle\times\left({\sin\xi\over\xi}+{3\cos\xi\over\xi^{2}}-{3\sin\xi\over\xi^{3}}\right)$
$\displaystyle\left.-\hat{\boldsymbol{D}}^{(a)}_{12}\cdot\hat{\boldsymbol{D}}^{(b)}_{12}\left({\sin\xi\over\xi}+{\cos\xi\over\xi^{2}}-{\sin\xi\over\xi^{3}}\right)\right]~{}~{}~{}$

for the situation depicted in Fig. 1(c). Here $\Gamma_{\rm free}$ denotes the
single-atom free space decay rate and $\hat{\boldsymbol{D}}^{(i)}_{12}$ and
$\hat{\boldsymbol{x}}$ are unit vectors in the direction of the (real) dipole
moment vector $\boldsymbol{D}^{(i)}_{12}$ of atom $i$ and in the direction of
the positive $x$ axis. Fig. 2 shows the real part of $\Gamma_{\rm mir}^{(ab)}$
for different orientations of the atomic dipole moments and for distances
$\|\boldsymbol{r}_{a}-\tilde{\boldsymbol{r}}_{b}\|$ of the order of the
wavelength of the emitted light. In the absence of a mirror interface,
$r_{b}=0$ and $\Gamma_{\rm mir}^{(ab)}$ tends to zero, as one would expect.

##  III Discussion

Having a closer look at the conditional Hamiltonian $H_{\rm cond}$ in Eq.
(II), we see that it contains a Hermitian and a non-Hermitian contribution.
The Hermitian contribution contains $H_{\rm A}$ and terms proportional to the
imaginary parts of $\Gamma_{\rm mir}^{(ab)}$ and $\Gamma_{\rm mir}^{(ba)}$
which describe the free dynamics of the atoms as well as interaction-induced
level shifts. As one can see from Eq. (35), these level shifts are essentially
the same as the level shifts cause by atomic dipole-dipole interactions in
free space [17, 18, 19, 20]. Apart from some additional terms, the only
difference is their dependence on the distance between atom $a$ and the mirror
image of atom $b$ instead of depending on the actual distance of the atoms.
This is why we refer to the interaction which we derive in this paper as a
dipole-dipole interaction.

The remaining terms in Eq. (II) describe the damping of population in excited
atomic states. Diagonalising the conditional Hamiltonian $H_{\rm cond}$, we
find that the spontaneous decay rate of the double-excited atomic state
$|22\rangle$ equals $2\Gamma_{\rm free}$, as usual. However, collective atomic
states which share only one excitation now have the spontaneous decay rates

$\displaystyle\Gamma_{\pm}$ $\displaystyle=$ $\displaystyle\Gamma_{\rm
free}\pm{\rm Re}\left(\Gamma^{(ab)}_{\rm mir}\right)$ (8)

which is different from their interaction-free spontaneous decay rate
$\Gamma_{\rm free}$. As we can see from Eq. (7) and Fig. 2, these changes are
analogous to the changes of spontaneous decay rates of two atoms with free-
space dipole-dipole interactions when the effective distance $\xi$ is
comparable to the wavelength of the emitted light. The atoms evolve as one
would expect in the presence of dipole-dipole interactions between atom $a$
and the mirror image of atom $b$ [17, 18, 19, 20].

Changes to spontaneous decay rates can be detected, for example, with the help
of fluorescence lifetime measurements. Moreover, when the atoms are driven by
a common laser field, we expect higher order photon correlation functions [8,
9, 10] to change and an interference pattern to emerge, if the spontaneously
emitted photons are collected on a far-away screen, as illustrated in Fig.
1(b). As described in Methods, the only assumptions regarding the distance of
the two atoms made in the derivation of the above equations are:

  1. 1.

The actual distance between the atoms should be relatively large, such that
$k_{0}|x_{a}-x_{b}|\gg 1$. Moreover, the time it takes light to travel from an
atom to the mirror should be much shorter than $\Delta t$.

  2. 2.

The actual distance between the two atoms in Fig. 1(c) should not be so large
that the time it takes light to travel from one atom to the other becomes
comparable to the lifetime of excited atomic states.

The first condition allows us to ignore direct atom-atom and atom-mirror
interactions which are relatively short-range. The second condition simplifies
the modelling of light propagation in the presence of the mirror interface and
is not very restrictive. For example, light can travel a $1\,$mm distance in
less than $3.4\cdot 10^{-12}\,$s which is much shorter than the typical
lifetime $1/\Gamma_{\rm free}$ of excited atomic states. However, when
analysing atomic interactions over very large distances, retardation effects
need to be taken into account and the dynamics of the two atoms can no longer
be described by a simple Markovian master equation, like the one in Eq. (5).

The main difference between the above-described and the usual dipole-dipole
interactions between two atoms in free space is that the former can be felt
over much longer distances. As we have shown above, mirror-mediated remote
dipole-dipole interaction can persist over distances which are several orders
of magnitude longer than the wavelength of the emitted light. They assume a
maximum when the effective relative distance $\xi$ between the position
$\boldsymbol{r}_{a}$ of atom $a$ and the position $\tilde{\boldsymbol{r}}_{b}$
of the mirror image of atom $b$ is of the order of one. This is the case even
when the actual distance of atom $a$ and atom $b$, i.e.
$\|\boldsymbol{r}_{a}-\boldsymbol{r}_{b}\|$, is several orders of magnitude
larger than the wavelength of the emitted light. The interaction which we
describe in this paper is therefore ultralong-range but also highly-targeted.

Another requirement for the mirror-mediated atomic interactions which we
predict in this paper is the presence of an asymmetric interface. As
illustrated in Fig. 1(c), we assume different surface roughness on both sides
of the thin reflecting layer. If both sides of the mirror surface would be
smooth and reflect light like an ideal mirror, the interaction constant
$\Gamma_{\rm mir}^{(ab)}$ in Eq. (35) would be proportional to
$r_{a}t_{b}+t_{a}r_{b}$. Unfortunately, we know from classical optics (from
Stoke’s law) that this factor equals zero for symmetric mirrors [46] and the
interaction which we predict here simply disappears in this case. However, in
the case of different surface roughness, Stokes theorem no longer applies and
constructive interference on one side of the mirror no longer implies the same
amount of destructive interference on the other.

In conclusion, in this paper we derived the quantum optical master equations
in Eq. (5) with the help of the usual dipole and the rotating wave
approximations. In addition, we evolved the atoms and the free radiation field
for a short time interval $\Delta t$ using a Dyson series expansion. This
approach allowed us to deduce the time derivative of the atomic density matrix
$\dot{\rho}_{\rm A}$ from the classical dynamics of light in the presence of a
two-sided partially-transparent mirror but in the absence of any atom-field
interactions. We then showed that the spontaneous emission rates of two atoms
on opposite sides of a partially-transparent asymmetric mirror interface can
produce an effective dipole-dipole interaction between atom $a$ and the mirror
image of atom $b$ and vice versa. The main result of this paper is the
prediction of a targeted, remote, mirror-mediated ultralong-range dipole-
dipole interaction which is likely to find a wide range of applications in the
design of novel photonic devices for quantum technology applications,
especially in quantum sensing.

##  IV Methods

###  IV.1 The derivation of quantum optical master equations

Our starting point for the derivation of the master equations in Eq. (5) of
the two atoms in Fig. 1(c) is their system Hamiltonian $H$. This Hamiltonian
can be written as $H_{\rm A}+H_{\rm F}+H_{\rm AF}$ where $H_{\rm A}$ and
$H_{\rm F}$ denote the free energy of the atoms and of the electromagnetic
field in the presence of the mirror interface. An expression for the
interaction Hamiltonian $H_{\rm AF}$ between the atoms and the surrounding
free radiation field within the usual dipole approximation [18, 19, 20] can be
found in Eq. (4). As we shall see below, in addition, we only need to know how
the atomic dipole moments and electric field observables evolve in the
Heisenberg picture in the absence of any atom-field interactions.

Suppose $\rho_{\rm A}(0)$ is the initial density matrix of the two atoms in
the Schrödinger picture, while the surrounding free radiation field is
initially in its vacuum state, $|0_{\rm F}\rangle$. We then evolve the atom-
field density matrix $|0_{\rm F}\rangle\rho_{\rm A}(0)\langle 0_{\rm F}|$ for
a time $\Delta t$ with the time evolution operator $U(\Delta t,0)$ of the
total system. Subsequently performing an absorbing measurement on the
surrounding free radiation field leads to the atomic density matrix

$\displaystyle\rho_{\rm A}(\Delta t)$ $\displaystyle=$ $\displaystyle{\rm
Tr}_{\rm F}\big{[}U(\Delta t,0)|0_{\rm F}\rangle\rho_{\rm A}(0)\langle 0_{\rm
F}|U^{\dagger}(\Delta t,0)\big{]}\,,~{}~{}~{}$ (9)

while placing the field again in its vacuum state. Here the trace over the
field degrees of freedom is taken to ensure that a measurement on the
surrounding electromagnetic field does not change the properties of the atomic
degrees of freedom. As requested by locality, the density matrices $\rho_{\rm
A}(\Delta t)$ and $U(\Delta t,0)|0_{\rm F}\rangle\rho_{\rm A}(0)\langle 0_{\rm
F}|U^{\dagger}(\Delta t,0)$ must have the same atomic expectation values, if
the outcome of the photon measurement is ignored [6, 7]. Next, we introduce
the time derivative

$\displaystyle\dot{\rho}_{\rm A}$ $\displaystyle=$ $\displaystyle{1\over\Delta
t}\left(\rho_{\rm A}(\Delta t)-\rho_{\rm A}(0)\right)$ (10)

which describes the dynamics of the atomic density ${\rho}_{\rm A}$ on the
coarse grained time scale $\Delta t$.

Since the time evolution operator $U(\Delta t,0)$ in Eq. (9) cannot be
calculated easily analytically, we write the total Hamiltonian $H$ of the
experimental setup in Fig. 1(c) in the following as the sum of the free
Hamiltonian $H_{0}=H_{\rm A}+H_{\rm F}$ and the interaction $H_{\rm AF}$. As
long as $\Delta t$ is neither too long nor too short, as described in Results,
we can analyse the dynamics of the system using a Dyson series expansion which
implies that

$\displaystyle U(\Delta t,0)$ $\displaystyle=$ $\displaystyle U_{0}(\Delta
t,0)-{{\rm i}\over\hbar}\int_{0}^{\Delta t}{\rm d}t\,U_{0}(\Delta t,t)\,H_{\rm
AF}\,U_{0}(t,0)-{1\over\hbar^{2}}\int_{0}^{\Delta t}{\rm d}t\int_{0}^{t}{\rm
d}t^{\prime}\,U_{0}(\Delta t,t)\,H_{\rm AF}\,U_{0}(t,t^{\prime})\,H_{\rm
AF}\,U_{0}(t^{\prime},0)$

to a very good approximation. Combining Eqs. (9) and (IV.1), while only taking
terms in zeroth order in $\Delta t$ into account, therefore leads to

$\displaystyle\rho_{\rm A}(\Delta t)$ $\displaystyle=$
$\displaystyle\frac{1}{\hbar^{2}}\int_{0}^{\Delta t}{\rm d}t\int_{0}^{\Delta
t}{\rm d}t^{\prime}\,{\rm Tr}_{\rm F}\left[U_{0}(\Delta t,t)\,H_{\rm
AF}\,U_{0}(t,0)\,|0_{\rm F}\rangle\rho_{\rm A}(0)\langle 0_{\rm
F}|\,U^{\dagger}_{0}(t^{\prime},0)\,H_{\rm AF}\,U^{\dagger}_{0}(\Delta
t,t^{\prime})\right]$ (12)
$\displaystyle-\frac{1}{\hbar^{2}}\int\limits_{0}^{\Delta t}{\rm
d}t\int\limits_{0}^{t}{\rm d}t^{\prime}\,\langle 0_{\rm F}|U_{0}(\Delta
t,t)\,H_{\rm AF}\,U_{0}(t,t^{\prime})\,H_{\rm AF}\,U_{0}(t^{\prime},0)|0_{\rm
F}\rangle\rho_{\rm A}(0)+{\rm c.c.}$ $\displaystyle+\langle 0_{\rm
F}|U_{0}(\Delta t,0)|0_{\rm F}\rangle\rho_{\rm A}(0)\langle 0_{\rm
F}|U^{\dagger}_{0}(\Delta t,0)|0_{\rm F}\rangle$

in first order in $\Delta t$. To obtain the above equation, we took into
account that $H_{\rm AF}$ either creates or annihilates a photon, while
$H_{0}$ preserves the number of excitations in the free radiation field.
Comparing this equation with Eq. (10), we see that Eq. (5) applies, if we
define

$\displaystyle{\cal L}(\rho_{\rm A})$ $\displaystyle=$
$\displaystyle\frac{1}{\hbar^{2}\Delta t}\int_{0}^{\Delta t}{\rm
d}t\int_{0}^{\Delta t}{\rm d}t^{\prime}\,{\rm Tr}_{\rm F}\left[U_{0}(\Delta
t,t)\,H_{\rm AF}\,U_{0}(t,0)\,|0_{\rm F}\rangle\rho_{\rm A}\langle 0_{\rm
F}|\,U^{\dagger}_{0}(t^{\prime},0)\,H_{\rm AF}\,U^{\dagger}_{0}(\Delta
t,t^{\prime})\right]\,,$ $\displaystyle H_{\rm cond}$ $\displaystyle=$
$\displaystyle H_{\rm A}-\frac{{\rm i}}{\hbar\Delta t}\int\limits_{0}^{\Delta
t}{\rm d}t\int\limits_{0}^{t}{\rm d}t^{\prime}\,\langle 0_{\rm F}|U_{0}(\Delta
t,t)\,H_{\rm AF}\,U_{0}(t,t^{\prime})\,H_{\rm AF}\,U_{0}(t^{\prime},0)|0_{\rm
F}\rangle\,.$ (13)

Here ${\cal L}(\rho_{\rm A})$ is the (unnormalised) density matrix of the
atoms conditional on the creation of a photon in $(0,\Delta t)$. Analogously,
the Hamiltonian $H_{\rm cond}$ describes their dynamics in the absence of a
photon [6, 7].

To further simplify the above expressions, we notice that $H_{0}=H_{\rm
A}+H_{\rm F}$ is the sum of two commuting Hamiltonians. Hence, $U_{0}=U_{\rm
A}\otimes U_{\rm F}$, where $U_{\rm A}$ and $U_{\rm F}$ denote the time
evolution operators associated with $H_{\rm A}$ and with $H_{\rm F}$,
respectively. In addition, we introduce the short hand notation

$\displaystyle\boldsymbol{D}^{(i)}(t)$ $\displaystyle=$ $\displaystyle
U^{\dagger}_{\rm A}(t,0)\,\boldsymbol{D}^{(i)}\,U_{\rm A}(t,0)$ (14)

and notice that the vacuum state is invariant under $U_{\rm F}$. Using Eqs.
(4) and (IV.1), one can then show that

$\displaystyle{\cal L}(\rho_{\rm A})$ $\displaystyle=$
$\displaystyle\frac{e^{2}}{\hbar^{2}\Delta t}\,\sum_{i,j=a,b}\int_{0}^{\Delta
t}{\rm d}t\int_{0}^{\Delta t}{\rm d}t^{\prime}\,{\rm Tr}_{\rm
F}\left[\boldsymbol{D}^{(i)}(t)\cdot U_{\rm F}(\Delta
t,t)\boldsymbol{E}(\boldsymbol{r}_{i})\,|0_{\rm F}\rangle\rho_{\rm A}\langle
0_{\rm
F}|\,\boldsymbol{D}^{(j)}(t^{\prime})^{\dagger}\cdot\boldsymbol{E}(\boldsymbol{r}_{j})^{\dagger}\,U^{\dagger}_{\rm
F}(\Delta t,t^{\prime})\right],$ $\displaystyle H_{\rm cond}$ $\displaystyle=$
$\displaystyle H_{\rm A}-\frac{{\rm i}e^{2}}{\hbar\Delta
t}\,\sum_{i,j=a,b}\int\limits_{0}^{\Delta t}{\rm d}t\int\limits_{0}^{t}{\rm
d}t^{\prime}\,\langle 0_{\rm
F}|\boldsymbol{D}^{(j)}(t)^{\dagger}\cdot\boldsymbol{E}(\boldsymbol{r}_{j})^{\dagger}\,U^{\dagger}_{\rm
F}(t^{\prime},0)\,\boldsymbol{D}^{(i)}(t^{\prime})\cdot U_{\rm
F}(t,0)\boldsymbol{E}(\boldsymbol{r}_{i})\,|0_{\rm F}\rangle$ (15)

in zeroth order in $\Delta t$. As one would expect, ${\cal L}(\rho_{\rm A})$
contains all the contributions of the atom-field density matrix which
correspond to the presence of a photon at $\Delta t$ in the free radiation
field. Moreover, the no-photon Hamiltonian contains contributions in which
excitation has been created within $(0,\Delta t)$ but is later re-absorbed by
the atoms.

####  IV.1.1 The free-space dynamics of atoms and field

Figure 3: Spacetime diagram for light propagating in the presence of a
partially-transparent mirror in the $x=0$ plane from the position of atom $a$
(green) and from the position of atom $b$ (blue), respectively, for some time
$t$ along the $x$ axis. Here we assume that the mirror interface (which should
be much thinner than the wavelength of the emitted light) is very smooth on
the left but rough on the right has left, thereby randomly altering the wave
vector of the reflected light originating from the position of atom $b$.
Suppose the time $t$ is larger than the time it takes to travel from an atom
to the mirror interface. In this case, the light arrives at time $t$ at the
indicated positions. Analogous observations in three dimensions can be used to
derive Eq. (IV.1.1).

Suppose $\hbar\omega_{0}$ is the energy gap between the ground and the excited
state of atom $i$, i.e. between $|1\rangle_{i}$ and $|2\rangle_{i}$. In
addition, we introduce the atomic raising and the lowering operators
$\sigma^{+}_{i}=|2\rangle_{ii}\langle 1|$ and
$\sigma^{-}_{i}=|1\rangle_{ii}\langle 2|$ and the complex dipole moment
operators $\boldsymbol{D}^{(i)}=\boldsymbol{D}_{12}^{(i)}\,\sigma^{-}_{i}$.
Using this notation, the atom Hamiltonian $H_{\rm A}$ can be written as

$\displaystyle H_{\rm A}$ $\displaystyle=$
$\displaystyle\sum_{i=a,b}\hbar\omega_{0}\,\sigma^{+}_{i}\sigma^{-}_{i}\,.$
(16)

Hence the time-dependent dipole moment operator $\boldsymbol{D}^{(i)}(t)$ in
Eq. (14) equals

$\displaystyle\boldsymbol{D}^{(i)}(t)$ $\displaystyle=$ $\displaystyle{\rm
e}^{-{\rm i}\omega_{0}t}\,\boldsymbol{D}_{12}^{(i)}\,\sigma^{-}_{i}\,.$ (17)

As we have seen above, this operator is the only atomic operator needed for
the derivation of Eq. (5). A closer look at Eq. (IV.1) shows that the only
other expressions that we need to calculate are states of the free radiation
field of the form $U_{\rm F}(t,0)\,\boldsymbol{E}(\boldsymbol{r})\,|0_{\rm
F}\rangle$ with $t\in(0,\Delta t)$. Fortunately, this state which is the state
obtained when creating a field excitation at time $t=0$ at position
$\boldsymbol{r}$ and then letting this radiation propagate for some time $t$
can be deduced relatively easily from classical electrodynamics. To do so, we
write $\boldsymbol{E}(\boldsymbol{r})$ as in Eq. (1) as a superposition of
field contributions $\boldsymbol{E}_{\boldsymbol{k}\lambda}(\boldsymbol{r})$
with well-defined directions of propagation
$\boldsymbol{s}=\boldsymbol{k}/\|\boldsymbol{k}\|$ and assume that [47]

$\displaystyle\boldsymbol{E}_{\boldsymbol{k}\lambda}(\boldsymbol{r})$
$\displaystyle=$ $\displaystyle-{\rm i}\left({\hbar ck\over
16\pi^{3}\varepsilon}\right)^{1/2}{\rm e}^{-{\rm
i}{\boldsymbol{k}}\cdot{\boldsymbol{r}}}\,a^{\dagger}_{{\boldsymbol{k}}\lambda}\,\boldsymbol{e}_{\boldsymbol{s}\lambda}~{}~{}$
(18)

with $k=\|\boldsymbol{k}\|$. The vectors
$\boldsymbol{e}_{{\boldsymbol{s}}\lambda}$ in this equation are normalised
polarisation vectors with
$\boldsymbol{e}_{\boldsymbol{s}1}\cdot\boldsymbol{e}_{\boldsymbol{s}2}=\boldsymbol{e}_{\boldsymbol{s}\lambda}\cdot\boldsymbol{s}=0$.
The complex electric field observables
$\boldsymbol{E}_{\boldsymbol{k}\lambda}(\boldsymbol{r})$ are functions of
bosonic annihilation operators $a_{{\boldsymbol{k}}\lambda}$ of free-space
photons with wave vectors ${\boldsymbol{k}}=k{\boldsymbol{s}}$ and
polarisations $\lambda$.

Fig. 3 illustrates the dynamics of light with the help of a space time diagram
for the case, where a field excitation is created at the position of atom $a$
(green) and atom $b$ (blue), respectively, at time $t=0$ and propagates until
a time $t$ in the $y$-$z$ plane along different light cones. Since $\Delta t$
is much larger than the time it takes light to travel from the position
$\boldsymbol{r}_{i}$ of atom $i$ towards the mirror surface, the same applies
to most of the times $t\in(0,\Delta t)$. This observation allows us to neglect
small times $t$ for which light that is eventually transmitted or reflected
has not yet reached the interface. Using the same arguments as in the
derivation of Eqs. (2) and (3), we therefore find that

$\displaystyle U_{\rm
F}(t,0)\,\boldsymbol{E}_{\boldsymbol{k}\lambda}(\boldsymbol{r}_{a})\,|0_{\rm
F}\rangle$ $\displaystyle=$
$\displaystyle\Theta(-s_{x})\left[r_{a}\,\boldsymbol{E}_{{\boldsymbol{R}}({\boldsymbol{k}})\lambda}({\boldsymbol{S}}({\boldsymbol{r}}_{a},t))+t_{a}\,\boldsymbol{E}_{\boldsymbol{k}\lambda}(\boldsymbol{r}_{a}+\boldsymbol{s}ct)\right]|0_{\rm
F}\rangle+\Theta(s_{x})\,\boldsymbol{E}_{\boldsymbol{k}\lambda}(\boldsymbol{r}_{a}+\boldsymbol{s}ct)\,|0_{\rm
F}\rangle\,,$ $\displaystyle U_{\rm
F}(t,0)\,\boldsymbol{E}_{\boldsymbol{k}\lambda}(\boldsymbol{r}_{b})\,|0_{\rm
F}\rangle$ $\displaystyle=$
$\displaystyle\Theta(s_{x})\left[r_{b}\,\boldsymbol{E}_{\tilde{\boldsymbol{k}}\lambda}(\tilde{\boldsymbol{r}}_{b}+\tilde{\boldsymbol{s}}ct)+t_{b}\,\boldsymbol{E}_{\boldsymbol{k}\lambda}(\boldsymbol{r}_{b}+\boldsymbol{s}ct)\right]|0_{\rm
F}\rangle+\Theta(-s_{x})\,\boldsymbol{E}_{\boldsymbol{k}\lambda}(\boldsymbol{r}_{b}+\boldsymbol{s}ct)\,|0_{\rm
F}\rangle$ (19)

for direction vectors $\boldsymbol{s}=(s_{x},s_{y},s_{z})$. Here the Heavyside
function $\Theta(x)$ equals 0 for $x<0$ and 1 otherwise, while
$\boldsymbol{r}=(x,y,z)$. Moreover, ${\boldsymbol{R}}({\boldsymbol{k}})$ and
${\boldsymbol{S}}({\boldsymbol{r}}_{a},t)$ denote the wave vector and the
position of a local electric field excitation at time $t$ which has been
created by at the position of atom $a$ at $t=0$ after reflection by the rough
side of the interface (cf. Fig. 1(c)).

####  IV.1.2 The conditional Hamiltonian $H_{\rm cond}$

Substituting Eq. (17) into Eq. (IV.1), we can show that the conditional
Hamiltonian $H_{\rm cond}$ can indeed be written as in Eq. (II), if we define
the constants $\Gamma_{\rm mir}^{(ij)}$ such that

$\displaystyle\Gamma_{\rm mir}^{(ij)}$ $\displaystyle=$
$\displaystyle{1\over\Delta t}\int_{0}^{\Delta t}{\rm d}t\int_{0}^{t}{\rm
d}t^{\prime}\,{e^{2}c\over 8\hbar\pi^{3}\varepsilon}\,{\rm e}^{{\rm
i}\omega_{0}(t-t^{\prime})}\,\gamma_{\rm mir}^{(ij)}(t,t^{\prime})$ (20)

with $\gamma_{\rm mir}^{(ij)}(t,t^{\prime})$ given by

$\displaystyle\gamma_{\rm mir}^{(ij)}(t,t^{\prime})$ $\displaystyle=$
$\displaystyle{16\pi^{3}\varepsilon\over\hbar c}\,\langle 0_{\rm
F}|\boldsymbol{D}_{12}^{(j)}\cdot\boldsymbol{E}(\boldsymbol{r}_{j})^{\dagger}\,U^{\dagger}_{\rm
F}(t^{\prime},0)\,\boldsymbol{D}^{(i)}_{12}\cdot U_{\rm
F}(t,0)\boldsymbol{E}(\boldsymbol{r}_{i})|0_{\rm F}\rangle\,.$ (21)

As we shall see below, only constants $\gamma_{\rm mir}^{(ij)}$ with
$t=t^{\prime}$ contribute to the no-photon dynamics of the atoms. This means,
the no-photon dynamics is only due to trajectories of the atom-field system
where a photon is created at a certain time $t$ that never reaches a detector,
since it gets immediately re-absorbed. More concretely, using Eqs. (1), (18)
and (IV.1.1) and performing one of the ${\boldsymbol{k}}$ integrations, we
find for example that

$\displaystyle\gamma_{\rm mir}^{(aa)}(t,t^{\prime})$ $\displaystyle=$
$\displaystyle\sum_{\lambda=1,2}\int{\rm d}^{3}\boldsymbol{k}\,k\,{\rm
e}^{-{\rm
i}ck(t-t^{\prime})}\left[\left(t_{a}^{2}\,\Theta(-s_{x})+\Theta(s_{x})\right)\big{(}\boldsymbol{D}^{(a)}_{12}\cdot\boldsymbol{e}_{\boldsymbol{s}\lambda}\big{)}^{2}+r_{a}^{2}\,\Theta(-s_{x})\,\big{(}\boldsymbol{D}^{(a)}_{12}\cdot\boldsymbol{e}_{R({\boldsymbol{s}})\lambda}\big{)}^{2}\right]$
$\displaystyle\gamma_{\rm mir}^{(bb)}(t,t^{\prime})$ $\displaystyle=$
$\displaystyle\sum_{\lambda=1,2}\int{\rm d}^{3}\boldsymbol{k}\,k\,{\rm
e}^{-{\rm
i}ck(t-t^{\prime})}\left[\left(\Theta(-s_{x})+t_{b}^{2}\,\Theta(s_{x})\right)\big{(}\boldsymbol{D}^{(b)}_{12}\cdot\boldsymbol{e}_{\boldsymbol{s}\lambda}\big{)}^{2}+r_{b}^{2}\,\Theta(s_{x})\,\big{(}\boldsymbol{D}^{(b)}_{12}\cdot\boldsymbol{e}_{\tilde{\boldsymbol{s}}\lambda}\big{)}^{2}\right].$
(22)

This equation does not contain any terms that are proportional to $r_{a}t_{a}$
and $r_{b}t_{b}$, since the wave vectors associated with reflection and
transmission are always orthogonal. To calculate $\Gamma_{\rm mir}^{(aa)}$ and
$\Gamma_{\rm mir}^{(bb)}$, we now introduce polar coordinates
$k\in(0,\infty)$, $\varphi\in(0,2\pi)$ and $\vartheta\in(0,\pi)$ such that

$\displaystyle\boldsymbol{s}=\left(\begin{array}[]{c}\cos\vartheta\\\
\cos\varphi\,\sin\vartheta\\\
\sin\varphi\,\sin\vartheta\end{array}\right),~{}~{}~{}\boldsymbol{e}_{\boldsymbol{s}1}=\begin{pmatrix}0\\\
\sin\varphi\\\
-\cos\varphi\end{pmatrix},~{}~{}~{}\boldsymbol{e}_{\boldsymbol{s}2}=\begin{pmatrix}\sin\vartheta\\\
-\cos\varphi\cos\vartheta\\\ -\sin\varphi\cos\vartheta\end{pmatrix}$ (26)

and replace ${\rm d}^{3}\boldsymbol{k}$ by ${\rm d}^{2}\boldsymbol{s}\,{\rm
d}k\,k^{2}$ with ${\rm d}^{2}\boldsymbol{s}={\rm d}\vartheta\,{\rm
d}\varphi\,\sin\vartheta$. Hence Eq. (IV.1.2) contains integrals of the form

$\displaystyle\int_{0}^{\infty}{\rm d}k\,k^{3}\,{\rm e}^{-{\rm i}ck\tau}$
$\displaystyle=$ $\displaystyle-{{\rm i}\over c^{3}}\,{{\rm d}^{3}\over{\rm
d}\tau^{3}}\,\int_{0}^{\infty}{\rm d}k\,{\rm e}^{-{\rm i}ck\tau}=-{{\rm
i}\pi\over c^{4}}\,{{\rm d}^{3}\over{\rm d}\tau^{3}}\,\delta(\tau)=-{{\rm
i}\pi\over c^{4}}\,\delta^{(3)}(\tau)$ (27)

with $\delta^{(3)}(\tau)$ denoting the third derivative of $\delta(\tau)$ with
respect to $\tau$. As usual, we here neglect an infinitely large imaginary
part which describes self-interactions. When avoiding approximations, these
become finite and can be absorbed into the free Hamiltonian $H_{\rm A}$ of the
atoms. Using Eq. (27), we can now show that

$\displaystyle{1\over\Delta t}\int_{0}^{\Delta t}{\rm d}t\int_{0}^{t}{\rm
d}t^{\prime}\,{\rm e}^{{\rm i}\omega_{0}(t-t^{\prime})}\int_{0}^{\infty}{\rm
d}k\,k^{3}\,{\rm e}^{-{\rm i}ck(t-t^{\prime})}$ $\displaystyle=$
$\displaystyle{1\over\Delta t}\int_{0}^{\Delta t}{\rm d}t\int_{0}^{t}{\rm
d}\tau\,{\rm e}^{{\rm i}\omega_{0}\tau}\int_{0}^{\infty}{\rm d}k\,k^{3}\,{\rm
e}^{-{\rm i}ck\tau}$ (28) $\displaystyle=$ $\displaystyle{{\rm i}\pi\over
c^{3}\Delta t}\int_{0}^{\Delta t}{\rm d}t\int_{0}^{t}{\rm d}\tau\,{\rm
e}^{{\rm i}\omega_{0}\tau}\,\delta^{(3)}(\tau)$ $\displaystyle=$
$\displaystyle{\pi\omega_{0}^{3}\over c^{4}}\,.$

Combining the above equations and assuming that there is no focussing and that
the vectors $R({\boldsymbol{s}})$ cover the wave vectors ${\boldsymbol{s}}$ on
the right hand side of the interface evenly, we then find that

$\displaystyle\Gamma_{\rm mir}^{(ii)}$ $\displaystyle=$
$\displaystyle{e^{2}\omega_{0}^{3}\over 8\pi^{2}\hbar\varepsilon
c^{3}}\sum_{\lambda=1,2}\int{\rm
d}^{2}\boldsymbol{s}\,\left[r_{i}^{2}\,\Theta(\mp
s_{x})\,\big{(}\boldsymbol{D}^{(i)}_{12}\cdot\boldsymbol{e}_{\tilde{\boldsymbol{s}}\lambda}\big{)}^{2}+\left(t_{i}^{2}\,\Theta(\mp
s_{x})+\Theta(\pm
s_{x})\right)\big{(}\boldsymbol{D}^{(i)}_{12}\cdot\boldsymbol{e}_{\boldsymbol{s}\lambda}\big{)}^{2}\right]\,.$
(29)

Introducing the notation
$\boldsymbol{D}^{(i)}_{12}=\|\boldsymbol{D}_{12}\|\,(d^{(i)}_{1},d^{(i)}_{2},d^{(i)}_{3})^{\rm
T}$ with $|d^{(i)}_{1}|^{2}+|d^{(i)}_{2}|^{2}+|d^{(i)}_{3}|^{2}=1$, one can
now show that the $\Gamma_{\rm mir}^{(ii)}$ both equal the free-space decay
rate of an atom with dipole moment
$\boldsymbol{D}^{(i)}_{12}=\boldsymbol{D}_{12}$,

$\displaystyle\Gamma_{\rm mir}^{(aa)}=\Gamma_{\rm mir}^{(bb)}$
$\displaystyle=$ $\displaystyle\Gamma_{\rm
free}={e^{2}\omega_{0}^{3}\,\|\boldsymbol{D}_{12}\|^{2}\over
3\pi\hbar\varepsilon c^{3}}\,,$ (30)

since $r_{i}^{2}+t_{i}^{2}=1$ [7]. The derivation of the two remaining
constants, $\Gamma_{\rm mir}^{(ab)}$ and $\Gamma_{\rm mir}^{(ba)}$, is
slightly more complicated. However, here we are only interested in the case
where the distance of the atoms from the mirror interface and the distance
between atom $a$ and atom $b$ are all much longer than the wavelength of the
emitted light. In this case, we can safely ignore terms which contain
exponential functions of
$\boldsymbol{k}\cdot(\boldsymbol{r}_{i}-\tilde{\boldsymbol{r}}_{i})$ and
$\boldsymbol{k}\cdot(\boldsymbol{r}_{a}-{\boldsymbol{r}}_{b})$ when evaluating
the above constants. These terms average away when one of the $\boldsymbol{k}$
integrations hidden in the electric field observables in Eq. (1) is performed,
since they describe interactions between an atom and its mirror image and the
usual dipole-dipole interaction between atom $a$ and $b$, respectively, which
are relatively short-range. Fortunately, exponentials of
$\boldsymbol{k}\cdot(\tilde{\boldsymbol{r}}_{a}-\boldsymbol{r}_{b})$ and
$\boldsymbol{k}\cdot(\boldsymbol{r}_{a}-\tilde{\boldsymbol{r}}_{b})$ must be
kept, when atom $a$ is placed near the position of the mirror image of atom
$b$ and $\boldsymbol{r}_{b}$ and $\tilde{\boldsymbol{r}}_{a}$ are close. In
this case

$\displaystyle\gamma_{\rm mir}^{(ab)}(t,t^{\prime})$ $\displaystyle=$
$\displaystyle\sum_{\lambda=1,2}\int{\rm
d}^{3}\boldsymbol{k}\,k\left[t_{a}r_{b}\,{\rm e}^{-{\rm
i}\boldsymbol{k}\cdot(\boldsymbol{r}_{a}+{\boldsymbol{s}}ct)}\,{\rm e}^{{\rm
i}\boldsymbol{k}\cdot(\tilde{\boldsymbol{r}}_{b}+{\boldsymbol{s}}ct^{\prime})}\,\big{(}\boldsymbol{D}^{(a)}_{12}\cdot\boldsymbol{e}_{\boldsymbol{s}\lambda}\big{)}\big{(}\boldsymbol{D}^{(b)}_{12}\cdot\boldsymbol{e}_{\boldsymbol{s}\lambda}\big{)}\right.$
(31) $\displaystyle\hskip 56.9055pt\left.+r_{a}t_{b}\,{\rm e}^{-{\rm
i}\boldsymbol{R}(\boldsymbol{k})\cdot\boldsymbol{S}(\boldsymbol{r}_{a},t^{\prime})}\,{\rm
e}^{{\rm
i}\boldsymbol{R}(\boldsymbol{k})\cdot(\boldsymbol{r}_{b}+{\boldsymbol{R}}(\boldsymbol{s})ct)}\big{(}\boldsymbol{D}^{(a)}_{12}\cdot\boldsymbol{e}_{\boldsymbol{R}(\boldsymbol{s})\lambda}\big{)}\big{(}\boldsymbol{D}^{(b)}_{12}\cdot\boldsymbol{e}_{\boldsymbol{R}(\boldsymbol{s})\lambda}\big{)}\right]\Theta(-s_{x})$

and $\gamma_{\rm mir}^{(ba)}(t,t^{\prime})=\gamma_{\rm
mir}^{(ab)}(t,t^{\prime})^{*}$, in analogy to Eq. (IV.1.2). If the atoms are
aligned such that ${\boldsymbol{r}}_{a}$ and ${\boldsymbol{r}}_{b}$ have the
same $y$ and $z$ coordinates, then

$\displaystyle{\rm e}^{-{\rm
i}\boldsymbol{k}\cdot(\boldsymbol{r}_{a}+{\boldsymbol{s}}ct)}\,{\rm e}^{{\rm
i}\boldsymbol{k}\cdot(\tilde{\boldsymbol{r}}_{b}+{\boldsymbol{s}}ct^{\prime})}$
$\displaystyle=$ $\displaystyle{\rm e}^{-{\rm i}ck(t-t^{\prime})}\,{\rm
e}^{-{\rm i}k\cos\vartheta(x_{a}+x_{b})}$ (32)

and the first term in the above equations can be further simplified relatively
easily. However, if the scattering operator ${\boldsymbol{S}}$ scrambles the
wave vectors of the reflected light randomly, the $r_{a}t_{b}$ term in Eq.
(31) averages away to zero and becomes negligible. In this case, the phase
difference between the reflected light originating from atom $a$ and the
transmitted light originating from atom $b$ differs for all directions. Using
again the above equations and taking into account that $s_{x}$ is only
negative when $\vartheta\in(\pi/2,\pi)$, we therefore find that

$\displaystyle\Gamma_{\rm mir}^{(ab)}$ $\displaystyle=$
$\displaystyle{3t_{a}r_{b}\Gamma_{\rm free}\over
8\pi}\sum_{\lambda=1,2}\int_{\pi/2}^{\pi}{\rm d}\vartheta\,{\rm e}^{-{\rm
i}\xi\cos\vartheta}\,\sin\vartheta\int_{0}^{2\pi}{\rm
d}\varphi\,\left(\hat{\boldsymbol{D}}^{(a)}_{12}\cdot\boldsymbol{e}_{\boldsymbol{s}\lambda}\right)\left(\hat{\boldsymbol{D}}^{(b)}_{12}\cdot\boldsymbol{e}_{\boldsymbol{s}\lambda}\right)\,.$
(33)

Here the hat symbols indicate that the vectors $\boldsymbol{D}^{(i)}_{12}$
have been normalised, $\xi=k_{0}(x_{a}+x_{b})$ with $k_{0}=\omega_{0}/c$ is an
effective distance and the polarisation vectors
$\boldsymbol{e}_{\boldsymbol{s}\lambda}$ can be found in Eq. (26). Hence
performing the $\varphi$ integration and substituting $u=-\cos\vartheta$
yields

$\displaystyle\Gamma_{\rm mir}^{(ab)}$ $\displaystyle=$
$\displaystyle{3t_{a}r_{b}\Gamma_{\rm free}\over 16}\int_{0}^{1}{\rm d}u\,{\rm
e}^{{\rm i}\xi
u}\left[2d_{1}^{(a)}d_{1}^{(b)}\left(1-u^{2}\right)+\left(d_{2}^{(a)}d_{2}^{(b)}+d_{3}^{(a)}d_{3}^{(b)}\right)\left(1+u^{2}\right)\right].$
(34)

Performing the final integration, the above constant simplifies to

$\displaystyle\Gamma_{\rm mir}^{(ab)}$ $\displaystyle=$
$\displaystyle{3t_{a}r_{b}\Gamma_{\rm free}\over
16}\left[\sum_{i=2}^{3}d_{i}^{(a)}d_{i}^{(b)}\left(2{\rm e}^{{\rm
i}\xi}\left({1\over{\rm i}\xi}+{1\over\xi^{2}}-{1\over{\rm
i}\xi^{3}}\right)-{1\over{\rm i}\xi}+{2\over{\rm
i}\xi^{3}}\right)-2d_{1}^{(a)}d_{1}^{(b)}\left(2{\rm e}^{{\rm
i}\xi}\left({1\over\xi^{2}}-{1\over{\rm i}\xi^{3}}\right)+{1\over{\rm
i}\xi}+{2\over{\rm i}\xi^{3}}\right)\right]\\!.~{}~{}~{}~{}~{}$ (35)

Analogously, one can show that $\Gamma_{\rm mir}^{(ba)}=\Gamma_{\rm
mir}^{(ab)*}$. The real part of both constants coincides with Eq. (7) in the
main text. If the $y$ and the $z$ coordinates of the position of atom $a$ and
atom $b$ are not the same, additional terms have to be taken into account in
the above derivation. However, our physical intuition tells us that the remote
interaction between atom $a$ and atom $b$ depends also in this case only on
the distance $\xi$ between the position of atom $a$ and the position of the
mirror image of atom $b$ and not on the actual distance of the atoms.

####  IV.1.3 The reset operator ${\cal L}(\rho_{\rm A})$

For completeness, we now also calculate the state of the atoms ${\cal
L}(\rho_{\rm A})$ in case of an emission. Substituting Eq. (17) into Eq.
(IV.1) and introducing the variables $\tau=\Delta t-t$ and
$\tau^{\prime}=\Delta t-t^{\prime}$ yields

$\displaystyle\mathcal{L}(\rho_{\rm A})$ $\displaystyle=$
$\displaystyle\sum_{i,j=a,b}\widetilde{\Gamma}_{\rm
mir}^{(ij)}\,\sigma_{i}^{-}\rho_{\rm A}\sigma_{j}^{+}$ (36)

with the constants $\widetilde{\Gamma}_{\rm mir}^{(ij)}$ given by

$\displaystyle\widetilde{\Gamma}_{\rm mir}^{(ij)}$ $\displaystyle=$
$\displaystyle{e^{2}\over\hbar^{2}\Delta t}\int_{0}^{\Delta t}{\rm
d}\tau\int_{0}^{\Delta t}{\rm d}\tau^{\prime}\,{\rm e}^{{\rm
i}\omega_{0}(\tau-\tau^{\prime})}\,\langle 0_{\rm
F}|\,\boldsymbol{D}^{(j)}_{12}\cdot\boldsymbol{E}(\boldsymbol{r}_{j})\,U^{\dagger}_{\rm
F}(\tau^{\prime},0)\boldsymbol{D}^{(i)}_{12}\cdot U_{\rm
F}(\tau,0)\boldsymbol{E}(\boldsymbol{r}_{i})^{\dagger}\,|0_{\rm F}\rangle\,.$
(37)

These constants have many similarities with the constants $\Gamma_{\rm
mir}^{(ij)}$ in Eqs. (20). The only differences are a missing factor 2 and a
different upper limit on the second time integral. Proceeding as in the
previous subsection, we find that evaluating Eq. (37) now leads to time
integrals of the form

$\displaystyle\int_{0}^{\Delta t}{\rm d}\tau\int_{0}^{\Delta t}{\rm
d}\tau^{\prime}\,{\rm e}^{{\rm i}\omega(\tau-\tau^{\prime})}$ $\displaystyle=$
$\displaystyle 2{\rm Re}\left(\int_{0}^{\Delta t}{\rm
d}\tau\int_{0}^{\tau}{\rm d}\tau^{\prime}\,{\rm e}^{{\rm
i}\omega(\tau-\tau^{\prime})}\right)$ (38)

with $\omega=\omega_{0}-ck$. Hence all the constants $\widetilde{\Gamma}_{\rm
mir}^{(ij)}$ are real and $\widetilde{\Gamma}_{\rm mir}^{(ij)}={\rm
Re}\left(\Gamma_{\rm mir}^{(ij)}\right)$ which yields the reset operator
$\mathcal{L}(\rho_{\rm A})$ in Eq. (II).

###  IV.2 Spontaneous decay of excited atomic states

Figure 4: (a) Time dependence of the photon emission rate $I(t)$ in Eq. (42)
for different initial populations $p$ of the excited atomic state in the
presence of mirror-mediated dipole-dipole interactions. Here ${\rm
Re}\big{(}\Gamma_{\rm mir}^{(ab)}\big{)}=0.05\,\Gamma_{\rm free}$. (b) The
rate $I(t)$ differs from the emission rate $I_{0}(t)$ of the atoms in the
absence of interactions for the same $p$. As one would expect in the case of a
broadening of spontaneous decay rates, the loss of atomic excitation happens
faster at relatively short times and slower at later times.

To determine the spontaneous decay rates of the two atoms, we now absorb all
the Hermitian terms of the conditional Hamiltonian $H_{\rm cond}$ in Eq. (II)
into their free Hamiltonian $H_{\rm A}$. This is possible since these terms
describe atomic level shift and allows us to write $H_{\rm cond}$ as

$\displaystyle H_{\rm cond}$ $\displaystyle=$ $\displaystyle H_{\rm A}-{{\rm
i}\hbar\over
2}\left[\Gamma_{+}\,L_{+}^{\dagger}L_{+}+\Gamma_{-}\,L_{-}^{\dagger}L_{-}\right]$
(39)

where the $\Gamma_{\pm}$ are the spontaneous decay rates of the two atoms in
Eq. (8) and where the $L_{\pm}=(\sigma_{a}^{-}\pm\sigma_{b}^{-})/\sqrt{2}$ are
atomic lowering operators. Hence the time evolution operator $U_{\rm
cond}(t,0)=\exp\left(-{\rm i}H_{\rm cond}t/\hbar\right)$ which describes the
dynamics of the two atoms on opposite sides of the mirror interface under the
condition of no photon emission in $(0,t)$ equals, in the interaction picture
with respect to $H_{0}=H_{\rm A}$ and $t=0$,

$\displaystyle U_{\rm cond}(t,0)$ $\displaystyle=$
$\displaystyle|11\rangle\langle 11|+|+\rangle\langle+|\,{\rm
e}^{-\Gamma_{+}t/2}+|-\rangle\langle-|\,{\rm e}^{-\Gamma_{-}t/2}$ (40)
$\displaystyle+|22\rangle\langle 22|\,{\rm e}^{-\Gamma_{\rm free}t}$

with $|\pm\rangle=(|12\rangle\pm|21\rangle/\sqrt{2}$ denoting the symmetric
and the antisymmetric collective states of the two atoms with ground states
$|1\rangle$ and excited states $|2\rangle$. Suppose an incoherent excitation
process prepares each atom with probability $p$ in its excited state, thereby
creating a statistical mixture of the atomic states $|11\rangle$,
$|12\rangle$, $|21\rangle$ and $|22\rangle$. In this case, the probability
$P_{0}(t)=\|U_{\rm cond}(t,0)\,|\psi_{\rm I}\rangle\|^{2}$ for no photon
emission in $(0,t)$ is the sum of three exponentials and equals [6, 7]

$\displaystyle P_{0}(t)$ $\displaystyle=$
$\displaystyle(1-p)^{2}+(1-p)p\left({\rm e}^{-\Gamma_{+}t}+{\rm
e}^{-\Gamma_{-}t}\right)$ (41) $\displaystyle+p^{2}\,{\rm e}^{-2\Gamma_{\rm
free}t}\,.$

In the absence of any dipole-dipole interactions, $\Gamma_{\pm}=\Gamma_{\rm
free}$ and the atoms decay exactly as they would in free space. However, in
general, this is no longer the case. For example, for $p\ll 1$, the
probability of finding both atoms in the excited state becomes negligible and
the probability density $I(t)$ for a photon emission at $t$ coincides with the
probability density $w_{1}(t)=-{\rm d}/{\rm d}t\,P_{0}(t)$ for the emission of
a first photon at $t$. Hence,

$\displaystyle I(t)$ $\displaystyle=$ $\displaystyle 2p\left[\Gamma_{\rm
free}\,{\rm cosh}\left({\rm Re}\big{(}\Gamma_{\rm
mir}^{(ab)}\big{)}t\right)\right.$ (42) $\displaystyle\left.-{\rm
Re}\big{(}\Gamma_{\rm mir}^{(ab)}\big{)}\,{\rm sinh}\left({\rm
Re}\big{(}\Gamma_{\rm mir}^{(ab)}\big{)}t\right)\right]{\rm e}^{-\Gamma_{\rm
free}t}~{}~{}$

up to first order in $p$. As illustrated in Fig. 4, this emission rate is
qualitatively different from the emission rate $I_{0}(t)$ of the atoms in the
absence of any dipole-dipole interactions. It is therefore possible to use
fluorescence lifetime measurements to detect the above described changes of
the decay rates and to obtain a signature of the remote mirror-mediated
dipole-dipole interactions which we predict in this paper.
Data availability statement: Statement of compliance with EPSRC policy
framework on research data: This publication is theoretical work that does not
require supporting research data.
Acknowledgement. We would like to thank Lewis Clark, Daniel Hodgson and Axel
Kuhn for stimulating discussions. We also acknowledge financial support from
the Oxford Quantum Technology Hub NQIT (grant number EP/M013243/1) and the
EPSRC (2115757). Statement of compliance with EPSRC policy framework on
research data: This publication is theoretical work that does not require
supporting research data.
Author Contributions: Conceptualisation, G.J. and A.B.; methodology, N.F.-W.,
B.D. and A.B.; formal analysis, N.F.-W., B.D. and A.B.; writing and original
draft preparation, N.F.-W., B.D. and A.B.; writing and review and editing, T.
M. and G.J.; visualization, B.D. and A.B.; supervision, A.B.; funding
acquisition, N.F.-W., G.J. and A.B. All authors have read and agreed to the
published version of the manuscript.
Conflicts of Interest: The authors declare no conflict of interest.

## References

  * [1] Scully, M. O. & Drühl, K. Quantum eraser: A proposed photon correlation experiment concerning observation and ”delayed choice” in quantum mechanics. Phys. Rev. A 25, 2208 (1982).
  * [2] Eichmann. U. et al. Young’s interference experiment with light scattered from two atoms. Phys. Rev. Lett. 70, 2359 (1993).
  * [3] Schön, C. & Beige, A. Analysis of a two-atom double-slit experiment based on environment-induced measurements. Phys. Rev. A 64, 023806 (2001).
  * [4] Beige, A., Schön, C. & Pachos, J. Interference of spontaneously emitted photons. Fortschr. Phys. 50, 594 (2002).
  * [5] Beige, A. A new take on the slit experiment. Nat. Photon. 10, 290 (2016).
  * [6] Hegerfeldt, G. C. How to reset an atom after a photon detection: Applications to photon counting processes. Phys. Rev. A 47, 449 (1993).
  * [7] Stokes, A. et al. Extending the validity range of quantum optical master equations. Phys. Rev. A 85, 053805 (2012).
  * [8] Beige, A. & Hegerfeldt, G. C. Transition from antibunching to bunching for two dipole-interacting atoms. Phys. Rev. A 58, 4133 (1998).
  * [9] Masson, S. J. et al. A. Many-Body Signatures of Collective Decay in Atomic Chains. Phys. Rev. Lett. 125, 263601 (2020).
  * [10] Wolf, S. et al. Light of Two Atoms in Free Space: Bunching or Antibunching? Phys. Rev. Lett. 124, 063603 (2020).
  * [11] Lim, Y. L., Beige, A. & Kwek, L. C. Repeat-Until-Success Linear Optics Distributed Quantum Computing. Phys. Rev. Lett. 95, 030505 (2005).
  * [12] Barrett S. D. & Kok, P. Efficient high-fidelity quantum computation using matter qubits and linear optics. Phys. Rev. A 71, 060310(R) (2005).
  * [13] Duan, L.-M. & Monroe, C. Colloquium: Quantum networks with trapped ions. Rev. Mod. Phys. 82, 1209 (2010).
  * [14] Hensen, B. et al. Loophole-free Bell inequality violation using electron spins separated by 1.3 kilometres. Nature 526, 682 (2015).
  * [15] Stephenson, L. et al. High-rate high-fidelity entanglement of qubits across an elementary quantum network. Phys. Rev. Lett. 124, 110501 (2020).
  * [16] Araneda, G. et al. Interference of single photons emitted by entangled atoms in free space. Phys. Rev. Lett. 120, 193603 (2018).
  * [17] Dicke, R. H. Coherence in Spontaneous Radiation Processes. Phys. Rev. 93, 94 (1954).
  * [18] Agarwal, G. S. Quantum statistical theories of spontaneous emission and their relation to other approaches. Springer Tracts in Modern Physics; Quantum Optics, Vol. 70 (Springer Verlag Berlin, 1974).
  * [19] Gross, M. & Haroche, S. Superradiance: An essay on the theory of collective spontaneous emission. Phys. Rep. 93, 301 (1982).
  * [20] Tanas R. & Ficek, Z. Entangled states and collective nonclassical effects in two-atom systems, Phys. Rep. 372, 369 (2002).
  * [21] Furtak-Wells, N. et al. Quantizing the electromagnetic field near two-sided partially-transparent mirrors. Phys. Rev. A 97, 043827 (2018).
  * [22] Klimov, V. V. & Letokhov, V. S. Resonant energy exchange at nanoscale curved interface. Chem. Phys. Lett. 285, 313 (1998).
  * [23] Klimov, V. V. & Letokhov, V. S. Resonant interaction between two atomic dipoles separated by the surface of a dielectric microsphere. Phys. Rev. A 58, 3235 (1998).
  * [24] Drexhage, K. H. Influence of a dielectric interface on fluorescence decay time. J. Lum. 1, 693 (1970).
  * [25] Worthing, P. T., Amos, R. M. & Barnes, W. L. Modification of the spontaneous emission rate of Eu3+ ions embedded within a dielectric layer above a silver mirror. Phys. Rev. A 59, 865 (1999).
  * [26] Eschner, J. et al. Light interference from single atoms and their mirror images. Nature 413, 495 (2001).
  * [27] Hoi, I. C. et al. Probing the quantum vacuum with an artificial atom in front of a mirror. Nat. Phys. 11, 1045 (2015).
  * [28] Palacino R. et al. Tuning the collective decay of two entangled emitters by means of a nearby surface. J. Phys. B 50, 1969 (2017).
  * [29] Sinha, K., Venkatesh, B. P. & Meystre, P. Collective effects in Casimir-Polder forces. Phys. Rev. Lett. 121, 183605 (2018).
  * [30] Zhou, W., Rizzuto, L. & Passante, R. Vacuum fluctuations and radiation reaction contributions to the resonance dipole-dipole interaction between two atoms near a reflecting boundary. Phys. Rev. A 97, 042503 (2018).
  * [31] Pendry, J. B. Negative Refraction Makes a Perfect Lens. Phys. Rev. Lett. 85, 3966 (2000).
  * [32] Andrew, P. & Barnes, W. L. Förster Energy Transfer in an Optical Microcavity. Science 290, 785 (2000).
  * [33] Andrew, P. & Barnes, W. L. Energy Transfer Across a Metal Film Mediated by Surface Plasmon Polaritons. Science 306, 1002 (2004).
  * [34] Kästel, J. & Fleischauer, M. Suppression of spontaneous emission and superradiance over macroscopic distances in media with negative refraction. Phys. Rev. A 71, 011804(R) (2005).
  * [35] Biehs, S.-A., Menon, V. M. & Agarwal, G. S. Long-range dipole-dipole interaction and anomalous Förster energy transfer across a hyperbolic metamaterial. Phys. Rev. B, 93, 245439 (2016).
  * [36] Biehs, S.-A. & Agarwal, G. S. Qubit entanglement across $\epsilon$-near-zero media. Phys. Rev. A 96, 022308 (2017).
  * [37] Deshmukh, R. et al. Long-Range Resonant Energy Transfer Using Optical Topological Transitions in metamaterials. ACS Photonics 5, 2737 (2018).
  * [38] Newman, R. D. Observation of long-range dipole-dipole interactions in hyperbolic metamaterials. Sci. Adv. 4, eaar5278 (2018).
  * [39] Karanikolos, V. Entanglement of quantum. emitters interacting through an ultra-thin noble metal nanodisk. Opt. Express 28, 24173 (2020).
  * [40] Gonzalez-Tudela, A. et al. Entanglement of two qubits mediated by one-dimensional plasmonic waveguides. Phys. Rev. Lett. 106, 020501 (2011).
  * [41] Le Kien, F. & Rauschenbeutel, A. Propagation of nanofiber-guided light through an array of atoms. Phys. Rev. A 90 063816 (2014).
  * [42] Hung, C.-L. et al. Quantum spin dynamics with pairwise-tunable, long-range interactions. Proc. Natl. Acad. Sci. U.S.A. 113, E4946 (2016).
  * [43] Solano, P. et al. Super- and sub-radiance reveal infinite-range interactions through a nanofiber. Nat. Comm. 8, 1857 (2017).
  * [44] Southall, J. et al. Locally-acting mirror Hamiltonians. J. Mod. Opt. 68, 647 (2021).
  * [45] Hodgson, D. et al. Local photons. Front. Photon. 3, 978855 (2022).
  * [46] Al Ghamdi, A. et al. A. Remote Non-Invasive Fabry-Perot Cavity Spectroscopy for Label-Free Sensing. Sensors 23, 385 (2023).
  * [47] Bennett, R., Barlow T. M. & Beige, A. A physically-motivated quantisation of the electromagnetic field. Eur. J. Phys. 37, 014001 (2016).