# Analytic solution to the nonlinear generation of squeezed states in a
thermal bath

Paul R. B. Hughes<EMAIL_ADDRESS>Marc M. Dignam  Department of Physics,
Engineering Physics and Astronomy, Queen’s University, Kingston, ON K7L 3N6,
Canada

###### Abstract

We model squeezed state generation in a lossy optical cavity in the presence
of a thermal bath using the Lindblad master equation. We show that the exact
solution is a squeezed thermal state, where thermal photons arise both from
loss and from the thermal bath. We derive an exact, closed-form solution for
the evolution of the quadrature uncertainty arising from pulsed degenerate
spontaneous parametric down conversion in the cavity. We apply this solution
under different pump conditions and show in detail how the thermal environment
reduces quadrature squeezing as well as the second order coherence function.

Introduction. Nonlinear optical processes such as spontaneous parametric down
conversion (SPDC) and spontaneous four wave mixing are often used to generate
nonclassical states of light, such as photon pairs, single-mode quadrature
squeezed states, multimode squeezed states, and entangled optical modes [1, 2,
3, 4, 5]. Photon pairs can be used as a heralded single-photon source or as
entangled two-photon states [6, 7]. Quadrature squeezed states can be used to
reduce the uncertainty in interferometric measurements [8, 9] or to create
continuous variable entanglement [10, 11]. Because the nonlinear interactions
are generally quite weak, one usually requires a resonator to form an optical
parametric oscillator (OPO) to enhance the process. Some common resonators are
microring resonators [12, 13, 14] and Fabry-Perot cavities [4, 15]. The use of
a resonator has been shown to lower the uncertainty in one quadrature below
vacuum fluctuations at the expense of the other, but with a steady-state
squeezing limit of 3dB in the resonator [1, 16]. Larger squeezing can been
obtained for light coupled out of the resonator [15, 17, 18, 19, 20, 21, 22],
but within it, alternative methods are required to overcome the 3dB limit.
These include pulsed excitation [12], quantum feedback [23, 24, 25], and
dissipation [26, 27].

The steady-state quantum fluctuations of the state in an OPO can be derived
using the Langevin equations for a stochastic process [28, 2, 29, 30]. When
loss and detuning are not considered in these systems, the signal field in the
OPO is a squeezed vacuum state. Other groups have considered the effects of
detuning from resonance on squeezing in the OPO in the steady state [31, 32].
The effects of loss and a thermal bath on the generation and evolution of the
density operator of the light in an OPO can be modeled using the Lindblad
master equation (LME). Recently it was shown that when there is no thermal
bath, the exact solution to the LME is a squeezed thermal state (STS) [16].

At optical or near infrared frequencies, the thermal effects of an environment
on the generation and nature of the squeezed states is negligible at or below
room temperature when employing SPDC in a resonator. At lower frequencies of a
few tens of terahertz or less, thermal noise can have a significant effect on
the generation, evolution, and final state. In particular quadrature squeezing
can be greatly reduced unless one cools the system to millikelvin temperatures
[27]. It is therefore important to be able to accurately and efficiently model
the evolution dynamics in such systems and to quantify the effects of
temperature on the final squeezed state for CW and pulsed pumping
configurations. To this end, in this work, we derive the exact solution to the
LME for SPDC in a lossy OPO coupled to a thermal bath and show that the
density operator is that of a STS. With this exact solution, we are able to
derive a closed-form solution for the evolution of the quadrature uncertainty
for an arbitrary, un-chirped pump pulse and to examine the evolution of the
squeezing parameter, the thermal photon number and the second-order quantum
coherence function.

The paper is organized as follows. We first outline the theory behind the
generation of the signal field in a resonator. We show that the solution to
the LME is a STS, where the thermal photon number and squeezing parameter
evolution are described by three coupled first-order differential equations.
Using these equations, we derive closed form solutions for the quadrature
uncertainties, which to the best of the authors’ knowledge have never been
derived previously. Next, we examine the transient and steady-state properties
of the system excited by a constant-amplitude pump pulse, presenting an exact
analytic solution for the uncertainties. Using the second order quantum
coherence and the quadrature uncertainty, we investigate the nonclassicallity
of the light and discuss the threshold where the quadrature is squeezed below
vacuum noise. Finally, we present the results for a Gaussian pump pulse and
examine the relationship between the pulse amplitude and the quadrature
squeezing as a function of the bath temperature.

Theory. We consider the generation of a squeezed state in a single mode of a
resonant cavity with frequency $\omega$. Shown schematically in Fig. 1, the
system consists of a resonant cavity that is coupled to a thermal bath of
photons and is excited by a coherent optical pump. The pump operates at a
frequency $\omega_{p}=2\omega$ and generates signal photons in the resonator
through SPDC. The pump is a coherent state with time-dependent coherent state
amplitude $\alpha(t)=\alpha_{0}(t)e^{-i\omega_{p}t}$, where $\alpha_{0}(t)$ is
the pump envelope. Thus, we treat the pump classically in the undepleted pump
approximation. When the interaction with the bath is neglected, the system
Hamiltonian is given by [16]

$H=\hbar\omega{b}^{\dagger}b+\alpha(t)\gamma{b^{\dagger}}^{2}+\alpha^{*}(t)\gamma^{*}b^{2},$
(1)

where $b^{\dagger}$ ($b$) is the creation (annihilation) operator of photons
in the cavity and $\gamma=\hbar\omega_{p}\chi_{eff}^{(2)}/n_{eff}^{2}$ is the
coupling coefficient of the pump field to the signal field for an effective
second order nonlinear susceptibility $\chi_{eff}^{(2)}$ and refractive index
$n_{eff}$ in the cavity.

The signal mode in the cavity is coupled to a thermal bath at temperature
$T_{b}$, which has a mean photon number,
$n_{b}=(\exp(\hbar\omega/kT_{b})-1)^{-1}$ at the signal frequency. The density
operator of the cavity $\rho(t)$ evolves according to the LME [33]

$\displaystyle\begin{aligned}
\frac{d}{dt}\rho(t)=-\frac{i}{\hbar}\commutator{H}{\rho(t)}&+\Gamma(n_{b}+1)D[{b}](\rho(t))\\\
&+\Gamma n_{b}D[{b}^{\dagger}](\rho(t)),\end{aligned}$ (2)

where $\Gamma$ is the power decay constant of the cavity photons into the
bath, while

$D[F](\rho)\equiv F\rho
F^{\dagger}-\frac{1}{2}\anticommutator{F^{\dagger}F}{\rho}$ (3)

is the dissipator, which accounts for the two-way coupling between bath and
cavity.

$\hbar\omega_{p}$$\alpha(t)$$\gamma$$T_{b}$$\hbar\omega$$\hbar\omega$$\Gamma$

Figure 1: Schematic diagram of the system. A coherent state $\alpha(t)$ at
frequency $\omega_{p}$ pumps the resonator which is coupled $(\Gamma)$ to the
environment at temperature $T_{b}$, generating pairs of signal photons at
$\omega$.

In a previous work, it was shown that for the special case where the bath is
at zero temperature ($n_{b}=0$), the exact solution to Eq. 2 is a STS [16]. In
this work, we prove that at non-zero temperatures, the solution is still a
STS, but that the evolution of the squeezing and thermal temperature depends,
in general, on the bath temperature.

We find that as long as the initial state is the vacuum, a thermal state, or a
STS, the exact solution to the above LME is the time-dependent STS:

$\displaystyle\rho(t)=S(\xi(t))\rho_{T}(n_{th}(t))S^{\dagger}(\xi(t)),$ (4)

where

$S(\xi)=\exp\left[\frac{1}{2}(\xi^{*}b^{2}-\xi{b^{\dagger}}^{2})\right]$ (5)

is the squeezing operator, with the time-dependent, complex squeezing factor
$\xi(t)=u(t)e^{i\phi(t)}$, and

$\displaystyle\rho_{T}(n_{th})=\frac{1}{1+n_{th}}\left(\frac{n_{th}}{1+n_{th}}\right)^{b^{\dagger}b}$
(6)

is a thermal state, with a time-dependent thermal population,
$n_{th}\left(t\right)$.

To show that this is the exact solution, we write the density operator in the
form

$\displaystyle\rho(t)=S(\xi)\rho_{T}^{1/2}(n_{th})O(t)\rho_{T}^{1/2}(n_{th})S^{\dagger}(\xi).$
(7)

We then need to prove that the operator,

$\displaystyle
O(t)=\rho_{T}^{-1/2}(n_{th})S^{\dagger}(\xi)\rho(t)S(\xi)\rho_{T}^{-1/2}(n_{th})$
(8)

is simply the identity operator for all time. In the supplementary material,
we show that this is indeed the case as long as the thermal photon number, the
squeezing amplitude, and phase evolve according to the following three coupled
first-order differential equations:

$\displaystyle\frac{dn_{th}}{dt}=\Gamma\left[n_{b}\cosh(2u)+\sinh^{2}(u)-n_{th}\right],$
(9)
$\displaystyle\frac{du}{dt}=-\frac{i}{\hbar}(\gamma^{*}\alpha^{*}e^{i\phi}-\gamma\alpha
e^{-i\phi})-\frac{\Gamma}{2}\sinh(2u)\frac{2n_{b}+1}{2n_{th}+1},$ (10)
$\displaystyle\frac{d\phi}{dt}=-2\omega+\frac{2}{\hbar}\frac{\cosh(2u)}{\sinh(2u)}(\gamma^{*}\alpha^{*}e^{i\phi}+\gamma\alpha
e^{-i\phi}).$ (11)

Eqs. 9, 10 and 11, are the dynamic equations valid for any initial STS and for
any $\alpha(t)$. In all that follows, we restrict ourselves to unchirped pump
pulses, such that $\alpha_{0}(t)\gamma=|\alpha_{0}(t)\gamma|e^{i\theta}$,
where $\theta$ is a time-independent phase and we assume that the initial
state is an unsqueezed thermal state, such that $u(0)=0$. To avoid a
divergence on the left hand side of Eq. 11 that arises at $t=0$ for an
unsqueezed state, we impose the initial condition on $\phi$ that
$(\gamma^{*}\alpha^{*}(0)e^{i\phi(0)}+\gamma\alpha(0)e^{-i\phi(0)})=0$, or
$\phi(0)=\theta+\pi/2$. Using this in Eq. 11, we see that for all time
$(\gamma^{*}\alpha^{*}(t)e^{i\phi(t)}+\gamma\alpha(t)e^{-i\phi(t)})=0$ and
$\phi(t)=\theta+\pi/2-2\omega t$.

We now define the pump function,

$g(t)\equiv 4|\alpha_{0}(t)\gamma|/\hbar\Gamma,$ (12)

which is the ratio of the pumping rate to the loss rate, such that $g(t)=1$ is
the critical pump rate at which the injection of photons is exactly balanced
by the loss. Using this definition and the above initial condition, the
equations of motion become

$\displaystyle\frac{dn_{th}}{dt}=\Gamma\left[n_{b}\cosh(2u)+\sinh^{2}(u)-n_{th}\right],$
(13) $\displaystyle\frac{du(t)}{dt}=\frac{\Gamma
g(t)}{2}-\frac{\Gamma\sinh(2u)}{2}\frac{2n_{b}+1}{2n_{th}+1},$ (14)

with

$\phi(t)=\theta+\pi/2-2\omega t.$ (15)

The above dynamic equations contain an explicit dependence on the temperature
of the bath; however, in the limit that $n_{b}=0$, they reduce to what we
obtained in our previous $T_{b}=0$ work [16]. Additionally, in the simple case
that there is no pump present, but the initial state is a thermal state that
is not at the bath temperature, these equations show that the system remains a
thermal state with the thermal photon number evolving as
$n_{th}(t)=n_{b}+(n_{th}(0)-n_{b})e^{-\Gamma t}$, a result that has been shown
previously using the LME [34].

We can see from Eq. 14 that the bath population increases the decay rate of
the squeezing factor $u(t)$. However, this contribution is not present for the
typical initial state in which the system is in equilibrium with the
environment. To see this, let $n_{th}^{0}$ be the thermal population that
would arise if $T_{b}=0$. We define it using the equation

$\displaystyle 2n_{th}+1=(2n_{b}+1)(2n_{th}^{0}+1).$ (16)

Using this in Eq. 13, we find that

$\displaystyle\frac{dn_{th}^{0}}{dt}=\Gamma[\sinh^{2}(u)-n_{th}^{0}],$ (17)

which is exactly the evolution of the thermal photon number when the bath
temperature is zero. We can also rewrite Eq. 14 using $n_{th}^{0}$ as

$\displaystyle\frac{du(t)}{dt}=\frac{\Gamma
g(t)}{2}-\frac{\Gamma\sinh(2u)}{2(2n_{th}^{0}+1)}.$ (18)

The bath population is still implicitly present in these equations, since from
Eq. 16, the initial value of $n_{th}^{0}$ is given by

$\displaystyle n_{th}^{0}(0)=\frac{n_{th}(0)-n_{b}}{2n_{b}+1}.$ (19)

However, in the case where the cavity begins in a thermal state in equilibrium
with the bath ($n_{th}(0)=n_{b}$), the evolution of $n_{th}^{0}$ and thus
squeezing factor $u$ are independent of the bath temperature, while the actual
thermal population only depends on $n_{b}$ through the prefactor of $2n_{b}+1$
(see Eq. 16) 111We can see from Eq. 19 that these properties arise not just
for an initial state in equilibrium with the thermal bath, but any initial
thermal population that satisfies $n_{th}(0)=n_{b}+a(2n_{b}+1)$. The dynamics
of the squeezing amplitude will then be independent of the bath, but the decay
will scale by the arbitrary factor $a$.. We will examine the dependence of the
thermal population, the total population, and the squeezing amplitude on the
bath temperature and initial thermal population in more detail later in this
letter.

We now define the two quadrature operators,

$\displaystyle X=b^{\dagger}e^{-i\beta(t)}+be^{i\beta(t)},$ (20)
$\displaystyle Y=-i(b^{\dagger}e^{-i\beta(t)}-be^{i\beta(t)}),$ (21)

where $\beta(t)$ is the local oscillator phase. For $\beta(t)\equiv\omega t$,
the system is squeezed in $X$ and antisqueezed in $Y$. For a STS, the
uncertainties in these quadratures can be shown to be given by $\Delta
X^{2}=(2n_{th}+1)e^{-2u}$ and $\Delta Y^{2}=(2n_{th}+1)e^{2u}$ [36], which
allows us to determine the evolution of the squeezing from the evolution of
$n_{th}$ and $u$. Alternatively, we can derive differential equations for the
quadrature uncertainties. Taking the time derivative of $\Delta X^{2}$ and
using Eqs. 13 and 14 and the hyperbolic identities, we find that

$\displaystyle\frac{d}{dt}\Delta
X^{2}=\Big{[}2\frac{dn_{th}}{dt}-2(2n_{th}+1)\frac{du}{dt}\Big{]}e^{-2u}$
$\displaystyle\begin{aligned}
=\Gamma\Big{[}&(2n_{b}+1)\left(\cosh(2u)+\sinh(2u)\right)e^{-2u}\\\
&-\left(1+g(t)\right)(2n_{th}+1)e^{-2u}\Big{]},\end{aligned}$ (22)

which can also be written as

$\displaystyle\frac{d}{dt}\Delta X^{2}=\Gamma\left[(2n_{b}+1)-(1+g(t))\Delta
X^{2}\right].$ (23)

Similarly,

$\displaystyle\frac{d}{dt}\Delta Y^{2}=\Gamma\left[(2n_{b}+1)-(1-g(t))\Delta
Y^{2}\right].$ (24)

These equations show that the squeezing dynamics depend only on the pumping
strength and the thermal bath population. Furthermore, the evaluation of the
quadrature squeezing only requires the solution of a single first-order
differential equation, which is directly solvable using standard techniques.
To the authors’ knowledge, this is the first time a closed-form solution has
been derived for the nonlinear generation and evolution of the quadrature
uncertainty in a lossy cavity.

Constant Pump and Steady State. In this section, we examine the early-time
evolution and steady-state solution of the system when it is excited by a pump
that has a constant strength, $g_{0}$ that is turned on at $t=0$.
Before analysing the quadrature uncertainties directly, we return to the
question of the dependence of the STS parameters on the bath temperature and
initial state of the system. Recall that for the special case where
$n_{th}(0)=n_{b}$, the squeezing factor $u(t)$ is independent of $n_{b}$. As
we now show, when the initial state is a thermal state at a different
temperature from the environment, $u(t)$ is still only weakly dependent on
both the initial thermal population and the bath temperature.

Figure 2: (a) Total population $n(t)$ and (b) squeezing factor $u(t)$ as a
function of time for an initial thermal state with $n_{th}(0)=1.5$. The cavity
is pumped with a continuous wave pulse of strength $g(t)=g_{0}=0.8$. The
different curves are for different thermal bath populations, $n_{b}$.

In Fig. 2, we plot the total photon number and squeezing amplitude as a
function of time for a cavity pumped by a continuous pulse excitation. It is
prepared in the same initial state each time but coupled to environments at
different bath temperatures. From Fig. 2(a), we see that when $n_{b}$ is
increased, the photon number increases by much more than simply $n_{b}$.
Meanwhile, Fig. 2(b) shows that $u(t)$ exhibits a small dependence on $n_{b}$
at early times and that even this dependence disappears as
$t\rightarrow\infty$. We can determine the steady-state characteristics of the
continuous wave pump by setting the derivatives in Eqs. 14 and 13 to zero. For
the steady-state squeezing amplitude, we obtain

$\displaystyle u^{ss}=\frac{1}{2}\tanh^{-1}(g_{0}),$ (25)

which is independent of the environment and only exists for pumping below the
critical pump strength, $g_{0}=1$. The steady-state thermal and total
populations are given by

$\displaystyle\begin{aligned} n_{th}^{ss}&=n_{b}+\sinh^{2}(u)(2n_{b}+1)\\\
&=\frac{1}{2}\left(\frac{2n_{b}+1}{\sqrt{1-g_{0}^{2}}}-1\right),\end{aligned}$
(26)

and

$\displaystyle n^{ss}=\frac{2n_{b}+g_{0}^{2}}{2(1-g_{0}^{2})},$ (27)

where, to determine the total population, we have used the relation for a STS
that $n=n_{th}\cosh(2u)+\sinh^{2}(u)$ [36]. Thus, as discussed earlier, the
thermal environment adds many more photons to the system than just $n_{b}$,
but it does nothing to the steady state squeezing factor.

We now consider the evolution of quadrature uncertainties. The dynamic
equations, Eqs. 23 and 24, can be solved exactly for an arbitrary time-
dependent pump $g(t)$, but we first consider solutions for a constant-pump
excitation, where $g(t)=g_{0}\neq 1$ for $t>0$. The exact solutions are

$\displaystyle\Delta X^{2}(t)=\frac{2n_{b}+1}{1+g_{0}}+\left(\Delta
X^{2}(0)-\frac{2n_{b}+1}{1+g_{0}}\right)e^{-\Gamma(1+g_{0})t},$ (28)
$\displaystyle\Delta Y^{2}(t)=\frac{2n_{b}+1}{1-g_{0}}+\left(\Delta
Y^{2}(0)-\frac{2n_{b}+1}{1-g_{0}}\right)e^{-\Gamma(1-g_{0})t}.$ (29)

When the system is prepared as a thermal state in equilibrium with the bath,
both quadratures start with a value of $2n_{b}+1$, and Eqs. 28 and 29 simplify
to

$\displaystyle\Delta
X^{2}(t)=\frac{2n_{b}+1}{1+g_{0}}\left[1+g_{0}e^{-\Gamma(1+g_{0})t}\right],$
(30) $\displaystyle\Delta
Y^{2}(t)=\frac{2n_{b}+1}{1-g_{0}}\left[1-g_{0}e^{-\Gamma(1-g_{0})t}\right].$
(31)

Increased squeezing in $X$ will usually be accompanied by increased anti-
squeezing in $Y$. Because of this, it is important to examine the merits of
using a short, strong pump pulse rather than a long, weak pulse to squeeze the
signal. To this end, we consider the antisqueezing at the time $\tau_{1}$ at
which the squeezing in $X$ reaches the threshold value of $\Delta X^{2}=1$.
For $n_{th}(0)=n_{b}>0$, we have from Eqs. 30 and 31 that

$\displaystyle\Delta
Y^{2}(\tau_{1})=\frac{2n_{b}+1}{1-g_{0}}\left[1-\left(\frac{1}{g_{0}}\frac{g_{0}-2n_{b}}{2n_{b}+1}\right)^{\frac{1-g_{0}}{1+g_{0}}}\right].$
(32)

In the limit of weak pumping ($g_{0}\ll 1$), $\Delta
Y^{2}(\tau_{1})\rightarrow 2n_{b}/g_{0}$, while in the strong pumping limit
($g_{0}>>1$), $\Delta Y^{2}(\tau_{1})\rightarrow
2n_{b}(2n_{b}+1)/(g_{0}-2n_{b})$. In both limits, increasing $g_{0}$ reduces
the uncertainty in $\Delta Y^{2}$ at $t=\tau_{1}$. Thus, in order to avoid
excess growth in the anti-squeezed uncertainty, it is always beneficial to
maximize the pumping strength if a specified squeezing is desired. This is
because the longer it takes to reach a desired squeezing level, the more
thermal photons will be generated due to loss.

We can see from Eq. 23 that the squeezed quadrature will reach a steady-state
value of

$\displaystyle\Delta X_{min}^{2}=\frac{2n_{b}+1}{1+g_{0}},$ (33)

for all values of $g_{0}$. The anti-squeezing, however, will only reach steady
state for $g_{0}<1$, diverging otherwise. Below critical pumping, this anti-
squeezing maximum is

$\displaystyle\Delta Y^{2}_{max}=\frac{2n_{b}+1}{1-g_{0}}.$ (34)

In the steady-state, the squeezing is limited to a minimum quadrature
uncertainty of $(2n_{b}+1)/2$, and therefore we cannot achieve any squeezing
if $n_{b}\geq 0.5$, which agrees with the results for thermalized squeezed
states by previous authors [37].

Figure 3: Second order simultaneous coherence function $g^{(2)}(t)$ as a
function of time for $n_{th}(0)=n_{b}>0$ for different bath populations,
$n_{b}$. The pump starts at $t=0$ and has a constant amplitude (a) $g_{0}=0.2$
and (b) $g_{0}=0.8$.

We now consider the evolution of the equal-time, second order quantum
coherence function [2]

$g^{(2)}\left(t\right)\equiv\frac{Tr\\{b^{\dagger}b^{\dagger}bb\rho\left(t\right)\\}}{n^{2}\left(t\right)},$
(35)

which quantifies the correlation between two simultaneous photon measurements
at time $t$, where $g^{(2)}(t)>1$ indicates super-Poissonian statistics and
photon bunching. In a thermal state $g^{(2)}=2$, while in a STS [36]

$\displaystyle
g^{(2)}=2+\frac{(2n_{th}+1)^{2}\sinh^{2}(2u)}{[(2n_{th}+1)\cosh(2u)-1]^{2}}.$
(36)

In the large population limit ($n\rightarrow\infty$), the STS value approaches
the squeezed vacuum state coherence $g^{(2)}=3$ [16]. In Fig. 3, we plot
$g^{(2)}(t)$ as a function of time for two different constant pumping
strengths and several different bath temperatures, all with the initial
condition $n_{th}(0)=n_{b}$. We see that $g^{(2)}$ peaks at early times before
settling down to a lower steady-state value. As the environmental population
is increased, the peak and steady-state values are decreased and the peak
occurs at a later time. When the pump is increased, the peak is larger and
occurs at an earlier time, while the steady state value is reduced. Using Eqs.
36, 25 and 26, we determine the steady-state coherence to be

$\displaystyle
g^{(2)}_{ss}=2+\left(\frac{(2n_{b}+1)g_{0}}{2n_{b}+g_{0}^{2}}\right)^{2}.$
(37)

The coherence peak above the steady-state value arises as the state
transitions from a thermal state to a STS with higher steady-state coherence.
As the pumping begins, if $n_{b}\ll 1$, the population is still small but many
squeezed pairs are being created before they can be removed by loss and before
the total population becomes too large; because $g^{(2)}$ is normalized to the
square of the total population, this results in a larger $g^{(2)}$ during this
time period. When the bath temperature is larger though, $g^{(2)}$ is
suppressed by the existing thermal population, so we find that the peak is
significantly reduced or even absent.

Figure 4: Time, $\tau_{p}$ for the coherence function $g^{2}(t)$ to reach a
maximum value as a function of $n_{b}$ and $g_{0}$ for a constant pump and
initial condition $n_{th}(0)=n_{b}$. The red dashed line indicates
$g_{0}=2n_{b}$. Note that the peak time is on a logarithmic scale.

We can determine when the transition to an STS fails to create a peak in
$g^{(2)}(t)$ by finding the time when the coherence is maximized. In Fig. 4,
we plot this peak time, $\tau_{p}$ as a function of the pump strength and bath
population. We note a clear distinction between the coherence peak times when
squeezing is possible ($g_{0}>2n_{b}$, upper left) and when it is not (lower
right). When squeezing is not possible, the coherence reaches the steady-state
value monotonically and does not peak.

We have seen from Eq. 37 that the steady state coherence will decrease with
larger pumping, while the peak value will increase. In Fig. 5, we compare the
maximum and steady-state coherence as a function of pump strength and bath
population. Again, we see that for $g_{0}<2n_{b}$, the peak value of $g^{(2)}$
is nearly identical to the steady-state value 222Direct comparison between the
theoretical steady-state value in Eq. 37 and the maximum of the numeric
simulation shows a difference of less than $10^{-8}$ for the region
$g_{0}<2n_{b}$, which can be attributed to the limits of the computational
precision..

Figure 5: Maximum (a) and steady-state (b) second order coherence function,
$g^{(2)}$ as a function of $n_{b}$ and $g_{0}$ for a constant pump and initial
condition $n_{th}(0)=n_{b}$. Note that the two plots are almost identical
below the red dashed line where $g_{0}<2n_{b}$.

Arbitrary Pump Pulse. We now examine the solutions for an arbitrary pump
envelope and for a Gaussian pulse. The linear, non-homogeneous ODE in Eq. 23
has closed-form solution [39],

$\Delta
X^{2}(t)=\left[\Gamma(2n_{b}+1)\int_{0}^{t}q(\tilde{t})d\tilde{t}+q(0)\Delta
X^{2}(0)\right]q^{-1}(t),$ (38)

where

$q(t)\equiv\exp\left(\Gamma\int(1+g(t))dt\right).$ (39)

If the initial state is a thermal state in equilibrium with the bath, this
simplifies to

$\displaystyle\Delta
X^{2}(t)=(2n_{b}+1)\left[\Gamma\int_{0}^{t}q(t^{\prime})dt^{\prime}+q(0)\right]q^{-1}(t),$
(40)

which includes the bath temperature only as a prefactor. This means that for
the usual initial condition of $n_{th}(0)=n_{b}$, one only needs to solve a
single equation to obtain the time evolution of quadrature variance for all
bath temperatures. We note that the expression for $\Delta Y^{2}(t)$ is
identical to the one for $\Delta X^{2}(t)$, but with $g(t)\rightarrow-g(t)$.
We now consider the particular example of excitation by a Gaussian pulse, with
an envelope given by

$\displaystyle
g(t)=g_{0}\exp\left[-\frac{1}{2}\frac{\Gamma^{2}(t-t_{o})^{2}}{\sigma^{2}}\right].$
(41)

Figure 6: Squeezed quadrature $\Delta X^{2}$ as a function of time for a
Gaussian pulse envelope with $g_{0}$ = 5, $\sigma=1/\sqrt{2}$, and $\Gamma
t_{o}=2.5$ for the initial condition $n_{th}(0)=n_{b}$.

In Fig. 6, we plot the pump envelope as a function of time, along with the
squeezed quadrature uncertainty for three different thermal bath populations.
We see that the different uncertainty profiles are identical apart from the
scaling factor $2n_{b}+1$ in Eq. 40. In particular, the uncertainty is
minimized at the same time, $\tau_{0}$ for all temperatures. Unless the pulse
is very short relative to the cavity lifetime, this minimum will occur very
close to the time, $\tau_{M}$ at which the pump reaches its maximum value.
Therefore, we can approximate the pump strength at the minimum uncertainty by
the maximum pump strength, as was done in Ref. [12]. Doing this, we obtain

$\displaystyle\Delta X^{2}_{min}\approx\frac{2n_{b}+1}{1+g(\tau_{M})}.$ (42)

Note that for the given pump pulse, the quadrature squeezing below shot noise
disappears for $n_{b}\gtrsim 2.5$. However, with a well-chosen pulse strength
and shaping, we see that substantially more squeezing can be achieved than is
possible in steady-state.

Conclusion. In this work, we derived a closed-form solution to the Lindblad
master equation, and used it to determine the effects of a thermal environment
on the generation of quadrature squeezing via SPDC in a lossy resonator. We
proved that the solution is a squeezed thermal state, with contributions to
the thermal population coming from loss to and photons from the thermal bath.
We derived a closed-form solution for the evolution of the quadrature
uncertainty for an arbitrary un-chirped classical pump pulse and applied it
for both a constant and a Gaussian pump pulse. We found that the thermal bath
reduces quadrature squeezing and the equal-time, second order coherence
function.

The results presented in this work can be used to help determine the
temperature, loss, and pump requirements for squeezed state generation in the
few-terahertz regime, where room temperature environmental photons can
significantly degrade squeezing. This will be particularly relevant for
squeezed state generation in microwave cavities or optomechanical systems [40,
41], where pulse optimization will be necessary to overcome thermal
limitations. In future work, we plan to extend these results to two-mode
cavity system to determine the effect of a thermal environment on the
entanglement correlation variance and to directly examine the generation of
squeezing in optomechanical systems.

Acknowledgements This work was supported by the Canada Foundation for
Innovation and the Natural Sciences and Engineering Research Council of Canada
(NSERC).

## References

  * Milburn and Walls [1981] G. Milburn and D. Walls, Production of squeezed states in a degenerate parametric amplifier, Optics Communications 39, 401 (1981).
  * Walls and Milburn [2007] D. F. Walls and G. J. Milburn,  _Quantum Optics_ (Springer Berlin, Heidelberg, 2007).
  * Walls [1983] D. F. Walls, Squeezed states of light, Nature (London) 306, 141 (1983).
  * Wu  _et al._ [1986] L.-A. Wu, H. J. Kimble, J. L. Hall, and H. Wu, Generation of squeezed states by parametric down conversion, Phys. Rev. Lett. 57, 2520 (1986).
  * Clerk  _et al._ [2010] A. A. Clerk, M. H. Devoret, S. M. Girvin, F. Marquardt, and R. J. Schoelkopf, Introduction to quantum noise, measurement, and amplification, Rev. Mod. Phys. 82, 1155 (2010).
  * Kim  _et al._ [2002] M. S. Kim, W. Son, V. Bužek, and P. L. Knight, Entanglement by a beam splitter: Nonclassicality as a prerequisite for entanglement, Phys. Rev. A 65, 032323 (2002).
  * Azuma  _et al._ [2024] H. Azuma, W. J. Munro, and K. Nemoto, Heralded single-photon source based on superpositions of squeezed states, Phys. Rev. A 109, 053711 (2024).
  * Grangier  _et al._ [1987] P. Grangier, R. E. Slusher, B. Yurke, and A. LaPorta, Squeezed-light enhanced polarization interferometer, Phys. Rev. Lett. 59, 2153 (1987).
  * Aasi  _et al._ [2013] J. Aasi _et al._ , Enhanced sensitivity of the ligo gravitational wave detector by using squeezed states of light, Nature Photonics 7, 613 (2013).
  * Furusawa  _et al._ [2015] A. Furusawa, G. Masada, Kazunori, _et al._ , Continuous-variable entanglement on a chip, Nature Photonics 9, 316 (2015).
  * Asavanant and Furusawa [2024] W. Asavanant and A. Furusawa, Multipartite continuous-variable optical quantum entanglement: Generation and application, Phys. Rev. A 109, 040101 (2024).
  * Vendromin and Dignam [2020] C. Vendromin and M. M. Dignam, Optimization of a lossy microring resonator system for the generation of quadrature-squeezed states, Phys. Rev. A 102, 023705 (2020).
  * Yang  _et al._ [2007] Z. Yang, P. Chak, A. D. Bristow, H. M. van Driel, R. Iyer, J. S. Aitchison, A. L. Smirl, and J. E. Sipe, Enhanced second-harmonic generation in algaas microring resonators, Opt. Lett. 32, 826 (2007).
  * Dutt  _et al._ [2016] A. Dutt, S. Miller, K. Luke, J. Cardenas, A. L. Gaeta, P. Nussenzveig, and M. Lipson, Tunable squeezing using coupled ring resonators on a silicon nitride chip, Opt. Lett. 41, 223 (2016).
  * Yurke [1984] B. Yurke, Use of cavities in squeezed-state generation, Phys. Rev. A 29 (1984).
  * Seifoory  _et al._ [2017] H. Seifoory, S. Doutre, M. M. Dignam, and J. E. Sipe, Squeezed thermal states: the result of parametric down conversion in lossy cavities, Journal of the Optical Society of America B 34, 1587 (2017).
  * Schneider  _et al._ [1988] K. Schneider, M. Lang, J. Mlynek, S. Schiller, F. Marin, A. Bramati, E. Giacobino, T. c Zhang, J.-P. Poizat, J. f Roch, and P. Grangier, Generation of strongly squeezed continuous-wave light at 1064 nm, Opt. Express 2, 59 (1988).
  * Wu  _et al._ [1987] L.-A. Wu, M. Xiao, and H. J. Kimble, Squeezed states of light from an optical parametric oscillator, J. Opt. Soc. Am. B 4 (1987).
  * Aoki  _et al._ [2006] T. Aoki, G. Takahashi, and A. Furusawa, Squeezing at 946nm with periodically poled ktipo(4), Optics Epress 14, 6930 (2006).
  * Takeno  _et al._ [2007] Y. Takeno, M. Yukawa, H. Yonezawa, and A. Furusawa, Observation of -9 db quadrature squeezing with improvement of phase stability in homodyne measurement, Opt. Express 15, 4321 (2007).
  * Vahlbruch  _et al._ [2008] H. Vahlbruch, M. Mehmet, S. Chelkowski, B. Hage, A. Franzen, N. Lastzka, S. Goßler, K. Danzmann, and R. Schnabel, Observation of squeezed light with 10-db quantum-noise reduction, Phys. Rev. Lett. 100, 033602 (2008).
  * Mehmet  _et al._ [2011] M. Mehmet, S. Ast, T. Eberle, S. Steinlechner, H. Vahlbruch, and R. Schnabel, Quantum optics; (270.6570) squeezed states; (190.4970) parametric oscillators and amplifiers, Optics Epress 19, 25763 (2011).
  * Caves  _et al._ [1980] C. M. Caves, K. S. Thorne, R. W. P. Drever, V. D. Sandberg, and M. Zimmermann, On the measurement of a weak classical force coupled to a quantum-mechanical oscillator. i. issues of principle, Rev. Mod. Phys. 52, 341 (1980).
  * Ruskov  _et al._ [2005] R. Ruskov, K. Schwab, and A. N. Korotkov, Squeezing of a nanomechanical resonator by quantum nondemolition measurement and feedback, Phys. Rev. B 71, 235407 (2005).
  * Lei  _et al._ [2016] C. U. Lei, A. J. Weinstein, J. Suh, E. E. Wollman, A. Kronwald, F. Marquardt, A. A. Clerk, and K. C. Schwab, Quantum nondemolition measurement of a quantum squeezed state beyond the 3 db limit, Phys. Rev. Lett. 117, 100801 (2016).
  * Kronwald  _et al._ [2013] A. Kronwald, F. Marquardt, and A. A. Clerk, Arbitrarily large steady-state bosonic squeezing via dissipation, Phys. Rev. A 88, 063833 (2013).
  * Dassonneville  _et al._ [2021] R. Dassonneville, R. Assouly, T. Peronnin, A. Clerk, A. Bienfait, and B. Huard, Dissipative stabilization of squeezing beyond 3 db in a microwave mode, PRX Quantum 2, 020323 (2021).
  * Drummond  _et al._ [1981] P. Drummond, K. McNeil, and D. Walls, Non-equilibrium transitions in sub/second harmonic generation, Optica Acta: International Journal of Optics 28, 211 (1981).
  * Anwar and Zubairy [1992] J. Anwar and M. S. Zubairy, Effect of squeezing on the degenerate parametric oscillator, Phys. Rev. A 45, 1804 (1992).
  * Collett and Gardiner [1984] M. J. Collett and C. W. Gardiner, Squeezing of intracavity and traveling-wave light fields produced in parametric amplification, Phys. Rev. A 30, 1386 (1984).
  * Dunlop  _et al._ [2006] A. E. Dunlop, E. H. Huntington, C. C. Harb, and T. C. Ralph, Generation of a frequency comb of squeezing in an optical parametric oscillator, Phys. Rev. A 73, 013817 (2006).
  * Jabri and Eleuch [2019] H. Jabri and H. Eleuch, Perfect squeezing of terahertz light by two quantum wells using a squeezed vacuum reservoir, J. Opt. Soc. Am. B 36, C1 (2019).
  * Breuer and Petruccione [2007] H.-P. Breuer and F. Petruccione,  _The Theory of Open Quantum Systems_ (Oxford University Press, 2007).
  * Carmicheal [1999] H. J. Carmicheal,  _Statistical Methods in Quantum Optics_ (Springer Berlin, Heidelberg, 1999).
  * Note [1] We can see from Eq. 19 that these properties arise not just for an initial state in equilibrium with the thermal bath, but any initial thermal population that satisfies $n_{th}(0)=n_{b}+a(2n_{b}+1)$. The dynamics of the squeezing amplitude will then be independent of the bath, but the decay will scale by the arbitrary factor $a$.
  * Kim  _et al._ [1989] M. S. Kim, F. A. M. de Oliveira, and P. L. Knight, Properties of squeezed number states and squeezed thermal states, Phys. Rev. A 40, 2494 (1989).
  * Fearn and Collett [1988] H. Fearn and M. J. Collett, Representations of squeezed states with thermal noise, Journal of Modern Optics 35, 553 (1988).
  * Note [2] Direct comparison between the theoretical steady-state value in Eq. 37 and the maximum of the numeric simulation shows a difference of less than $10^{-8}$ for the region $g_{0}<2n_{b}$, which can be attributed to the limits of the computational precision.
  * Constanda [2017] C. Constanda,  _Differential equations a primer for scientists and engineers_ , 2nd ed. (Springer Cham, 2017).
  * Han  _et al._ [2021] K. Han, Y. Wang, and G.-Q. Zhang, Enhancement of microwave squeezing via parametric down-conversion in a superconducting quantum circuit, Opt. Express 29, 13451 (2021).
  * Moon and Girvin [2005] K. Moon and S. M. Girvin, Theory of microwave parametric down-conversion and squeezing using circuit qed, Phys. Rev. Lett. 95, 140504 (2005).

## Supplementary Material: Solving the squeezed thermal state equations

In this supplementary material, we use the technique of Hossein et. al. [16]
to solve equation Eq. 8 with the LME Eq. 2, including the effect of the
environmental temperature. We begin by taking the derivative of Eq. 8 and use
the chain rule to split it in five parts:

$\displaystyle\dot{O}(t)=\dot{O}_{T}(t)+\dot{O}_{S}(t)+\dot{O}_{0}+\dot{O}_{V}+\dot{O}_{L}.$
(S1)

The first four terms, containing contributions from the thermal state density
operator, the squeeze operator, the unperturbed Hamiltonian, and the pump
laser are identical to zero-temperature case, derived in the Appendix of Ref.
[16]. Defining $x\equiv e^{-\hbar\omega/kT},s\equiv\sinh(u),$ and
$c\equiv\cosh(u)$, these expressions for the first four terms can be shown to
be

$\displaystyle\dot{O}_{T}=\anticommutator{J}{O},$ (S2)

where

$J\equiv\frac{d\rho_{T}^{-1/2}}{dt}\rho_{T}^{1/2}=\frac{1}{2x}\frac{dx}{dt}(n_{th}-b^{\dagger}b);$
(S3) $\dot{O}_{S}=\anticommutator{M}{O}+i\commutator{N}{O},$ (S4)

where

$\displaystyle\begin{aligned}
M+iN=&(-is^{2}\dot{\phi})(b^{\dagger}b+\frac{1}{2})\\\
+&\frac{1}{2}\dot{u}(x^{-1}{b^{\dagger}}^{2}e^{i\phi}-xb^{2}e^{-i\phi})\\\
+&\frac{1}{2}ics\dot{\phi}(x^{-1}{b^{\dagger}}^{2}e^{i\phi}+xb^{2}e^{-i\phi});\end{aligned}$
(S5) $\displaystyle\frac{\dot{O}_{0}}{-i\omega}=GO-
OG^{\dagger}=\commutator{P}{O}+i\anticommutator{Q}{O},$ (S6)

where

$\displaystyle\begin{aligned}
G&\equiv\rho_{T}^{-1/2}S^{\dagger}b^{\dagger}bS\rho_{T}^{1/2}\\\
&=s^{2}+(c^{2}+s^{2})b^{\dagger}b-cs(x^{-1}{b^{\dagger}}^{2}e^{i\phi}+xb^{2}e^{-i\phi}),\end{aligned}$
(S7)

thus

$\displaystyle\begin{aligned} P\equiv&s^{2}+(c^{2}+s^{2})b^{\dagger}b\\\
-&\frac{1}{2}cs(x^{-1}+x)({b^{\dagger}}^{2}e^{i\phi}+b^{2}e^{-i\phi}),\end{aligned}$
(S8) $\displaystyle\begin{aligned}
Q\equiv\frac{1}{2}ics(x^{-1}-x)({b^{\dagger}}^{2}e^{i\phi}-b^{2}e^{-i\phi});\end{aligned}$
(S9)

and

$\displaystyle\dot{O}_{V}=-\frac{i}{\hbar}\commutator{\bar{P}}{O}+\frac{1}{\hbar}\anticommutator{\bar{Q}}{O},$
(S10)

where

$\displaystyle\begin{aligned} \bar{P}\equiv&-cs(\alpha\gamma
e^{-i\phi}+\alpha^{*}\gamma^{*}e^{i\phi})\\\ &-2cs(\alpha\gamma
e^{-i\phi}+\alpha^{*}\gamma^{*}e^{i\phi})b^{\dagger}b\\\
&+\frac{1}{2}(\alpha\gamma(x^{-1}+x)c^{2}+\alpha^{*}\gamma^{*}(x^{-1}+x)s^{2}e^{2i\phi}){b^{\dagger}}^{2}\\\
&+\frac{1}{2}(\alpha\gamma(x^{-1}+x)s^{2}e^{-i\phi}+\alpha^{*}\gamma^{*}(x^{-1}+x)c^{2})b^{2}\end{aligned}$
(S11)

and

$\displaystyle\begin{aligned}
\bar{Q}=&-\frac{i}{2}(x^{-1}-x)(\alpha^{*}\gamma^{*}s^{2}e^{2i\phi}+\alpha\gamma
c^{2}){b^{\dagger}}^{2}\\\ &+\frac{i}{2}(x^{-1}-x)(\alpha\gamma
s^{2}e^{2i\phi}+\alpha^{*}\gamma^{*}c^{2}){b^{2}}.\end{aligned}$ (S12)

The final term $\dot{O}_{L}$ has new components related to the bath
population. We see from the master equation that

$\displaystyle\begin{aligned}
\dot{O}_{L}=\Gamma(n_{b}+1)\Big{(}&\rho_{T}^{-1/2}S^{\dagger}b\rho
b^{\dagger}S\rho_{T}^{-1/2}\\\
&-\frac{1}{2}\rho_{T}^{-1/2}S^{\dagger}(\anticommutator{b^{\dagger}b}{\rho})S\rho_{T}^{-1/2}\Big{)}\\\
+\Gamma n_{b}\Big{(}&\rho_{T}^{-1/2}S^{\dagger}b^{\dagger}\rho
bS\rho_{T}^{-1/2}\\\
&-\frac{1}{2}\rho_{T}^{-1/2}S^{\dagger}(\anticommutator{bb^{\dagger}}{\rho})S\rho_{T}^{-1/2}\Big{)}.\end{aligned}$
(S13)

Now, defining the operators

$\displaystyle T=\rho_{T}^{-1/2}S^{\dagger}bS\rho_{T}^{1/2},$ (S14)
$\displaystyle\tilde{T}=\rho_{T}^{-1/2}S^{\dagger}b^{\dagger}S\rho_{T}^{1/2},$
(S15)

and solving

$\displaystyle\begin{aligned}
\tilde{G}&\equiv\rho_{T}^{-1/2}S^{\dagger}bb^{\dagger}S\rho^{1/2}\\\
&=G+\rho_{T}^{-1/2}S^{\dagger}\commutator{b}{b^{\dagger}}S\rho_{T}^{1/2}=G+1,\end{aligned}$
(S16)

Eq. S13 becomes

$\displaystyle\begin{aligned}
\dot{O}_{L}=&\Gamma(n_{b}+1)\left[TOT^{\dagger}-\frac{1}{2}(GO+OG^{\dagger})\right]\\\
+&\Gamma
n_{b}\left[\tilde{T}O{\tilde{T}}^{\dagger}-\frac{1}{2}(\tilde{G}O+O{\tilde{G}}^{\dagger})\right]\\\
=&\Gamma(n_{b}+1)\left[TOT^{\dagger}-\frac{1}{2}(GO+OG^{\dagger})\right]\\\
+&\Gamma
n_{b}\left[\tilde{T}O{\tilde{T}}^{\dagger}-\frac{1}{2}(GO+OG^{\dagger})-O\right]\\\
=&\Gamma(n_{b}+1)\left[TOT^{\dagger}-\frac{1}{2}\anticommutator{P}{O}-\frac{i}{2}\commutator{Q}{O}\right]\\\
&+\Gamma
n_{b}\left[\tilde{T}O{\tilde{T}}^{\dagger}-\frac{1}{2}\anticommutator{P}{O}-\frac{i}{2}\commutator{Q}{O}-O\right].\end{aligned}$
(S17)

Inserting Eqs. S2, S4, S6, S10 and S17 into Eq. S1 and setting $O$ to the
identity for all time yields

$\displaystyle\begin{aligned} 0&=2J+2M+2\omega Q+\frac{2}{\hbar}\bar{Q}\\\
&+\Gamma(n_{b}+1)(TT^{\dagger}-P)+\Gamma
n_{b}(\tilde{T}{\tilde{T}}^{\dagger}-P-1).\end{aligned}$ (S18)

solving for the remaining terms,

$\displaystyle
TT^{\dagger}=xc^{2}+(xc^{2}+x^{-1}s^{2})b^{\dagger}b-cs(b^{2}e^{-i\phi}+{b^{\dagger}}^{2}e^{i\phi})$
$\displaystyle\tilde{T}{\tilde{T}}^{\dagger}=x^{-1}c^{2}b^{\dagger}b-cse^{i\phi}{b^{\dagger}}^{2}-cse^{-i\phi}b^{2}+xs^{2}(b^{\dagger}b+1)$

we use Eqs. S3, S5, S9 and S12 in Eq. S18:

$\displaystyle\begin{aligned}
0&=\frac{1}{x}\frac{dx}{dt}(n_{th}-b^{\dagger}b)+\frac{1}{2}\dot{u}(x^{-1}-x)({b^{\dagger}}^{2}e^{i\phi}+b^{2}e^{-i\phi})\\\
&+\frac{1}{2}ics\dot{\phi}(x^{-1}-x)({b^{\dagger}}^{2}e^{i\phi}-b^{2}e^{-i\phi})\\\
&+i\omega cs(x^{-1}-x)({b^{\dagger}}^{2}e^{i\phi}-b^{2}e^{-i\phi})\\\
&+\frac{i}{\hbar}(x^{-1}-x)(\gamma\alpha
s^{2}e^{-2i\phi}+\gamma^{*}\alpha^{*}c^{2})b^{2}\\\
&-\frac{i}{\hbar}(x^{-1}-x)(\gamma^{*}\alpha^{*}s^{2}e^{2i\phi}+\gamma\alpha
c^{2}){b^{\dagger}}^{2}\\\ &\begin{aligned}
+\Gamma(n_{b}+1)\Big{(}xc^{2}&+(xc^{2}+x^{-1}s^{2})b^{\dagger}b\\\
&-cs(b^{2}e^{-i\phi}+{b^{\dagger}}^{2}e^{i\phi})\Big{)}\end{aligned}\\\
&\begin{aligned} -\Gamma(n_{b}+1)\Big{(}s^{2}&+(c^{2}+s^{2})b^{\dagger}b\\\
&-\frac{1}{2}cs(x^{-1}+x)({b^{\dagger}}^{2}e^{i\phi}+b^{2}e^{-i\phi})\Big{)}\end{aligned}\\\
&\begin{aligned} +\Gamma
n_{b}\Big{(}xs^{2}&+(x^{-1}c^{2}+xs^{2})b^{\dagger}b\\\
&-cs(e^{i\phi}{b^{\dagger}}^{2}+e^{-i\phi}b^{2})\Big{)}\end{aligned}\\\
&\begin{aligned} -\Gamma n_{b}\Big{(}s^{2}&+(c^{2}+s^{2})b^{\dagger}b\\\
&-\frac{1}{2}cs(x^{-1}+x)({b^{\dagger}}^{2}e^{i\phi}+b^{2}e^{-i\phi})+1\Big{)}.\end{aligned}\end{aligned}$
(S19)

We now introduce the Hermitian operators

$\displaystyle\chi_{1}={b^{\dagger}}^{2}e^{i\phi}+{b}^{2}e^{-i\phi}$ (S20)
$\displaystyle\chi_{2}=i({b^{\dagger}}^{2}e^{i\phi}-{b}^{2}e^{-i\phi}).$ (S21)

Substituting these into Eq. S19, we obtain

$\displaystyle\begin{aligned}
0&=\frac{1}{x}\frac{dx}{dt}(n_{th}-b^{\dagger}b)+\frac{1}{2}\dot{u}(x^{-1}-x)\chi_{1}\\\
&+\frac{1}{2}cs\dot{\phi}(x^{-1}-x)\chi_{2}+\omega cs(x^{-1}-x)\chi_{2}\\\
&+\frac{i}{\hbar}(x^{-1}-x)(\gamma\alpha
s^{2}e^{-i\phi}+\gamma^{*}\alpha^{*}c^{2}e^{i\phi})\frac{1}{2}(\chi_{1}+i\chi_{2})\\\
&-\frac{i}{\hbar}(x^{-1}-x)(\gamma^{*}\alpha^{*}s^{2}e^{i\phi}+\gamma\alpha
c^{2}e^{-i\phi})\frac{1}{2}(\chi_{1}-i\chi_{2})\\\
&+\Gamma(n_{b}+1)\left(xc^{2}+(xc^{2}+x^{-1}s^{2})b^{\dagger}b-cs\chi_{1}\right)\\\
&-\Gamma(n_{b}+1)\left(s^{2}+(c^{2}+s^{2})b^{\dagger}b-\frac{1}{2}cs(x^{-1}+x)\chi_{1}\right)\\\
&+\Gamma
n_{b}\left(xs^{2}+(x^{-1}c^{2}+xs^{2})b^{\dagger}b-cs\chi_{1}\right)\\\
&-\Gamma
n_{b}\left(s^{2}+(c^{2}+s^{2})b^{\dagger}b-\frac{1}{2}cs(x^{-1}+x)\chi_{1}\right)\\\
&-\Gamma n_{b}.\end{aligned}$ (S22)

This allows us to form four independent equations using

$\displaystyle F_{1}\chi_{1}+F_{2}\chi_{2}+F_{3}b^{\dagger}b+F_{4}=0,$ (S23)

where

$\displaystyle\begin{aligned}
F_{1}&=\frac{1}{2}\dot{u}(x^{-1}-x)+\frac{i}{2\hbar}(x^{-1}-x)(\gamma^{*}\alpha^{*}e^{i\phi}-\gamma\alpha
e^{-i\phi})\\\
&+\Gamma(2n_{b}+1)cs\big{(}\frac{1}{2}(x^{-1}+x)-1\big{)},\end{aligned}$ (S24)
$\displaystyle\begin{aligned} F_{2}&=\frac{1}{2}cs\dot{\phi}(x^{-1}-x)+\omega
cs(x^{-1}-x)\\\
&-\frac{1}{2\hbar}(x^{-1}-x)(c^{2}+s^{2})(\gamma^{*}\alpha^{*}e^{i\phi}+\gamma\alpha
e^{-i\phi}),\end{aligned}$ (S25) $\displaystyle\begin{aligned}
F_{3}=&-\frac{1}{x}\frac{dx}{dt}+\Gamma(n_{b}+1)(xc^{2}+x^{-1}s^{2}-c^{2}-s^{2})\\\
&+\Gamma n_{b}(x^{-1}c^{2}+xs^{2}-c^{2}-s^{2}),\end{aligned}$ (S26)

and

$\displaystyle\begin{aligned} F_{4}=&\frac{1}{x}\frac{dx}{dt}n_{th}\\\
&+\Gamma(n_{b}+1)(xc^{2}-s^{2})+\Gamma n_{b}(xs^{2}-s^{2}-1).\end{aligned}$
(S27)

The equations for $F_{3}$ and $F_{4}$ both generate equations in $x$, which
are consistent and can be reduced to

$\displaystyle\frac{dx}{dt}=-(1-x)\Gamma\left[n_{b}(2s^{2}+1)(x-1)+s^{2}(x-1)+x\right].$
(S28)

Using $\frac{dn_{th}}{dt}=\frac{1}{(1-x)^{2}}\frac{dx}{dt}$, this translates
to the following equation for $n_{th}$:

$\displaystyle\begin{aligned}
\frac{dn_{th}}{dt}=&\Gamma\left[n_{b}(2s^{2}+1)+s^{2}-\frac{x}{1-x}\right]\\\
=&\Gamma\left[n_{b}(2s^{2}+1)+s^{2}-n_{th}\right],\end{aligned}$ (S29)

which is equivalent to Eq. 9. Meanwhile, from $F_{2}=0$, we obtain

$\displaystyle\dot{\phi}=-2\omega+\frac{1}{\hbar}\frac{c^{2}+s^{2}}{cs}(\gamma^{*}\alpha^{*}e^{i\phi}+\gamma\alpha
e^{-i\phi}),$ (S30)

which is Eq. 11. Finally, from $F_{1}=0$, we obtain

$\displaystyle\dot{u}=-\frac{i}{\hbar}(\gamma^{*}\alpha^{*}e^{i\phi}-\gamma\alpha
e^{-i\phi})-\frac{\Gamma(2n_{b}+1)cs}{2n_{th}+1},$ (S31)

which is simply Eq. 10.


# Geometric algorithms for predicting resilience and recovering damage in
neural networks

Guruprasad Raghavan
Caltech
Pasadena, CA 91106
<EMAIL_ADDRESS>
&Jiayi Li
UCLA
Los Angeles, CA 90095
<EMAIL_ADDRESS>&Matt Thomson
Caltech
Pasadena, CA 91106
<EMAIL_ADDRESS>

###### Abstract

Biological neural networks have evolved to maintain performance despite
significant circuit damage. To survive damage, biological network
architectures have both intrinsic resilience to component loss and also
activate recovery programs that adjust network weights through plasticity to
stabilize performance. Despite the importance of resilience in technology
applications, the resilience of artificial neural networks is poorly
understood, and autonomous recovery algorithms have yet to be developed. In
this paper, we establish a mathematical framework to analyze the resilience of
artificial neural networks through the lens of differential geometry. Our
geometric language provides natural algorithms that identify local
vulnerabilities in trained networks as well as recovery algorithms that
dynamically adjust networks to compensate for damage. We reveal striking
weight perturbation vulnerabilities in common image analysis architectures,
including MLP’s and CNN’s trained on MNIST and CIFAR-10 respectively. We also
uncover high-performance recovery paths that enable the same networks to
dynamically re-adjust their parameters to compensate for damage. Broadly, our
work provides procedures that endow artificial systems with resilience and
rapid-recovery routines to enable their deployment for critical applications.

##  1 Introduction

Brains are remarkable machines whose computational capabilities have inspired
many breakthroughs in machine learning [1, 2, 3, 4]. However, the resilience
of the brain, its ability to maintain computational capabilities in harsh
conditions and following circuit damage, remains poorly developed in current
artificial intelligence paradigms [5] . Biological neural networks are known
to implement redundancy and other architectural features that allow circuits
to maintain performance following loss of neurons or lesion to sub-circuits
[6, 7, 8, 9, 10]. In addition to architectural resilience, biological neural
networks execute recovery programs that allow circuits to repair themselves
through the activation of network plasticity following damage [11, 12, 13].
For example, recovery algorithms reestablish olfactory and visual behaviors in
mammals following sensory specific cortical circuit lesions [14, 15]. Through
resilience and recovery mechanisms, biological neural networks can maintain
steady performance in the face of dynamic challenges like changing external
environments, cell damage, partial circuit loss as well as catastrophic
injuries like the loss of large sections of the cortex. [16, 17, 18, 19].

Like brains, artificial neural networks must increasingly execute critical
applications that require robustness to both hardware component damage and
memory errors that could corrupt network weights. Recent studies have
highlighted the importance of network robustness to soft errors that can lead
to weight corruption and network failure [20] in applications including (i)
decision-making in the healthcare industry, (ii) image and sensor analysis in
self-driving cars and (iii) robotic control systems. Errors in dynamic access
memory can occur due to malicious attacks (the RowHammer), but a particular
focus has been on errors induced by high energy particles [21] that occur at
surprising rates [22]. Further, the rising implementation of neural networks
on physical hardware (like neuromorphic, edge devices) [23, 24], where
networks can be disconnected from the internet and are under control of an end
user, necessitates the need for damage-resilient and dynamically recovering
artificial neural networks.

The resilience of living neural networks motivates theoretical and practical
efforts to understand the resilience of artificial neural networks and to
design new algorithms that reverse engineer resilience and recovery into
artificial systems [25]. Recent studies [26, 27] have demonstrated empirically
that MLP and CNN architectures can be surprisingly robust to large scale node
deletion. However, there is currently little understanding of the empirically
observed resilience or what ultimately causes networks to fail. Mathematical
frameworks will be important for understanding the resilience neural networks
and for developing recovery procedures that can maintain network performance
during damage.

We propose a mathematical framework grounded in differential geometry for
studying the resilience and the recovery of artificial neural nets. We
formalize damage/response behavior as dynamic movement on a curved pseudo-
Riemannian manifold. Our geometric language provides new procedures for
identifying network vulnerabilities by predicting local perturbations that
adversely impact the functional performance of the network. Further, we
demonstrate that geodesics, minimum length paths, on the weight manifold
provide high performance recovery paths that the network can traverse to
maintain performance while damaged. Our algorithms allow networks to maintain
high-performance during rounds of damage and repair through computationally
efficient weight-update algorithms that do not require conventional
retraining. Broadly, our work provides procedures that will help endow
artificial systems with resilience and autonomous recovery policies to emulate
the properties of biological neural networks.

##  2 Analyzing network resilience with differential geometry

We develop a geometric framework for understanding how artificial neural
networks respond to damage using differential geometry to analyze changes in
functional performance given changes in network weights. Two recent papers
have highlighted intrinsic robustness properties of layered neural networks
[26, 27]. We provide a geometric approach for understanding robustness as
arising from underlying geometric properties of the weight manifold that are
quantified by the metric tensor. The geometric approach allows us to identify
vulnerabilities in common neural network architectures as well as define new
strategies for repairing damaged networks.

We represent a feed-forward neural network as a smooth,
$\mathbb{C}^{\infty}$function $f(\mathbf{x},\mathbf{w})$, that maps an input
vector, $\mathbf{x}\in\mathbb{R}^{\text{k}}$, to an output vector,
$f(\mathbf{x},\mathbf{w})=\mathbf{y}\in\mathbb{R}^{\text{m}}$. The function,
$f(\mathbf{x},\mathbf{w})$, is parameterized by a vector of weights,
$\mathbf{w}\in\mathbb{R}^{\text{n}}$, that are typically set in training to
solve a specific task. We refer to $W=\mathbb{R}^{n}$ as the weight space
($W$) of the network, and we refer to $\mathcal{F}=\mathbb{R}^{m}$ as the
functional manifold [28]. In addition to $f$, we will sometimes be interested
in considering a loss function,
$L:\mathbb{R}^{\text{m}}\times\mathbb{R}\rightarrow\mathbb{R}$, that provides
a scalar measure of network performance for a given task (Figure 1).

We ask how the performance of a trained neural network, $\mathbf{w_{t}}$, will
change when subjected to weight perturbation, shifting
$\mathbf{w_{trained}}\rightarrow\mathbf{w_{damaged}}$. We use differential
geometry to develop a mathematical theory, rooted in a functional notion of
distance, to analyze how arbitrary weight perturbations
$\mathbf{w_{t}}\rightarrow\mathbf{w_{d}}$ impact functional performance of a
network. Specifically, we construct a local distance metric, $\mathbf{g}$,
that can be applied at any point in $W$ to measure the functional impact of an
arbitrary network perturbation.

To construct a metric mathematically, we fix the input, $\mathbf{x}$, into a
network and ask how the output of the network, $f(\mathbf{w},\mathbf{x})$,
moves on the functional manifold, $\mathcal{F}$, given an infinitesimal weight
perturbation, $\mathbf{du}$, in $W$ where
$\mathbf{w_{d}}=\mathbf{w_{t}}+\mathbf{du}$. For an infinitesimal perturbation
$\mathbf{du}$,

$f(\mathbf{x},\mathbf{w_{t}}+\mathbf{du})\approx
f(\mathbf{x},\mathbf{w_{t}})+\mathbf{J_{w_{t}}}\ \mathbf{du},$ (1)

where $\mathbf{J_{w_{t}}}$ is the Jacobian of $f(\mathbf{x},\mathbf{w})$ for a
fixed $\mathbf{x}$, $J_{i,j}=\frac{\partial f_{i}}{\partial w^{j}}$, evaluated
at $\mathbf{w_{t}}$. We measure the change in functional performance given
$\mathbf{du}$ as the mean squared error

$\displaystyle
d(\mathbf{w_{t}},\mathbf{w_{d}})=|f(\mathbf{x},\mathbf{w_{t}})-f(\mathbf{x},\mathbf{w_{d}})|^{2}$
$\displaystyle=\mathbf{du}^{T}\ (\mathbf{J_{w_{t}}}^{T}\ \mathbf{J_{w_{t}}})\
\mathbf{du}$ (2) $\displaystyle=\mathbf{du}^{T}\ \mathbf{g_{w_{t}}}\
\mathbf{du},$ (3)

where $\mathbf{g_{w_{t}}}=\mathbf{J_{w_{t}}}^{T}\mathbf{J_{w_{t}}}$ is the
metric tensor evaluated at the point $\mathbf{w_{t}}\in W$. The metric tensor
$\mathbf{g}$ is an $n\times n$ symmetric matrix that defines an inner product
and local distance metric,
$\langle\mathbf{du},\mathbf{du}\rangle_{\mathbf{w}}=\mathbf{du^{T}}\
\mathbf{g_{w}}\ \mathbf{du}$, on the tangent space of the manifold, $T_{w}(W)$
at each $\mathbf{w}\in W$.

Figure 1: Geometric framework for analyzing neural network resilience (A)
Three networks ($N_{1},N_{2},N_{3}$) in weights space $W$ and their relative
distance in functional space and loss space. Damage is analyzed by asking how
movement in weight space changes functional performance and loss through
introduction of a pullback metric $\mathbf{g}$. (B) We consider local damage
to a network as an infinitesimal perturbation that can be analyzed in the
tangent space of a trained network. (C) Global damage is modeled as long range
movement of network weights along a path, $\gamma(t)$, in weight space.

Explicitly,

$g_{ij}=\sum_{k=1}^{m}\frac{\partial
f_{k}(\textbf{x},\mathbf{w})}{\partial\mathbf{w^{i}}}\frac{\partial
f_{k}(\textbf{x},\mathbf{w})}{\partial\mathbf{w^{j}}},$ (5)

where the partial derivatives $\frac{\partial
f_{k}(\textbf{x},\mathbf{w})}{\partial\mathbf{w^{i}}}$ measure change in
functional output of a network given a change in weight. In the appendix, we
extend the metric formulation to cases where we consider a set, $\mathbf{X}$,
of training data and view $\mathbf{g}$ as the average of metrics derived from
individual training examples. The metric, $\mathbf{g}$, provides a local
measure of functional distance on the pseudo-Riemmanian manifold
$(W,\mathbf{g})$. At each point in weight space, the metric defines the
length, $\langle\mathbf{du},\mathbf{du}\rangle_{\mathbf{w}}$, of a local
perturbation by its impact on the functional output of the network (Figure
1b).

Globally, we can use the metric to determine the functional performance change
across a path connected set of networks. Mathematically, the metric changes as
we move in $W$ due to the curvature of the ambient space that reflects changes
in the vulnerability of a network to weight perturbation (Figure 1c). As a
network moves along a path, $\gamma(t)\in W$ from a given trained network
$\gamma(0)=\mathbf{w_{t}}$ to a damaged network $\gamma(1)=\mathbf{w_{d}}$, we
can analyze the integrated impact of damage on network performance along
$\gamma(t)$ by using the metric to calculate the length of the path
$\gamma(t)$ as:

$L(\gamma)=\int^{1}_{0}{\langle\frac{d\gamma(t)}{dt},\frac{d\gamma(t)}{dt}\rangle_{\gamma(t)}}\
dt,$ (6)

where
$\langle\frac{d\gamma(t)}{dt},\frac{d\gamma(t)}{dt}\rangle_{\gamma(t)}=\frac{d\gamma(t)}{dt}^{T}\mathbf{g}_{\gamma(t)}\frac{d\gamma(t)}{dt}$
is the infinitesimal functional change accrued while traversing path
$\gamma(t)\in W$.

In what follows, we study the resilience of neural networks by analyzing the
structure of the metric tensor along paths in weight space. We show that the
metric tensor can be used to develop recovery procedures by finding ‘geodesic
paths’, minimum length paths, in the pseudo-Riemannian manifold that allow
networks to respond to damage while suffering minimal performance degradation.

##  3 The geometry of local damage and network vulnerability

We, first, apply our mathematical framework to analyze the response of trained
neural networks to small, local weight perturbations. Empirical studies have
demonstrated that trained networks are often robust to small, local weight
perturbation [26, 27]. We connect local resilience to the spectral properties
of the metric tensor, $\mathbf{g}$, at a given position, $\mathbf{w_{t}}$, in
weight space. We find that networks are typically robust to random local
weight perturbations but also have catastrophic vulnerabilities to specific
low magnitude weight perturbations that dramatically alter network
performance.

To understand local damage, we consider a trained network, $\mathbf{w_{t}}$,
and we subject the network to an infinitesimal weight perturbations in a
direction $\mathbf{du}=c_{i}\ \mathbf{dw^{i}}$ yielding the perturbed weights
$\mathbf{w^{\prime}}=\mathbf{w_{t}}+\mathbf{du}$. We use $\mathbf{dw^{i}}$ to
indicate an infinitesimal displacement vector in the direction $w_{i}$.
Formally, we view $\mathbf{du}$ as a vector in the tangent space of W at
$\mathbf{w_{t}}$, $T_{w_{t}}(W)$ (Figure 1B). The metric tensor, evaluated at
the point $\mathbf{w_{t}}$ provides a local measure of functional performance
change induced by the perturbation along $\mathbf{du}$ through Equation 8.

Figure 2: Metric tensor explains local resilience and predicts catastrophic
vulnerabilities. (A) Spectra of the metric tensor for MLP-1, MLP-2 and VGG-11
(B,C) Test performance of networks perturbed within a unit-ball in $W$(B)
perturbed VGG-11 trained on CIFAR-10, (C) perturbed MLP-2 trained on MNIST
(D,E) Designing adversarial perturbations to destroy trained networks’
performance. (D) adversarial perturbation within unit-ball in $W$ lowers
accuracy to $13\%$ in VGG-11, (E) adversarial weight perturbation within unit-
ball in $W$ lowers accuracy to $70\%$ in MLP-2.

As a positive semi-definite, symmetric matrix, $\mathbf{g}$ (evaluated at
$\mathbf{w_{t}}$) has an orthonormal eigenbasis $\\{\mathbf{v_{i}}\\}$ with
eigenvalues $\lambda_{i}$, $\lambda_{i}\geq 0$. The eigenvalue $\lambda_{i}$
locally determines how a perturbation along the eigenvector $\mathbf{v_{i}}$
will alter functional performance. Expanding an arbitrary perturbation,
$\mathbf{du}$ in the basis $\\{\mathbf{v_{i}}\\}$, as
$\mathbf{du}=\sum_{i}c_{i}\ \mathbf{v_{i}}$, the functional performance change
of the network is

$\displaystyle d(\mathbf{w_{t}},\mathbf{w_{t}}+\mathbf{du})$
$\displaystyle=\mathbf{du^{T}}\ \mathbf{g_{w_{t}}}\ \mathbf{du}$ (7)
$\displaystyle=\sum_{i}c^{2}_{i}\lambda_{i}$ (8)

where $c_{i}=\langle\mathbf{du},\mathbf{v_{i}}\rangle$ quantifies the
contribution of vector $\mathbf{v_{i}}$ to $\mathbf{du}$. Thus, the
performance change, $d(\mathbf{w_{t}},\mathbf{w_{t}}+\mathbf{du})$, incurred
by a network, following perturbation $\mathbf{du}$ is determined by the
magnitude of each $\lambda_{i}$ and the projection of $\mathbf{du}$ onto
$\mathbf{v_{i}}$. The eigenvalues $\lambda$ convert weight changes into change
in functional performance and so have units of $\frac{\text{performance
change}}{\text{weight change}}$. A network will be resilient to weight
perturbations directed along eigenvectors, $\mathbf{v_{i}}$, with small
eigenvalues ($\lambda_{i}<10^{-3}$). Alternately, networks are vulnerable to
perturbations along directions with larger eigenvalues
($\lambda_{i}>10^{-3}$). Our definition of resilient directions,
$\lambda_{i}<10^{-3}$, is an operational direction that selects directions
where a unit of weight change will produce a performance change of less than
$10^{-3}$ or $.1\%$.

Mathematically, we can understand the resilience of networks to randomly
distributed weight perturbations by calculating the average response of a
network to Gaussian weight perturbations, $\mathbf{du}\sim P(\mathbf{du})$,
where $P(du_{i})=\mathcal{N}(0,\tfrac{\sigma}{d})$ ($n=\text{dim}(W)$ and
$\mathbb{E}[||\mathbf{du}||_{2}]\ =\sigma)$ . The expectation of the induced
performance change for such a Gaussian perturbation is

$\displaystyle\mathbb{E}_{du_{i}\sim\mathcal{N}(0,\tfrac{\sigma}{d})}[d(\mathbf{w},\mathbf{w}+\mathbf{du})]$
$\displaystyle=\frac{\sigma}{d}\sum_{i}\lambda_{i}$ (9) $\displaystyle<\
\sigma\ \rho\ \lambda_{1},$ (10)

where $\rho$ indicates the fraction of vulnerable directions, and
$\lambda_{1}$ is the largest eigenvalue of $\mathbf{g}$.

Empirically, we find that trained networks are, perhaps as expected, robust to
‘random’ local perturbation (Figure 2) due to a large fraction of resilient
eigendirections ($\rho<10^{-3}$). Such local network robustness holds for a
series of trained network architectures including (i) Multi-layer perceptrons
(MLP-1, MLP-2) trained on MNIST and (ii) Convolutional neural networks
(VGG-11) trained on CIFAR-10. MLP-1 is a single hidden layer network, with
variable number of hidden nodes, while MLP-2 is the LeNet architecture
borrowed from [29] (2 hidden layers, with 300 and 100 hidden nodes
respectively). VGG-11 for CIFAR-10 is adapted from [30]111The network
architecture, pre-trained models and optimization algo’s are specified in the
appendix.

Consistent with their eigenspectra (VGG-11: $\rho<10^{-4}$, MLP-1, MLP-2:
$\rho<10^{-3}$) , both MLP and CNN architectures exhibit minimal performance
degradation for unit-ball perturbations 222Unit-ball perturbations, due to the
high dimensionality of the space, induce an average weight change of
$<10^{-6}$ for individual weights  ($\sigma=1$, Figure 2A) . When perturbed
along 1000 directions of unit-norm, the trained MLP-2 (initial test accuracy
of 98$\%$) maintains accuracy of 97.2-97.6$\%$ (Figure 2C). Perturbation of
VGG-11 trained on CIFAR-10 (initial test accuracy of 91$\%$) yields networks
with test accuracy between 88-91$\%$ (Figure 2B).

Resilience to such small local perturbations might be expected, but our
framework also exposes hidden catastrophic vulnerabilities to perturbations of
the same order in both networks. By designing adversarial weight perturbations
to lie along the ‘vulnerable’ eigenvectors of $\mathbf{g}$ ($\mathbf{v_{i}}$
with large $\lambda_{i}$), we can induce sharp performance declines across
architectures (Figure 2D,E). For the VGG-11 network trained on CIFAR-10, an
adversarial weight perturbation decreases accuracy from 91$\%$ to 13$\%$
(Figure 2D). Similarly, adversarial perturbation reduces the performance of
MLP-2 network trained on MNIST from 98$\%$ to 70$\%$ (Figure 2E). For the
CIFAR-10 network, a relatively small perturbation causes the network to make
critical classification errors making the erroneous inference of most CIFAR-10
images to being in the class of ‘automobiles’. In this way, the local geometry
of the weight manifold allows us to discover subtle weight perturbations that
cause catastrophic changes in network performance for small change in network
weights.

##  4 Acceleration identifies global break-down points in a network

From earlier studies [26, 27], we know that trained MLP’s and CNN’s can be
surprisingly robust to much more profound global damage including large scale
node deletion. In this section, we develop a concept of break-down
acceleration using the covariant derivative of a network along paths
connecting the trained network and damaged network in $W$. Break-down
acceleration predicts failure points that emerge in weight space through rapid
changes in the curvature of the weight space, and ultimately allows us to
develop procedures to thwart break-down by avoiding acceleration.

(a)  Figure 3:  Break-down acceleration characterizes network break-down
points following damage Performance of an (A) MLP-1 network (1 hidden-layer,
variable hidden nodes) and (B) VGG-11 during simulated damage to distinct
layers. Both networks experience sharp performance break down when network
damage exceeds (A) $\sim$90$\%$ of hidden-nodes for MLP-1 and (B) $\sim$60$\%$
of nodes in any layer for VGG-11. (C,D,E) Damage paths in manifold ($W,g$).
(C) A cartoon of the loss landscape showing multiple break-down paths from the
trained network to the damaged network (D,E) The covariant derivative of the
accuracy along multiple damage paths for (D) MLP-2 and (E) VGG-11 are shown. A
steep increase in the covariant derivative (acceleration) along damage paths
corresponds to the networks’ sharp break-down to global damage.(F) Multiple
damage paths (colored lines) shown from trained MLP-2 (N1) to its damaged
counter-part (N2). The z-axes is the test-accuracy of the networks, while x,y
axes are the isomap embedding of networks in a 2D space.

Mathematically, we represent global damage as a path in weight space,
$\gamma{(t)}\in W$ with $t\in[0,1]$, that connects a trained network,
$\gamma(0)=\mathbf{w_{t}}$, to its damaged counterpart
$\gamma(1)=\mathbf{w_{d}}$ (Figure 3C). Practically, global damage might
emerge as a discrete event (node deletion), our analysis provides a continuous
approximation to discrete network damage. As a network moves along a path from
$\mathbf{w_{t}}$ to $\mathbf{w_{d}}$, the metric tensor itself changes,
changing its spectra and its vulnerability.

Along a path, $\gamma(t)\in W$ the velocity vector, $v(t)=\frac{d\gamma}{dt}$,
quantifies the change in the functional performance of a network per unit
time. Mathematically, we define the break-down speed ($s$) of a network along
a path in weight space as the norm of the network’s velocity vector computed
using the metric tensor
$s(t)=\langle\frac{d\gamma}{dt},\frac{d\gamma}{dt}\rangle_{\gamma(t)}=\sum_{ij}g_{ij}w_{ti}\
w_{tj}$. Non-linear break-down points emerge along paths in $W$ when break-
down speed undergoes a rapid acceleration, so that $\frac{ds}{dt}>>0$. We can
calculate the break-down speed and acceleration explicitly for a network
following simple straight or Euclidean path from a trained to damaged
configuration. Taking $w_{d}=0$, $\gamma(t)=\mathbf{w_{t}}(1-t)$ and
$\frac{d\gamma(t)}{dt}=-\mathbf{w_{t}}$, we have

$\displaystyle\frac{ds}{dt}=\sum_{i,j}\sum_{k}\frac{dg_{ij}}{dw_{k}}w_{tk}\
w_{ti}\ w_{tj}$ (11)

where $g_{ij}$ is evaluated along $\gamma(t)$. The change in the metric tensor
$\frac{dg_{ij}}{dw_{k}}$ along a path $\gamma(t)$, thus, determines whether
performance decays at a constant $\frac{ds}{dt}=0$, $\frac{dg_{ij}}{dk}=0$ or
at an accelerating $\frac{ds}{dt}>0$, $\frac{dg_{ij}}{dw_{k}}>0$ rate. For
curved paths break-down acceleration can be analyzed using an object known as
the covariant derivative, $\mathbf{\nabla}_{\gamma(t)}v(t)$ (Appendix).

In practice, calculation of the break-down acceleration identifies, damage
failure points in real neural networks. For example, both MLP-1 and VGG-11
architectures tolerate considerable node deletion (Figure 3A). MLP-1, (one
hidden layer, 400 hidden units) trained on MNIST, tolerates damage to 80$\%$
of the network nodes reducing functional performance of the network by merely
$\sim$10$\%$. Similarly, VGG-11 trained on CIFAR-10 tolerates $\sim 60\%$ node
damage in any layer without performance degradation. However, both networks
exhibit drastic break-down in functional performance beyond these node damage
thresholds (Figure 3A,B). Mathematically, break-down points occur where the
acceleration of the network, as measured by the covariant derivative along the
damage path, rapidly increases (Figure 3D-F) Steep increases in the covariant
derivative identify points of loss acceleration corresponding to the
functional breakdown of both the networks analyzed.

##  5 Geodesic paths enable network recovery

Thus, globally, network break-down occurs along a damage path in $W$due to
abrupt changes in the curvature of the underlying functional landscape that
result in abrupt change in the metric. The mathematical connection between
break-down and curvature suggests a strategy for designing recovery protocols
that can adapt a neural network’s weights to compensate for damage. Inspired
by recovery mechanisms in neuroscience that compensate for damage by altering
the weights of undamaged nodes. We apply the concept of break-down
acceleration to develop recovery procedures for artificial neural networks
that compensate for damage through continuous adjustment of the undamaged
weights by minimizing the acceleration along the path.

Mathematically, minimum acceleration paths in weight space are known as
geodesic paths. Geodesic paths, by definition, provide both minimum length and
minimum acceleration paths in weight space. Specifically, we consider a
trained network, $\mathbf{w}$, subjected to weight damage that zeros a subset
of weights, $w_{i}=0$, for $i\in n_{\text{damaged}}$. Our strategy responds to
damage by adjusting undamaged weights, $w_{i}$ for $i\notin
n_{\text{damaged}}$ to maximize network performance by moving the network
along a geodesic in $W$. Geodesic paths can be computed directly using our
metric $\mathbf{g}$ and also represent the minimum distance paths (with
distance defined in equation-6) between two points on $W$.

(a)  Figure 4: Geodesic paths allow damage compensation through weight
adjustment: (A) Test accuracy of geodesic recovery paths (blue) versus naive
damage paths (red) for VGG-1 network while $30$ convolution filters and 1000
nodes from fully-connected layers are damaged. While naive path exhibits sharp
break-down, geodesic procedure adjusts undamaged weights to maintain test
accuracy. (B) Magnitude of the covariant derivative (break-down acceleration)
for geodesic (blue) and naive damage paths (red). C Test accuracy (C) and (D)
number of network update epochs (D) for geodesic recovery (blue) vs fine-
tuning (green) while 50 (out of 60) conv-filters are deleted from layer1 in
VGG11. Geodesic recovery requires $\leq$10 total update epochs .(E) A
depiction of multiple recovery paths on the loss landscape from trained
network (N1) to networks on the damage hyper-plane (N2, N3, N4, N5). The
z-axes is network-loss, while the x,y axes are neural net weights. (F)
Geodesic strategy (blue) allows networks to dynamically transition between
configurations: C1, trained VGG11 network ; C2, 50 conv-filters removed from
C1; C3, 1000 additional nodes removed from classifier-layers in C2; C4, 30
conv-filters in conv-layer1 restored to C3. Dynamic transitioning enabled
within 5 epochs. Naive strategy is red).

In general, geodesic paths are typically calculated using the geodesic
equation (Appendix) an ordinary differential equation that uses derivatives of
the metric tensor to identify minimum acceleration paths in a space given an
initial velocity. However, solutions to the geodesic equation are
computationally prohibitive for large neural networks as they require
evaluation of the Christoffel symbols which scale as a third order polynomial
of the number of parameters in the neural network ($\mathcal{O}(n^{3})$).

Therefore, we developed an approximation to the geodesic equation using a
first order expansion of the loss function. Given a trained network, our
procedure updates the weights of the network to optimize performance given a
direction of damage. To discover a geodesic path $\gamma(t)$, we begin at a
trained network and iteratively solve for the tangent vector , $\theta(w)$, at
every point, $\mathbf{w}=\gamma(t)$, along the path, starting from
$\mathbf{w_{t}}$ and terminating at the damage hyperplane, $W_{d}$. The damage
hyperplane is the set of all networks, $\mathbf{w}\in W$, such that $w_{i}=0$,
for $i\in n_{\text{damaged}}$. We specifically solve

$\text{argmin}_{\theta(\mathbf{w})}\
\langle\theta(\mathbf{w}),\theta(\mathbf{w})\rangle_{\mathbf{w}}-\beta\
\theta(\mathbf{w})^{T}v_{w}\\\ \text{ subject to:
}\theta(\mathbf{w})^{T}\theta(\mathbf{w})\leq 0.01.$ (12)

The tangent vector $\theta(\mathbf{w})$ is obtained by simultaneously
optimizing two objective functions: (1) minimizing the increase in functional
distance along the path measured by the metric tensor ($\mathbf{g_{w}}$) [min:
($\langle\theta(\mathbf{w}),\theta(\mathbf{w})\rangle_{\mathbf{w}}$) =
$(\theta(\mathbf{w})^{T}\mathbf{g_{w}}\theta(\mathbf{w}))$ ] and (2)
maximizing the dot-product between the tangent vector and $v_{\mathbf{w}}$,
vector pointing in the direction of the hyperplane [max:
($\theta(\mathbf{w})^{T}v_{\mathbf{w}}$)] to enable movement towards the
damage hyperplane. By finding geodesic paths to the damage hyperplane, we find
weight adjustments that can be made within a network during damage to maintain
performance (Figure 4E).

Our optimization procedure is a quadratic program that trades off, through the
hyper-parameter $\beta$, motion towards the damage hyper-plane and the
maximization of the functional performance of the intermediate networks along
the path (optimization procedure elaborated in the appendix). The strategy
discovers multiple paths from the trained network $\mathbf{w_{t}}$ to $W_{d}$,
damage hyper-plane, (depicted as path-1 to path-5 in Figure 4a) where networks
maintain high functional performance during damage. Of the many paths
obtained, we can select the path with the shortest total length (with respect
to the metric $\mathbf{g}$) as the best approximation to the geodesic in the
manifold.

The geodesic strategy enables damage compensation through continuously
updating weights in the network. We apply the geodesic strategy to discover
recovery paths from a trained network (VGG11) to a pre-defined damage
hyperplane333Fig-4A, 4B: Damage hyperplane for VGG11 is defined by := deletion
of 30 conv-filters from layer1,2 and 1000 nodes from fully-connected
layer1,2.. The recovery path is a high-performance path with all networks
performing above 87$\%$ test accuracy, and the recovery path maintains low
break-down acceleration when compared to the naive (linear) path (Figure 4B).
A similar analysis for MLP’s is presented in the appendix.

While high-performance paths can also be discovered through heuristic fine-
tuning, the geodesic procedure is both rationale and computationally
efficient. Specifically, an iterative prune-train cycle achieved through
structured pruning of a single node at a time, coupled with SGD re-training
[31, 32] (Figure 4C) requires 70 training epochs to identify a recovery path.
In comparison, the geodesic strategy finds paths that quantitatively out-
perform the iterative prune-train procedure and obtains these paths with only
10 training epochs (figure 4C,D).

Additionally, the same geodesic strategy enables us to dynamically shift
networks between different weight configurations (eg from a dense to sparse or
vice-versa) while maintaining performance (Figure 4F). The rapid shifting of
networks is relevant for networks on neuromorphic hardware to ensure that the
real-time functionality of the hardware isn’t compromised while transitioning
between different power configurations.

##  6 Discussion

We have established a mathematical framework to analyze resilience of neural
networks through the lens of differential geometry. We introduce a functional
distance metric on a Riemmanian weight manifold and apply the metric tensor,
covariant derivative, and the geodesic to predict the response of networks to
local and global damage. Mathematically, our work forms new connections
between machine learning and differential geometry. Practically, we develop
new procedures for (i) identifying vulnerabilities in neural networks and (ii)
compensating for network damage in real-time through computationally efficient
weight updates, enabling their rapid recovery. As neural networks are
increasingly deployed on edge devices with increased susceptibility to damage,
we believe these methods could be useful in a variety of practical
applications.

##  7 Broader Impact

The field of AI has grown by leaps and bounds in the last few years. As a
result, AI is increasingly being built into many critical applications across
the society. Additionally, to cater to the rising need of AI systems for real-
time applications, AI systems have been transitioning from cloud-
implementation to edge devices and neuromorphic hardware. Some of the real-
time critical applications that have actively adopted AI systems include (1)
decision making in the health-care industry, (2) real-time image and sensor
analysis in self-driving cars, (3) incorporation into IoT sensors and devices
installed in most households and (4) robotic control systems.

The failure of AI in any of these applications could be catastrophic. For
instance, errors committed by AI systems while classifying radiology reports
in the health-care industry, or the faulty real-time analysis of stream of
images being processed by AI systems in self-driving cars could lead to human
casualties. Hence, it has become extremely important for us to understand how
neural network architectures (performing critical applications) react to
perturbations, that could arise from many sources. AI implemented on the cloud
are a victim of DRAM (dynamic random access memory) errors that can occur at
surprising rates, either due to malicious attack or induced by high energy
particles. Additionally, the growing implementation of AI networks on physical
hardware (for instance, neuromorphic, edge devices) has made the need for
discovering damage-resilient networks and rapidly recovery damaged networks a
necessity.

Our paper lays down the mathematical framework to study resilience and
robustness of neural networks to damage and proposes algorithms to rapidly
recover networks experiencing damage. Our research will be extremely important
for AI systems implemented across many applications, as damage of systems is
inevitable and needs to be protected against. Although resilience and
robustness of AI systems is very important, there aren’t very many principled
studies on the same. To reduce the gap in our knowledge on the resilience of
AI, we propose a principled framework to understand the vulnerabilites of AI
networks. One of the immediate applications of our contribution would be the
design of damage-resilient networks and rapid recovery algorithms implemented
on neuromorphic hardware. We believe this research will be foundational as
neural networks are becoming ubiquitous across many applications, ranging from
rovers sent to mars to radiology applications.

## References

  * [1] Lazaros S Iliadis, Vera Kurkova and Barbara Hammer  “Brain-inspired computing and machine learning”  Springer, 2020
  * [2] Kunihiko Fukushima  “Neocognitron: A self-organizing neural network model for a mechanism of pattern recognition unaffected by shift in position”  In _Biological cybernetics_ 36.4 Springer, 1980, pp. 193–202
  * [3] Maciej Wołczyk, Jacek Tabor, Marek Śmieja and Szymon Maszke  “Biologically-Inspired Spatial Neural Networks”  In _arXiv preprint arXiv:1910.02776_ , 2019
  * [4] Anthony M Zador  “A critique of pure learning and what artificial neural networks can learn from animal brains”  In _Nature communications_ 10.1 Nature Publishing Group, 2019, pp. 1–7
  * [5] Dan Hendrycks and Thomas Dietterich  “Benchmarking neural network robustness to common corruptions and perturbations”  In _arXiv preprint arXiv:1903.12261_ , 2019
  * [6] Walter G Gonzalez, Hanwen Zhang, Anna Harutyunyan and Carlos Lois  “Persistence of neuronal representations through time and damage in the hippocampus”  In _Science_ 365.6455 American Association for the Advancement of Science, 2019, pp. 821–825
  * [7] Naomie Castor and Farid El Massioui  “Resilience after a neurological pathology: What impact on the cognitive abilities of patients with brain damage?”  In _Neuropsychological rehabilitation_ Taylor & Francis, 2018, pp. 1–19
  * [8] Eve Marder and Jean-Marc Goaillard  “Variability, compensation and homeostasis in neuron and network function”  In _Nature Reviews Neuroscience_ 7.7 Nature Publishing Group, 2006, pp. 563–574
  * [9] Shyam Srinivasan and Charles F Stevens  “Robustness and fault tolerance make brains harder to study”  In _BMC biology_ 9.1 Springer, 2011, pp. 46
  * [10] Anja Richter, Bernd Krämer, Esther K Diekhof and Oliver Gruber  “Resilience to adversity is associated with increased activity and connectivity in the VTA and hippocampus”  In _NeuroImage: Clinical_ 23 Elsevier, 2019, pp. 101920
  * [11] Paola Arlotta and Benedikt Berninger  “Brains in metamorphosis: reprogramming cell identity within the central nervous system”  In _Current opinion in neurobiology_ 27 Elsevier, 2014, pp. 208–214
  * [12] Monte A Gates, Rosemary A Fricker-Gates and Jeffrey D Macklis  “Reconstruction of cortical circuitry”  In _Progress in brain research_ 127 Elsevier, 2000, pp. 115–156
  * [13] Timothy H Murphy and Dale Corbett  “Plasticity during stroke recovery: from synapse to behaviour”  In _Nature Reviews Neuroscience_ 10.12 Nature Publishing Group, 2009, pp. 861–872
  * [14] Tara Keck et al.  “Synaptic scaling and homeostatic plasticity in the mouse visual cortex in vivo”  In _Neuron_ 80.2 Elsevier, 2013, pp. 327–334
  * [15] Y Kate Hong, Clay O Lacefield, Chris C Rodgers and Randy M Bruno  “Sensation, movement and learning in the absence of barrel cortex”  In _Nature_ 561.7724 Nature Publishing Group, 2018, pp. 542–546
  * [16] Roger Lewin  “Is your brain really necessary?”  In _Science_ 210.4475 JSTOR, 1980, pp. 1232–1234
  * [17] William Beecher Scoville and Brenda Milner  “Loss of recent memory after bilateral hippocampal lesions”  In _Journal of neurology, neurosurgery, and psychiatry_ 20.1 BMJ Publishing Group, 1957, pp. 11
  * [18] Kurt Hammerschmidt, Gabriela Whelan, Gregor Eichele and Julia Fischer  “Mice lacking the cerebral cortex develop normal song: insights into the foundations of vocal learning”  In _Scientific reports_ 5 Nature Publishing Group, 2015, pp. 8808
  * [19] Charles D Gilbert and Wu Li  “Adult visual cortical plasticity”  In _Neuron_ 75.2 Elsevier, 2012, pp. 250–264
  * [20] Guanpeng Li et al.  “Understanding error propagation in deep learning neural network (DNN) accelerators and applications”  In _Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis_ , 2017, pp. 1–12
  * [21] Austin P Arechiga and Alan J Michaels  “The robustness of modern deep learning architectures against single event upset errors”  In _2018 IEEE High Performance extreme Computing Conference (HPEC)_ , 2018, pp. 1–6  IEEE
  * [22] Bianca Schroeder, Eduardo Pinheiro and Wolf-Dietrich Weber  “DRAM errors in the wild: a large-scale field study”  In _ACM SIGMETRICS Performance Evaluation Review_ 37.1 ACM New York, NY, USA, 2009, pp. 193–204
  * [23] Don Monroe  “Neuromorphic computing gets ready for the (really) big time”  ACM New York, NY, USA, 2014
  * [24] Pavel Mach and Zdenek Becvar  “Mobile edge computing: A survey on architecture and computation offloading”  In _IEEE Communications Surveys & Tutorials_ 19.3 IEEE, 2017, pp. 1628–1656
  * [25] Donald Firesmith  “System Resilience: What Exactly is it?”  In _System Resilience: What Exactly is it?_ , 2019  URL: https://insights.sei.cmu.edu/sei_blog/2019/11/system-resilience-what-exactly-is-it.html
  * [26] Ari S Morcos, David GT Barrett, Neil C Rabinowitz and Matthew Botvinick  “On the importance of single directions for generalization”  In _arXiv preprint arXiv:1803.06959_ , 2018
  * [27] Nicholas Cheney, Martin Schrimpf and Gabriel Kreiman  “On the robustness of convolutional neural networks to internal architecture and weight perturbations”  In _arXiv preprint arXiv:1703.08245_ , 2017
  * [28] Detlef H Mache, József Szabados and Marcel G Bruin  “Trends and Applications in Constructive Approximation”  Springer Science & Business Media, 2006
  * [29] Yann LeCun, Léon Bottou, Yoshua Bengio and Patrick Haffner  “Gradient-based learning applied to document recognition”  In _Proceedings of the IEEE_ 86.11 Ieee, 1998, pp. 2278–2324
  * [30] Karen Simonyan and Andrew Zisserman  “Very deep convolutional networks for large-scale image recognition”  In _arXiv preprint arXiv:1409.1556_ , 2014
  * [31] Song Han, Jeff Pool, John Tran and William Dally  “Learning both weights and connections for efficient neural network”  In _Advances in neural information processing systems_ , 2015, pp. 1135–1143
  * [32] Jonathan Frankle and Michael Carbin  “The lottery ticket hypothesis: Finding sparse, trainable neural networks”  In _arXiv preprint arXiv:1803.03635_ , 2018