# Theory of phase-adaptive parametric cooling

Alekhya Ghosh  Max Planck Institute for the Science of Light, Staudtstraße 2,
D-91058 Erlangen, Germany  Department of Physics, Friedrich-Alexander-
Universität Erlangen-Nürnberg, Staudtstraße 7, D-91058 Erlangen, Germany
Pardeep Kumar  Max Planck Institute for the Science of Light, Staudtstraße 2,
D-91058 Erlangen, Germany  Fidel Jimenez  Pontificia Universidad Católica del
Perú, Av. Universitaria 1801, San Miguel 15088, Peru  Vivishek Sudhir  LIGO
Laboratory, Massachusetts Institute of Technology, Cambridge, Massachusetts
02139, USA  Department of Mechanical Engineering, Massachusetts Institute of
Technology, Cambridge, MA 02139, USA  Claudiu Genes  Max Planck Institute for
the Science of Light, Staudtstraße 2, D-91058 Erlangen, Germany  Department of
Physics, Friedrich-Alexander-Universität Erlangen-Nürnberg, Staudtstraße 7,
D-91058 Erlangen, Germany

###### Abstract

We propose an adaptive phase technique for the parametric cooling of
mechanical resonances. This involves the detection of the mechanical
quadratures, followed by a sequence of periodic controllable adjustments of
the phase of a parametric modulation. The technique allows the preparation of
the quantum ground state with an exponential loss of thermal energy, similarly
to the case of cold-damping or cavity self-cooling. Analytical derivations are
presented for the cooling rate and final occupancies both in the classical and
quantum regimes.

###### pacs:

42.50.Ar, 42.50.Lc, 42.72.-g

Mechanical resonators are oftentimes used in displacement Anetsberger _et al._
(2009), force Fogliano _et al._ (2021), acceleration Krause _et al._ (2012) or
mass sensing Sansa _et al._ (2020) applications. Bringing them to their
quantum ground state is of fundamental interest for research studying the
classical to quantum physics transition Romero-Isart _et al._ (2011). A
particularly successful path in achieving control over motion has been taken
in optomechanics where photons are interfaced with phonons Aspelmeyer _et al._
(2014). Cooling close to the quantum ground state of an isolated mechanical
resonance has been achieved both via cavity self-cooling Delić _et al._ (2019,
2020); Windey _et al._ (2019); Gonzalez-Ballestero _et al._ (2019) and via
externally imposed feedback Li _et al._ (2011); Millen _et al._ (2020); Rossi
_et al._ (2018); Tebbenjohanns _et al._ (2021); Magrini _et al._ (2021);
Whittle and et al (2021). In the former case, external driving of a cavity
optomechanical system leads to the possibility of resolved sideband cooling
Schliesser _et al._ (2006); Kippenberg and Vahala (2007). In the latter case,
detection of the light scattered from an optomechanical system allows for the
implementation of a cold-damping mechanism Dania _et al._ (2021), where a
viscous force proportional to the momentum is provided. Extensions to the
simultaneous feedback cooling of many mechanical resonances are also possible
Sommer and Genes (2019); Sommer _et al._ (2020). Both these two paths have
been extensively utilized theoretically and experimentally and their
advantages and limitations are well understood Genes _et al._ (2008);
Aspelmeyer _et al._ (2014); Komori _et al._ (2022). More recently, a third
option, dubbed parametric cooling, has arisen, as an efficient option for
cooling optically levitated particles Gieseler _et al._ (2012); Gonzalez-
Ballestero _et al._ (2021), atoms in cavities Sames _et al._ (2018) or nano-
electromechanical resonators Villanueva _et al._ (2012). For a generic
oscillator at resonance frequency $\Omega$, the main ingredient of this
technique is the periodic modulation of the spring constant at $2\Omega$. In
the particular case of levitated nanoparticles, this is achieved by feedback
control of the trapping potential Gieseler _et al._ (2013). The cooling
mechanism is then similar to cold-damping but with a non-exponential loss rate
of energy Gieseler _et al._ (2012); Rodenburg _et al._ (2016).

Figure 1: Parametric cooling via phase adaptive feedback. Generic oscillator
with natural frequency $\Omega$. Detection of the mechanical quadratures $q$
and $p$ allows for the design of an adaptive phase feedback strategy, based on
the parametric modulation of the spring constant at frequency $2\Omega$ and
adjustable modulation phase $\phi$.

Here, we propose a variation of the parametric cooling technique, which is not
based on a controlled cold-damping feedback loop Tebbenjohanns _et al._ (2021)
but, instead, on the adjustment of the modulation phase: we dub this technique
phase adaptive parametric cooling. The technique solely requires the detection
of the momentum and position quadratures of the oscillator, from which an
optimal modulation phase is deduced and fed back into the system (as
illustrated in Fig. 1). Parametric driving of an oscillator of displacement
$x$ with natural frequency $\Omega$ and mass $m$, consists in supplementing
the bare restoring force $\mathcal{F}=-m\Omega^{2}x$ with an extra small
modulation $\delta\mathcal{F}=2m\Gamma\Omega\cos{(2\Omega t+\phi)}x$ (under
the condition $\Gamma\ll\Omega$). Alternatively, this can be understood as a
periodic and phase dependent variation of the spring contant $k(t,\phi)$. The
key point in our treatment is the observation that the phase $\phi$ is a
crucial tuning knob: when properly adjusted, it can steer the system to a
cooling regime characterized by an exponential loss of energy at rate
$\Gamma$. For a thermal environment at average occupancy $n_{\text{th}}\gg 1$,
we analyze the classical competition between parametric cooling and thermal
reheating to predict a final occupancy as low as $\gamma
n_{\text{th}}/\Gamma$, where $\gamma$ is the intrinsic mechanical damping. The
same result holds in the quantum regime and can be generalized as well to the
simultaneous cooling of a number of oscillators Brandão _et al._ (2021);
Rudolph _et al._ (2020); Chauhan _et al._ (2020) via a single parametric
modulation drive.
The emergence of damping behavior is analytically shown by noting that a
parametric driven oscillator satisfies a Mathieu-like equation, under
perturbative conditions. This analytical solution allows for the
identification of the optimal phase of the parametric modulation and the
subsequent derivation of an efficient phase adaptive feedback procedure.
Classical and quantum limits for cooling efficiency can be derived and
compared to simulation of a set of stochastic equations. Finally, the
technique can easily be extended to the simultaneous cooling of many
mechanical modes.

Model and equations \- We consider a mechanical resonance (along some
quantized direction $\hat{x}$ of effective mass $m$ subjected to the standard
restoring force $\mathcal{F}=-m\Omega^{2}\hat{x}$ and to a weak parametric
modulation $\delta\mathcal{F}=2m\Gamma\Omega\cos{(2\Omega t+\phi)}\hat{x}$;
here $\Omega$ is the natural oscillation frequency and a modulation amplitude
$\Gamma\ll\Omega$ and phase $\phi$ are assumed. With the definition of the
zero point motional amplitude $x_{\text{zpm}}=\sqrt{\hbar/(m\Omega)}$, one can
introduce the dimensionless position quadrature
$\hat{q}=\hat{x}/x_{\text{zpm}}$ and the corresponding momentum quadrature
$\hat{p}$ such that $[\hat{q},\hat{p}]=i$. The Heisenberg-Langevin equations
of motion then can be derived as

$\displaystyle\frac{d\hat{q}}{dt}$ $\displaystyle=\Omega\hat{p},$ (1a)
$\displaystyle\frac{d\hat{p}}{dt}$
$\displaystyle=-\gamma\hat{p}-\Omega\hat{q}+2\Gamma\cos(2\Omega
t+\phi)\hat{q}+\hat{\zeta}(t),$ (1b)

describing a parametrically driven, quantum mechanical oscillator subject to
thermal noise via the term $\hat{\zeta}(t)$. In the absence of external
forces, the mechanical mode is in equilibrium with a thermal bath at
temperature $T_{\text{th}}$ to which an average occupancy $n_{\text{th}}$
corresponds. The action of the bath onto the mechanical resonance is modeled
via quantum stochastic noise operators $\hat{\zeta}(t)$ of zero average and
with two-time correlations
$\braket{\hat{\zeta}(t)\hat{\zeta}(t^{\prime})}=1/(2\pi)\textstyle\int d\omega
S_{\text{th}}(\omega)e^{-i\omega(t-t^{\prime})}$. The thermal power spectrum
is given by
$S_{\text{th}}(\omega)=\gamma\omega/\Omega\left\\{\coth{[\hbar\omega/(2k_{B}T_{\text{th}})]}+1\right\\}$.
For large mechanical quality factors $\mathcal{Q}_{\text{m}}=\Omega/\gamma\gg
1$, the thermal damping rate can be defined in terms of sidebands:
$\gamma=[S_{\text{th}}(\Omega)-S_{\text{th}}(-\Omega)]/2$ while the average
thermal occupancy is
$n_{\text{th}}=[S_{\text{th}}(\Omega)+S_{\text{th}}(-\Omega)]/(2\gamma)$.
Notice that the equipartition theorem implies that, in a thermal state the
variance in the two quadrature is the same and derivable as
$\braket{\hat{q}^{2}}=\braket{\hat{p}^{2}}=n_{\text{th}}+1/2$.
In a first step, we take an average over the Langevin equations and construct
a second order differential equation for the expectation value
$q=\braket{\hat{q}}$. We then follow a perturbative method which requires that
the modulation parameter $b=\Gamma/\Omega\ll 1$ is small. Eliminating the
trivial exponential damping by the transformation $q=\bar{q}e^{-\gamma t/2}$
allows for the derivation of a Mathieu-like equation
$\ddot{\bar{q}}+\left[\Omega^{\prime 2}-2b\Omega^{2}\cos(2\Omega
t+\phi)\right]\bar{q}=0$ where
$\Omega^{\prime}=\sqrt{\Omega^{2}-\gamma^{2}/4}$ (as
$\mathcal{Q}_{\text{m}}\gg 1$, one can safely approximate this term in the
following with $\Omega$). Further simplifications are then obtained via the
transformation to a dimensionless time variable $\bar{t}=\Omega t+\phi/2$. The
dynamics can now be exactly mapped onto the standard Mathieu equation

$\frac{d^{2}q_{M}}{d\bar{t}^{2}}+\left[1-2b\cos(2\bar{t})\right]q_{M}=0.$ (2)

In the following we will use approximate solution of this equation, obtained
in the limit $b\ll 1$. Notice that the reverse transformation from $\bar{t}$
to the real time variable is then done via the identification
$q(t)=q_{\text{M}}(\Omega t+\phi/2)e^{-\gamma t/2}$.

Figure 2: Single shot cooling a) Time evolution of the occupancy for three
randomly selected initial conditions $(q_{0}^{(j)},p_{0}^{(j)})$ for
$j=1,2,3$, each with an optimized phase $\phi_{0}^{(j)}$. All trajectories
show exponential cooling up to the turning points $\tau_{j}$ where the damped
solution is comparable to the growing one. The parameters are
$\gamma/\Omega=10^{-6}$, $n_{\text{th}}=10^{4}$ and $b=0.05$. The time is
expressed in inverse units of $\Omega$. b) Time dynamics of occupancy for
varying $b$ (other parameters same as in (a)) shows the direct dependence of
the cooling rate on $b$ and that the validity of the analytically predicted
exponential cooling holds even for relatively large $b=0.25$.

Parametric self-cooling solution – The general procedure for analytically
solving a differential equation in Mathieu form as in Eq. (2) is well known
Gradshteyn and Ryzhik (2007) and famously used, for example, in the theory of
harmonic trapping of ions in Paul traps Leibfried _et al._ (2003) (albeit in a
very different regime, where the modulation is much faster than the natural
frequency). In the limit of small modulation $b\ll 1$, one can show that, to a
very good approximation (see Appendix) the (classical) solution can be written
in terms of negative and positive damping components

$\displaystyle q=\mathcal{A}_{-}e^{-\frac{(\gamma+\Gamma)t}{2}}\sin{(\Omega
t+\phi^{\prime})}+\mathcal{A}_{+}e^{-\frac{(\gamma-\Gamma)t}{2}}\cos{(\Omega
t+\phi^{\prime})}$ (3)

with $\phi^{\prime}=\phi/2-\pi/4$. The solution above shows that parametric
driving brings in an additional optical damping rate $\Gamma$ (we have assumed
$b>0)$ that adds to the intrinsic thermal damping $\gamma$. The negative
solution indicates increased damping, while the positive solution leads to a
heating instability as soon as $\mathcal{\Gamma}$ becomes larger than
$\gamma$. The coefficients $\mathcal{A}_{\pm}$ are derived from the initial
conditions $q_{0}$ and $p_{0}$ (as detailed in Appendix) and their ratio
strongly depends on the driving phase $\phi$. For a fixed modulation phase
$\phi$, the initial conditions then dictate whether the system will follow a
cooling or heating dynamics in the transient regime; in the long time limit,
independently of the initial conditions, the system will blow up as the
heating solution will always prevail.
Let us ask for the conditions for which $\mathcal{A}_{-}\gg\mathcal{A}_{+}$,
such that the damped solution dominates at early times. Under the
approximation that $b\ll 1$, a simplified analytical expression can be found
(see Appendix)

$\phi_{\text{opt}}^{(0)}=\frac{\pi}{2}+2\tan^{-1}\left[\frac{q_{0}}{p_{0}}\right],$
(4)

that sets the proper, optimized value for the modulation phase. This allows
one to formulate the strategy for the single shot cooling mechanism: detection
of scattered light allows to infer the initial conditions, which in turn are
used to optimize the modulation phase. The resulting cooling dynamics occurs
for any initial conditions as illustrated in Fig. 2. The loss of energy is
exponential at the analytically predicted rate $\Gamma$ and dominates up to a
time $\tau$ where (roughly)
$\tau\approx(1/\Gamma)\ln(\mathcal{A}_{-}/\mathcal{A}_{+})$ (where the
positive and negative solutions are comparable). This means that the occupancy
reached at time $\tau$ is roughly $\mathcal{A}_{-}/\mathcal{A}_{+}$ times
smaller than the initial one. Analytically one can deduce
$\mathcal{A}_{-}/\mathcal{A}_{+}\approx 2(q_{0}^{2}-p_{0}^{2})/(bq_{0}^{2})$
(see Appendix) which indicates on average a reduction in thermal energy by a
factor of roughly $1/b$.
This result already suggests the mechanism to achieve control over the cooling
dynamics at arbitrarily long times: at regular times $j\delta\tau$ (with
$j=1...\mathcal{N}$ and the repetition time interval $\delta\tau<\tau$),
before heating starts to dominate, one detects the instantaneous values of
$q(j\delta\tau)$ and $p(j\delta\tau)$ and updates the phases
$\phi(j\delta\tau)$ to the values indicated by Eq. 4. Following this procedure
for $\mathcal{N}$ steps suggests the loss of energy roughly by a factor
$1/b^{\mathcal{N}}$. For values of $b$ around $0.1$, in only $6$ feedback
steps one can reduce an initial $n_{\text{th}}=10^{4}$ average occupancy to
well below unity.

Figure 3: Feedback phase adaptive cooling. a) Time evolution of the occupancy
for a given set of initial conditions assuming initially fixed phase (red
line) as opposed to adapted phase with $\delta\tau<\tau$ (blue line, fast
feedback) and $\delta\tau>\tau$ (magenta line, slow delayed feedback). b) Two
trajectories with different initial conditions randomly picked from an thermal
distribution with $n_{\text{th}}=10^{4}$. In the presence of thermal noise,
the cooling dynamics is purely exponential at the theoretically predicted loss
rate $\Gamma$. Close to the ground state, the competition between Brownian
noise thermal heating and parametric cooling shows the reach of an equilibrium
final occupancy at $n_{\text{final}}\approx\gamma n_{\text{th}}/\Gamma$.
Fluctuations around this value are owed to the competition between stochastic
processes involving thermal reheating and continuous optical cooling action.
c) Phase-space illustration of an initial thermal state ($10^{3}$ points in
orange) and the corresponding final cold state (in blue). d) Zoom in of the
final occupancy distribution in phase space, showing that the system reaches a
thermal state. The parameters are chosen as $\gamma/\Omega=10^{-6}$ and
$b=0.05$. The time is expressed in inverse units of $\Omega$.

Feedback cooling versus reheating – The performance of any cooling technique
involves the competition between the externally induced damping dynamics and
the inherent reheating due to the thermal environment in which the oscillator
is embedded. We first analyze the performance of an adaptive phase feedback
procedure in the presence of classical thermal noise (under the Markovian
approximation) modelled by as a Wiener process $\braket{dW(t)}$ of zero
average and variance $\braket{dW(t)^{2}}=dt$ equal to the infinitesimal time
increment. To this end one writes a set of coupled difference equations

$\displaystyle dq$ $\displaystyle=\Omega pdt,$ (5a) $\displaystyle dp$
$\displaystyle=-\gamma pdt-\Omega qdt+\mathcal{F}_{\text{mod}}dt+\sqrt{2\gamma
n_{\text{th}}}dW(t),$ (5b)

where the Brownian noise is normalized such that, in the absence of the trap
modulation force, the system thermalizes at rate $\gamma$ to the temperature
$T_{\text{th}}$ associated with an average quanta occupancy $n_{\text{th}}$.
For numerical simulations, one can model $dW(t)=\sqrt{dt}\mathcal{N}(0,1)$,
where $\mathcal{N}(0,1)$ describes a normally distributed random variable of
unit variance. A semi-analytical solution can be instead found by turning the
difference equations into a set of recurrence equations. We discretize the
time interval $[0,t]$ into $n$ steps of duration $dt=t/n$. The equations above
can be rewritten as
$\mathbf{v}_{n}-\mathcal{M}_{n}\mathbf{v}_{n-1}=\sqrt{2\gamma
n_{\text{th}}}\mathbf{u}dW_{n}$ where $\mathbf{v}=(q,p)^{\top}$ and
$\mathbf{u}=(0,1)^{\top}$ and the evolution matrix is defined as

$\mathcal{M}_{n}=\begin{bmatrix}1&\Omega dt\\\ -\Omega
dt+2\Gamma\cos\left[2\Omega(n-1)dt+\phi\right]dt&1-\gamma dt\end{bmatrix}.$
(6)

A formal analytical solution can be then found easily and can be written in
terms of time ordered matrices
$\mathcal{T}_{nj}=\mathcal{M}_{n}\mathcal{M}_{n-1}...\mathcal{M}_{j}$ in the
following form

$\mathbf{v}_{n}=\mathcal{T}_{n1}\mathbf{v}_{0}+\sqrt{2\gamma
n_{\text{th}}}\sum_{j=1}^{n-1}\mathcal{T}_{nj}\mathbf{u}dW_{n-j}.$ (7)

The first part in the equation above describes the deterministic evolution
from the initial conditions while the last part is the long term behavior
dominated by thermal noise. Notice that for $n_{\text{th}}=0$, the evolution
matrix is time independent and the $\mathcal{T}_{nj}$ is simply equal to
$\mathcal{M}^{(n-j)}$ which allows for analytical solutions.
The phase adaptive feedback algorithm is then as chosen as follows: at the
initial time $t=0$, the trapped modulation phase is picked at
$\phi_{\text{opt}}^{(0)}$, from the detected quadratures, to insure an initial
damping period. Monitoring of the $q$ and $p$ quadratures at regular time
intervals $j\delta\tau$ (with $j=1,...\mathcal{N}$) is then followed by an
update of the modulation phases $\phi_{\text{opt}}^{(j)}$ fixed by Eq. (4)
with the replacement of the instantaneous quadratures $q_{j}=q(j\delta\tau)$
and $p_{j}=p(j\delta\tau)$. In Fig. 3a, the performance of the technique is
exemplified on three different trajectories with identical initial conditions
but different feedback times (at zero environment temperature). For an
optimized $\phi_{\text{opt}}^{(0)}$ the red line shows the heating in the
absence of feedback, while regular interval feedback show either exponential
loss of energy when $\delta\tau<\tau$ (blue line) or imperfect exponential
loss for $\delta\tau>\tau$ (magenta line).
For large thermal occupancies $n_{\text{th}}=10^{4}$, the two trajectories
(picked randomly from a thermal distribution) shown in Fig. 3b show that for
quick feeedback, exponential loss is achieved towards the same final
occupancy. From equilibrium considerations Genes _et al._ (2008), the final
occupancy can be deduced as the ratio of the reheating rate of the ground
state $\gamma n_{\text{th}}$ and the total damping rate $\gamma+\Gamma$, such
that $n_{\text{final}}=\gamma n_{\text{th}}/(\gamma+\Gamma)$.

Quantum ground state cooling – The classical model of thermalization does not
include the characteristics of the quantum harmonic oscillator, i.e. that the
commutations between $\hat{q}$ and $\hat{p}$ which automatically impose a
minimal variance of $1/2$ in both quadratures even for zero temperature. To
this end, we instead consider the quantum dynamics of the system described by
Eqs. (1) in differential form and with a quantum white noise input that
includes the zero-point energy. Let us assume dynamics around the steady
state, characterized by vanishing average momentum and position and by an
optimally adjusted modulation phase $\phi$ that gives rise to a constant
cooling rate $\Gamma$. We can perform a Fourier transform of Eqs. (1) to
derive

$\chi^{-1}(\omega)\hat{q}(\omega)=\Gamma\left[e^{i\phi}\hat{q}(-2\Omega+\omega)+e^{-i\phi}\hat{q}(2\Omega+\omega)\right]+\hat{\zeta}(\omega),$
(8)

Figure 4: Simultaneous cooling. a) Final occupancy (lower panel, red - initial
occupancy, black - isolated cooling, blue - simultaneous cooling) of eight
equidistant modes (upper panel). b) Same as in a) but considering the case
where two modes are degenerate in frequency. c) Time dynamics of the cooling
process for all modes showing fluctuations associated with the competition
between the imposed cooling and the inherent heating dynamics owed to the
thermal bath. d) Same as in c) but with the clear message that degenerate
modes are decoupled from the cooling dynamics. The parameters are chosen as
$\gamma/\Omega=10^{-6}$, $b=0.05$ and the initial occupancy is close to
$10^{5}$ for all modes. The time is expressed in inverse units of $\Omega$.

where the mechanical susceptibility is defined as
$\chi(\omega)=\Omega\left(\Omega^{2}-\omega^{2}-i\gamma\omega\right)^{-1}$ and
the last term describes a sum of noise stemming from both the thermal bath and
the effective parametric cooling bath (introduced phenomenologically to
account for the fluctuation-dissipation theorem). Similar to the procedure
illustrated in Ref. Genes _et al._ (2008), we compute the power spectrum of
the position quadrature and integrate it over all frequencies to obtain the
final occupancy as $n_{\text{th}}=\braket{\hat{q}^{2}}+1/2$ (as we have
assumed that thermal equipartition requires
$\braket{\hat{q}^{2}}=\braket{\hat{p}^{2}}$). The power spectrum contains
Lorentzians of linewidth $\gamma+\Gamma$. Under the assumption that
$\gamma+\Gamma\ll\Omega$, one can estimate the final occupancy by only
considering power spectra components around $\pm\Omega$. To this end we write
$\omega=\pm\Omega+\delta$ where $\delta$ where is a small variation
$\delta\ll\Omega$. After integration and under the assumption that
$\Gamma>\gamma$ one can show that (see Appendix)

$\displaystyle\braket{\hat{q}^{2}}=\frac{\gamma(\Gamma+\gamma\sin{\phi})}{\Gamma^{2}-\gamma^{2}}\left(n_{\text{th}}+\frac{1}{2}\right)+\frac{\Gamma(\Gamma+\gamma\sin{\phi})}{2(\Gamma^{2}-\gamma^{2})}.$
(9)

The first contribution describes the effect of the reduction of the thermal
bath effect onto the oscillator, while the second term is necessary to secure
a larger than $1/2$ variance stemming from quantum fluctuations. A more
concise expression can be obtained under the more stringent conditions where
$\Gamma\gg\gamma$ and final occupancy
$n_{\text{final}}=\braket{\hat{q}^{2}}-1/2$ is expressed in a very simple form
as

$\displaystyle
n_{\text{final}}=\frac{\gamma}{\Gamma}\left(n_{\text{th}}+\frac{1}{2}\right)$
(10)

The result is not surprising and agrees with the stochastic classical analysis
previously described: the final temperature is roughly reduced by a factor
$\Gamma/\gamma$. This reduction is bounded from above by the mechanical
quality factor $\mathcal{Q}_{\text{m}}$ (as $\Gamma<\Omega$ is the required
condition for the validity of the perturbative solution).

Simultaneous cooling of multiple modes – The phase adaptive mechanism can be
extended to simultaneously cool a number $n_{\text{res}}$ of adjacent
mechanical resonances Brandão _et al._ (2021); Rudolph _et al._ (2020);
Chauhan _et al._ (2020). Denoting their frequencies with $\Omega_{j}$ with
$j=1,...n_{\text{res}}$, we assume a generalized modulation force
$\sum_{j}2\Gamma_{j}\cos(2\Omega_{j}t+\phi_{j})q_{j}$. Assuming spectrally
resolved detection of all quadratures, the set of modulation phases can be
extracted and used for the periodic adjustment of the modulation force at
intervals $j\delta\tau$. The results are presented in Fig. 4 both for
equidistant, well separated $n_{\text{res}}=8$ resonances, as well as for the
particular case of two degenerate modes. The final occupancy in steady state
is calculated for individual, isolated cooling (ignoring the presence of
adjacent resonances) and compared with the performance of the simultaneous
cooling technique. The results are consistent with previous treatments of
cold-damping Sommer and Genes (2019), for example, showing that cooling is
efficient as long as the modes are frequency separated by more than the
effective cooling rate $\Gamma_{j}$.

Discussions and outlook – As opposed to previous theoretical treatments
Rodenburg _et al._ (2016) of parametric cooling Gieseler _et al._ (2012), the
technique introduced here does not need external engineering of a cold-damping
feedback force and only makes use of a minimally invasive detection loop in
the process. The predicted dynamics is purely exponential in energy loss and
is captured in a fully analytical model both at the classical and at the
quantum level. While parametric cooling has been experimentally introduced for
the motional control of optically levitated particles Gieseler _et al._
(2012); Gonzalez-Ballestero _et al._ (2021), atoms in cavities Sames _et al._
(2018) or nano-electromechanical resonators Villanueva _et al._ (2012), our
description here is quite general involving solely the modulation of the
spring constant of the resonator. Therefore, possible applications could
extend to optomechanical systems such as a membrane-in-the-middle
optomechanical systems, to cooling of ions in traps or to the refrigeration of
phonon modes in solid-state systems.

Acknowledgments – We acknowledge financial support from the Max Planck Society
and the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) –
Project-ID 429529648 – TRR 306 QuCoLiMa (“Quantum Cooperativity of Light and
Matter”).

## References

  * Anetsberger _et al._ (2009) G. Anetsberger, O. Arcizet, Q. P. Unterreithmeier, R. Rivière, A. Schliesser, E. M. Weig, J. P. Kotthaus, and T. J. Kippenberg, “Near-field cavity optomechanics with nanomechanical oscillators,” Nature Physics 5, 909 (2009).
  * Fogliano _et al._ (2021) F. Fogliano, B. Besga, A. Reigue, L. Mercier de Lépinay, P. Heringlake, C. Gouriou, E. Eyraud, W. Wernsdorfer, B. Pigeau, and O. Arcizet, “Ultrasensitive nano-optomechanical force sensor operated at dilution temperatures,” Nature Communications 12, 4124 (2021).
  * Krause _et al._ (2012) A. G. Krause, M. Winger, T. D. Blasius, Q. Lin, and O. Painter, “A high-resolution microchip optomechanical accelerometer,” Nature Photonics 6, 768 (2012).
  * Sansa _et al._ (2020) M. Sansa, M. Defoort, A. Brenac, M. Hermouet, L. Banniard, A. Fafin, M. Gely, C. Masselon, I. Favero, G. Jourdan, and S. Hentz, “Optomechanical mass spectrometry,” Nature Communications 11, 3781 (2020).
  * Romero-Isart _et al._ (2011) O. Romero-Isart, A. C. Pflanzer, F. Blaser, R. Kaltenbaek, N. Kiesel, M. Aspelmeyer, and J. I. Cirac, “Large quantum superpositions and interference of massive nanometer-sized objects,” Phys. Rev. Lett. 107, 020405 (2011).
  * Aspelmeyer _et al._ (2014) M. Aspelmeyer, T. J. Kippenberg, and F. Marquardt, “Cavity optomechanics,” Rev. Mod. Phys. 86, 1391 (2014).
  * Delić _et al._ (2019) U. Delić, M. Reisenbauer, D. Grass, N. Kiesel, V. Vuletić, and M. Aspelmeyer, “Cavity cooling of a levitated nanosphere by coherent scattering,” Phys. Rev. Lett. 122, 123602 (2019).
  * Delić _et al._ (2020) U. Delić, M. Reisenbauer, K. Dare, D. Grass, V. Vuletić, N. Kiesel, and M. Aspelmeyer, “Cooling of a levitated nanoparticle to the motional quantum ground state,” Science 367, 892 (2020).
  * Windey _et al._ (2019) D. Windey, C. Gonzalez-Ballestero, P. Maurer, L. Novotny, O. Romero-Isart, and R. Reimann, “Cavity-based 3D cooling of a levitated nanoparticle via coherent scattering,” Phys. Rev. Lett. 122, 123601 (2019).
  * Gonzalez-Ballestero _et al._ (2019) C. Gonzalez-Ballestero, P. Maurer, D. Windey, L. Novotny, R. Reimann, and O. Romero-Isart, “Theory for cavity cooling of levitated nanoparticles via coherent scattering: Master equation approach,” Phys. Rev. A 100, 013805 (2019).
  * Li _et al._ (2011) T. Li, S. Kheifets, and M. G. Raizen, “Millikelvin cooling of an optically trapped microsphere in vacuum,” Nature Physics 7, 527 (2011).
  * Millen _et al._ (2020) J. Millen, T. S. Monteiro, R. Pettit, and A. N. Vamivakas, “Optomechanics with levitated particles,” Reports on Progress in Physics 83, 026401 (2020).
  * Rossi _et al._ (2018) M. Rossi, D. Mason, J. Chen, Y. Tsaturyan, and A. Schliesser, “Measurement-based quantum control of mechanical motion,” Nature 563, 53 (2018).
  * Tebbenjohanns _et al._ (2021) F. Tebbenjohanns, M. L. Mattana, M. Rossi, M. Frimmer, and L. Novotny, “Quantum control of a nanoparticle optically levitated in cryogenic free space,” Nature 595, 378 (2021).
  * Magrini _et al._ (2021) L. Magrini, P. Rosenzweig, C. Bach, A. Deutschmann-Olek, S. G. Hofer, S. Hong, N. Kiesel, A. Kugi, and M. Aspelmeyer, “Real-time optimal quantum control of mechanical motion at room temperature,” Nature 595, 373 (2021).
  * Whittle and et al (2021) C. Whittle and et al, “Approaching the motional ground state of a 10-kg object,” Science 372, 1333 (2021).
  * Schliesser _et al._ (2006) A. Schliesser, P. Del’Haye, N. Nooshi, K. J. Vahala, and T. J. Kippenberg, “Radiation pressure cooling of a micromechanical oscillator using dynamical backaction,” Phys. Rev. Lett. 97, 243905 (2006).
  * Kippenberg and Vahala (2007) T. Kippenberg and K. Vahala, “Cavity opto-mechanics,” Opt. Express 15, 17172 (2007).
  * Dania _et al._ (2021) L. Dania, D. S. Bykov, M. Knoll, P. Mestres, and T. E. Northup, “Optical and electrical feedback cooling of a silica nanoparticle levitated in a paul trap,” Phys. Rev. Research 3, 013018 (2021).
  * Sommer and Genes (2019) C. Sommer and C. Genes, “Partial optomechanical refrigeration via multimode cold-damping feedback,” Phys. Rev. Lett. 123, 203605 (2019).
  * Sommer _et al._ (2020) C. Sommer, A. Ghosh, and C. Genes, “Multimode cold-damping optomechanics with delayed feedback,” Phys. Rev. Research 2, 033299 (2020).
  * Genes _et al._ (2008) C. Genes, D. Vitali, P. Tombesi, S. Gigan, and M. Aspelmeyer, “Ground-state cooling of a micromechanical oscillator: Comparing cold damping and cavity-assisted cooling schemes,” Phys. Rev. A 77, 033804 (2008).
  * Komori _et al._ (2022) K. Komori, D. Ďurovčíková, and V. Sudhir, “Quantum theory of feedback cooling of an anelastic macromechanical oscillator,” Phys. Rev. A 105, 043520 (2022).
  * Gieseler _et al._ (2012) J. Gieseler, B. Deutsch, R. Quidant, and L. Novotny, “Subkelvin parametric feedback cooling of a laser-trapped nanoparticle,” Phys. Rev. Lett. 109, 103603 (2012).
  * Gonzalez-Ballestero _et al._ (2021) C. Gonzalez-Ballestero, M. Aspelmeyer, L. Novotny, R. Quidant, and O. Romero-Isart, “Levitodynamics: Levitation and control of microscopic objects in vacuum,” Science 374, 3027 (2021).
  * Sames _et al._ (2018) C. Sames, C. Hamsen, H. Chibani, P. A. Altin, T. Wilk, and G. Rempe, “Continuous parametric feedback cooling of a single atom in an optical cavity,” Phys. Rev. A 97, 053404 (2018).
  * Villanueva _et al._ (2012) L. G. Villanueva, R. B. Karabalin, M. H. Matheny, E. Kenig, M. C. Cross, and M. L. Roukes, “A nanoscale parametric feedback oscillator,” Nano Lett. 11, 5054 (2012).
  * Gieseler _et al._ (2013) J. Gieseler, L. Novotny, and R. Quidant, “Thermal nonlinearities in a nanomechanical oscillator,” Nature Physics 9, 806 (2013).
  * Rodenburg _et al._ (2016) B. Rodenburg, L. P. Neukirch, A. N. Vamivakas, and M. Bhattacharya, “Quantum model of cooling and force sensing with an optically trapped nanoparticle,” Optica 3, 318 (2016).
  * Brandão _et al._ (2021) I. Brandão, D. Tandeitnik, and G. T, “Coherent scattering-mediated correlations between levitated nanospheres,” Quantum Science and Technology 6, 045013 (2021).
  * Rudolph _et al._ (2020) H. Rudolph, K. Hornberger, and B. A. Stickler, “Entangling levitated nanoparticles by coherent scattering,” Phys. Rev. A 101, 011804 (2020).
  * Chauhan _et al._ (2020) A. K. Chauhan, O. Černotík, and R. Filip, “Stationary gaussian entanglement between levitated nanoparticles,” New Journal of Physics 22, 123021 (2020).
  * Gradshteyn and Ryzhik (2007) I. S. Gradshteyn and I. M. Ryzhik, _Table of integrals, series, and products_ , seventh ed. (Elsevier/Academic Press, Amsterdam, 2007).
  * Leibfried _et al._ (2003) D. Leibfried, R. Blatt, C. Monroe, and D. Wineland, “Quantum dynamics of single trapped ions,” Rev. Mod. Phys. 75, 281 (2003).

##  Appendix A Appendix A: Classical dynamics of a parametrically drive
mechanical resonator

The standard solution to the Mathieu equation is written as an infinite sum

$q_{M}(\bar{t})=\mathcal{E}_{-}e^{i\beta\bar{t}}\sum_{n=-\infty}^{n=\infty}C_{2n}e^{i2n\bar{t}}+\mathcal{E}_{+}e^{-i\beta\bar{t}}\sum_{n=-\infty}^{n=\infty}C_{2n}e^{-i2n\bar{t}}.$
(11)

The coefficients of all harmonics (see Ref. Leibfried _et al._ (2003)) can be
found to satisfy the following recursive equation

$C_{2n+2}-D_{2n}C_{2n}+C_{2n-2}=0,$ (12)

where

$D_{2n}=\frac{1-(2n+\beta)^{2}}{b}$ (13)

In the limit of small $b$, the first three terms in the expansion above
suffice to properly describe the trajectory: we reduce the analysis to
$n=0,\pm 1$. For $\beta=-1+x$, where $x$ is a complex number with amplitude
much smaller than unity, we find that $D_{0}=2x/b$ and $D_{2}=-2x/b$ while all
other $D$s are very large. Fixing (without loss of generality) $C_{0}=1$ and
truncating all coefficients for $|n|>1$ we obtain $x=ib$ and $C_{2}=i$, as the
only non-vanishing coefficient. Finally, we can write the solution as

$\displaystyle
q(t)=2\mathcal{A}_{-}e^{-\frac{(\gamma+\Gamma)t}{2}}\cos{\left(\Omega
t+\phi^{\prime}\right)}+2\mathcal{A}_{+}e^{-\frac{(\gamma-\Gamma)t}{2}}\sin{\left(\Omega
t+\phi^{\prime}\right)},$ (14)

where $\mathcal{A}_{-}=2\mathcal{E}_{-}e^{i\pi/4}e^{-b\phi/4}$ and
$\mathcal{A}_{+}=2\mathcal{E}_{+}e^{i\pi/4}e^{b\phi/4}$ and the newly
introduced phase is $\phi^{\prime}=\phi/2+\pi/4$. The initial conditions ask
that $q_{M}(\bar{t}=\frac{\phi}{2})=q(t=0)=q_{0}$ and
$\dot{q}_{M}(\bar{t}=\frac{\phi}{2})=p(t=0)=p_{0}$, which leads to the
following set of equations

$\displaystyle\mathcal{A}_{-}\cos{\phi^{\prime}}+\mathcal{A}_{+}\sin{\phi^{\prime}}$
$\displaystyle=q_{0},$ (15a)
$\displaystyle-\mathcal{A}_{-}\sin{\phi^{\prime}}+\mathcal{A}_{+}\cos{\phi^{\prime}}-\frac{b}{2}\mathcal{A}_{-}\cos{\phi^{\prime}}+\frac{b}{2}\mathcal{A}_{+}\sin{\phi^{\prime}}$
$\displaystyle=p_{0}.$ (15b)

In a first step, we neglect $\mathcal{A}_{+}$ assuming it is much smaller than
the damped solution and obtain from above

$\phi_{\text{opt}}^{(0)}=\frac{\pi}{2}+2\tan^{-1}\left[\frac{1}{\frac{p_{0}}{q_{0}}+\frac{b}{2}}\right].$
(16)

With the observation that, for the chosen optimal phase, the signs of the
$\sin$ and $\cos$ terms are always the same, we can obtain
$|\mathcal{A}_{-}|=\sqrt{q_{0}^{2}+p_{0}^{2}}$ thus equal to the initial
variance of the thermal state.

##  Appendix B Appendix B: Classical stochastic evolution

The coupled difference equations

$\displaystyle dq$ $\displaystyle=\Omega pdt,$ (17a) $\displaystyle dp$
$\displaystyle=-\gamma pdt-\Omega qdt+\mathcal{F}_{\text{mod}}dt+\sqrt{2\gamma
n_{\text{th}}}dW(t),$ (17b)

can be directly numerically simulated by modelling
$dW(t)=\sqrt{dt}\mathcal{N}(0,1)$, where $\mathcal{N}(0,1)$ describes a
normally distributed random variable of unit variance. However, more
analytical insight can be obtained by turning the difference equations into a
set of recurrence equations. We discretize the time interval $[0,t]$ into $n$
steps of duration $dt=t/n$ and can then rewrite the equations above as

$\mathbf{v}_{n}-\mathcal{M}_{n}\mathbf{v}_{n-1}=\sqrt{2\gamma
n_{\text{th}}}\mathbf{u}dW_{n}$ (18)

where $\mathbf{v}=(q,p)^{\top}$ and $\mathbf{u}=(0,1)^{\top}$ and the
evolution matrix is defined as

$\mathcal{M}_{n}=\begin{bmatrix}1&\Omega dt\\\ -\Omega
dt+2\Gamma\cos\left[2\Omega(n-1)dt+\phi\right]dt&1-\gamma dt\end{bmatrix}.$
(19)

A solution can be then found easily and can be written in terms of time
ordered matrices
$\mathcal{T}_{nj}=\mathcal{M}_{n}\mathcal{M}_{n-1}...\mathcal{M}_{j}$ in the
following form

$\mathbf{v}_{n}=\mathcal{T}_{n1}\mathbf{v}_{0}+\sqrt{2\gamma
n_{\text{th}}}\sum_{j=1}^{n-1}\mathcal{T}_{nj}\mathbf{u}dW_{n-j}.$ (20)

As a simple check, let us describe solely the thermalization dynamics of an
unmodulated oscillator (setting $\Gamma=0$). The time ordered matrices are
much simpler now: $\mathcal{T}_{n1}=\mathcal{M}^{n}$ and
$\mathcal{T}_{nj}=\mathcal{M}^{n-j}$. Under the assumption that
$\gamma\ll\Omega$, diagonalization of the matrix $\mathcal{M}=S\Lambda S^{-1}$
is straightforward in terms of the two eigenvalues of $\mathcal{M}$ equal to
$\lambda_{1}=1-\gamma dt/2-i\Omega dt,\,\lambda_{2}=1-\gamma dt/2+i\Omega dt$.
Notice that the two eigenvalues can be rewritten as $\lambda_{1}=(1-\gamma
dt/2)e^{-i\Omega dt}=re^{i\theta}$ and $\lambda_{2}=re^{-i\theta}$ where
$r=1-\gamma t/2n$ and $\theta=\Omega t/n$. The resulting quadratures after $n$
steps are written as

$\displaystyle q_{n}$
$\displaystyle=\frac{\lambda_{1}^{n}+\lambda_{2}^{n}}{2}q_{0}+i\,\frac{\lambda_{1}^{n}-\lambda_{2}^{n}}{2}p_{0}+\sqrt{2\gamma
n_{\text{th}}}\,\sum_{j=0}^{n-1}\frac{\lambda_{1}^{j}-\lambda_{2}^{j}}{2}dW_{j},$
(21a) $\displaystyle p_{n}$
$\displaystyle=\frac{\lambda_{1}^{n}+\lambda_{2}^{n}}{2}p_{0}-i\,\frac{\lambda_{1}^{n}-\lambda_{2}^{n}}{2}q_{0}+\sqrt{2\gamma
n_{\text{th}}}\,\sum_{j=0}^{n-1}\frac{\lambda_{1}^{j}+\lambda_{2}^{j}}{2}dW_{j},$
(21b)

where the deterministic parts describe simply the oscillatory weakly damped
transient evolution and the last terms are the effect of the thermal
environment. In the large $n$ limit we find a closed expression

$\displaystyle q(t)$ $\displaystyle=e^{-\gamma t/2}\left[q_{o}\cos(\Omega
t)+p_{0}\sin(\Omega t)\right]+\lim_{n\to\infty}\sqrt{2\gamma
n_{\text{th}}}\,\sum_{j=0}^{n-1}r^{j}\sin(j\theta)dW_{j},$ (22a)
$\displaystyle p(t)$ $\displaystyle=e^{-\gamma t/2}\left[p_{0}\cos(\Omega
t)-q_{0}\sin(\Omega t)\right]+\lim_{n\to\infty}\sqrt{2\gamma
n_{\text{th}}}\,\sum_{j=0}^{n-1}r^{j}\cos(j\theta)dW_{j},$ (22b)

where we have used that $\lim_{n\to\infty}(1-\gamma t/2n)^{n}=e^{-\gamma
t/2}$. From these expressions, one can estimate that in steady state
$\braket{q^{2}}_{\text{ss}}=\braket{p^{2}}_{\text{ss}}=n_{\text{th}}$ by using
$\braket{dW_{j}dW_{j^{\prime}}}=\delta_{jj^{\prime}}t/n$ and evaluating the
limit $\lim_{n\to\infty}\sum_{j=0}^{n-1}r^{2j}\sin^{2}(j\theta)=1/(2\gamma
t)$.

##  Appendix C Appendix C: Quantum ground state cooling

Let us assume steady state and fix the modulation phase to $\phi$ such that
the cooling rate is exponential at rate $\Gamma$ and both the position and
momentum expectation values vanish. We will compute the final occupancy from
the variance in position (assuming thermal equipartition between momentum and
position quadratures) obtained from the power spectrum in the Fourier domain
$S_{q}(\omega)$ via an integration
$\braket{\hat{q}^{2}}=1/(2\pi)\int_{-\infty}^{\infty}d\omega S_{q}(\omega)$.
The Fourier transform of Eqs. (1) leads to

$\displaystyle-i\omega\hat{q}(\omega)$ $\displaystyle=\Omega\hat{p}(\omega),$
(23a) $\displaystyle-i\omega\hat{p}(\omega)$
$\displaystyle=-\gamma\hat{p}(\omega)-\Omega\hat{q}(\omega)+\Gamma\left[e^{i\phi}\hat{q}(\omega-2\Omega)+e^{-i\phi}\hat{q}(\omega+2\Omega)\right]+\hat{\zeta}(\omega).$
(23b)

Notice that, for $\Gamma=0$, the above equations lead to
$\hat{q}(\omega)=\chi(\omega)\hat{\zeta}(\omega)$ where the mechanical
susceptibility is

$\chi(\omega)=\frac{\Omega}{(\Omega^{2}-\omega^{2})-i\gamma\omega}.$ (24)

This indicates a simple way of computing the variance by using the
correlations of the thermal environment
$\braket{\zeta(\omega)\zeta(\omega^{\prime})}=S_{\text{th}}(\omega)\delta(\omega+\omega^{\prime})$
where $S_{\text{th}}(\omega)$ is the thermal power spectrum with sidebands
$S_{\text{th}}(-\Omega)=2\gamma n_{\text{th}}$ and
$S_{\text{th}}(\Omega)=2\gamma(n_{\text{th}}+1)$. As for very high mechanical
quality factors $\mathcal{Q}_{\text{m}}\gg 1$, the susceptibility is a very
sharply peaked function around $\pm\Omega$ one can approximate it around the
two poles at $\pm\Omega$ by expanding it in terms of a small quantity
$\gamma\ll|\Delta|\ll\Omega$. The expansion is then
$\chi(-\Omega+\Delta)\simeq 1/(2\Delta+i\gamma)$ and
$\chi(+\Omega+\Delta)\simeq 1/(-2\Delta-i\gamma)$. The variance can then
approximated by
$\braket{\hat{q}^{2}}(t)=1/(2\pi)\int_{-\infty}^{\infty}d\Delta(4\Delta^{2}+\gamma^{2})^{-1}(S_{\text{th}}(-\Omega)+S_{\text{th}}(\Omega))=n_{\text{th}}+1/2$
as expected. We made use of the integral
$1/(2\pi)\int_{-\infty}^{\infty}d\Delta
1/(4\Delta^{2}+\gamma^{2})=\pi/(2\gamma)$.
We can now rewrite Eqs. (23) as a recursive equation

$\hat{q}(\omega)=\Gamma\chi(\omega)\left[e^{i\phi}\hat{q}(-2\Omega+\omega)+e^{-i\phi}\hat{q}(2\Omega+\omega)\right]+\chi(\omega)\hat{\zeta}(\omega).$
(25)

We the proceed as above making small variations around $\pm\Omega$ and
assuming that the only components which contribute to the power spectrum of
the position quadrature are $\hat{q}(\pm\Omega+\Delta)$. We can now separate
two coupled equations

$\displaystyle\hat{q}(-\Omega+\Delta)$
$\displaystyle=\Gamma\chi(-\Omega+\Delta)e^{-i\phi}\hat{q}(\Omega+\Delta)+\chi(-\Omega+\Delta)\hat{\zeta}(-\Omega+\Delta),$
(26a) $\displaystyle\hat{q}(\Omega+\Delta)$
$\displaystyle=\Gamma\chi(\Omega+\Delta)e^{i\phi}\hat{q}(-\Omega+\Delta)+\chi(\Omega+\Delta)\hat{\zeta}(\Omega+\Delta).$
(26b)

and invert them to find the solutions for $q(-\Omega+\Delta)$ and
$q(\Omega+\Delta)$ expressed as

$\displaystyle\hat{q}(-\Omega+\Delta)$
$\displaystyle=-\bar{\chi}(\Delta)\left[-(2\Delta+i\gamma)\zeta(-\Omega+\Delta)+\Gamma
e^{-i\phi}\zeta(\Omega+\Delta)\right],$ (27a)
$\displaystyle\hat{q}(\Omega+\Delta)$
$\displaystyle=-\bar{\chi}(\Delta)\left[\Gamma
e^{i\phi}\zeta(-\Omega+\Delta)+(2\Delta+i\gamma)\zeta(\Omega+\Delta)\right].$
(27b)

The modified mechanical susceptibility is approximated by the following
expression (under the assumption that $|\Delta|\ll\Omega$)

$\bar{\chi}(\Delta)=\frac{1}{4\Delta^{2}+\Gamma^{2}-\gamma^{2}+4i\gamma\Delta}.$
(28)

The denominator shows the presence of the optical damping rate $\Gamma$ but
not as a broadening adding to $\gamma$ as is the case for cold damping or
cavity cooling. We now again use the correlations of the noise term in the
frequency domain however supplemented with the bath responsible with damping
at rate $\Gamma$ and with zero temperature. Notice that, for small $\Delta$,
the only contributing terms are
$\braket{\zeta(-\Omega+\Delta)\zeta(\Omega+\Delta^{\prime})}\simeq
S(-\Omega)\delta(\Delta+\Delta^{\prime})$ and
$\braket{\zeta(\Omega+\Delta)\zeta(-\Omega+\Delta^{\prime})}\simeq
S(\Omega)\delta(\Delta+\Delta^{\prime})$ where now $S(-\Omega)=2\gamma
n_{\text{th}}$ and $S(\Omega)=2\gamma(n_{\text{th}}+1)+2\Gamma$. This leads to
the following contributions

$\displaystyle\braket{\hat{q}(-\Omega+\Delta)\hat{q}(-\Omega+\Delta^{{}^{\prime}})}$
$\displaystyle=$ $\displaystyle|\bar{\chi}(\Delta)|^{2}\Gamma
e^{-i\phi}\left[2\Delta\left(S(\Omega)-S(-\Omega)\right)-i\gamma\left(S(-\Omega)+S(\Omega)\right)\right]\delta(\Delta+\Delta^{{}^{\prime}}),$
(29a)
$\displaystyle\braket{\hat{q}(\Omega+\Delta)\hat{q}(\Omega+\Delta^{{}^{\prime}})}$
$\displaystyle=$ $\displaystyle|\bar{\chi}(\Delta)|^{2}\Gamma
e^{i\phi}\left[2\Delta\left(S(\Omega)-S(-\Omega)\right)+i\gamma\left(S(-\Omega)+S(\Omega)\right)\right]\delta(\Delta+\Delta^{{}^{\prime}}),$
(29b)
$\displaystyle\braket{\hat{q}(-\Omega+\Delta)\hat{q}(\Omega+\Delta^{{}^{\prime}})}$
$\displaystyle=$
$\displaystyle|\bar{\chi}(\Delta)|^{2}\left[(4\Delta^{2}+\gamma^{2})S(-\Omega)+\Gamma^{2}S(\Omega))\right]\delta(\Delta+\Delta^{{}^{\prime}}),$
(29c)
$\displaystyle\braket{\hat{q}(\Omega+\Delta)\hat{q}(-\Omega+\Delta^{{}^{\prime}})}$
$\displaystyle=$
$\displaystyle|\bar{\chi}(\Delta)|^{2}\left[\Gamma^{2}S(-\Omega)+(4\Delta^{2}+\gamma^{2})S(\Omega))\right]\delta(\Delta+\Delta^{{}^{\prime}}).$
(29d)

We can now add all the contributions to find the position power spectrum via
the following integral

$\braket{\hat{q}^{2}}=\frac{1}{2\pi}\int_{-\infty}^{\infty}d\Delta|\bar{\chi}(\Delta)|^{2}\left\\{4\gamma(2\gamma\Gamma\sin{\phi}+4\Delta^{2}+\Gamma^{2}+\gamma^{2})(n_{\text{th}}+\frac{1}{2})+2\Gamma\left(4\Delta(\gamma+\Gamma)\cos{\phi}+2\gamma\Gamma\sin{\phi}+4\Delta^{2}+\Gamma^{2}+\gamma^{2}\right)\right\\}$
(30)

The result for the relevant case $\Gamma>\gamma$ is

$\displaystyle\braket{\hat{q}^{2}}=\frac{\gamma(\Gamma+\gamma\sin{\phi})}{\Gamma^{2}-\gamma^{2}}\left(n_{\text{th}}+\frac{1}{2}\right)+\frac{\Gamma(\Gamma+\gamma\sin{\phi})}{2(\Gamma^{2}-\gamma^{2})}.$
(31)

Notice that in the limit where the additional damping dominates
$\Gamma\gg\gamma$ one can simplify the expression above and compute the final
occupancy $n_{\text{final}}=\braket{\hat{q}^{2}}-1/2$ to lead to

$\displaystyle
n_{\text{final}}=\frac{\gamma}{\Gamma}(n_{\text{th}}+\frac{1}{2})$ (32)

In the opposite case where $\Gamma<\gamma$ the variance reads

$\displaystyle\braket{\hat{q}^{2}}=\frac{\gamma(\gamma+\Gamma\sin{\phi})}{\gamma^{2}-\Gamma^{2}}\left(n_{\text{th}}+\frac{1}{2}\right)+\frac{\Gamma(\gamma+\Gamma\sin{\phi})}{2(\gamma^{2}-\Gamma^{2})}.$
(33)

For $\Gamma=0$ we recover the expected result
$\braket{\hat{q}^{2}}=n_{\text{th}}+1/2$.


that $\pi$ is preserved from $s,h$ to $s',h'$. The cases for $\cmdlock{}/\cmdunlock{}$ are trivial.

We call the assumption that $\forall \e\in \alpha(\set{\tr})_W. \e\not\properprefix\pi$, the stability assumption.

Case: $c=(x:=y)$ or $c=(x:=\pi')$ or $c=(x := \cmdalloc)$.

Clearly, $\lval{\pi}{s',h'}$ is defined.

By the stability assumption and
Definition <ref> ($\exec$), it follows that $x\neq \rootvar{\pi}$.
Therefore, by the concrete semantics of $c$, $s(\rootvar{\pi})=s'(\rootvar{\pi})$,
and in all cases, $h\subseteq h'$.
Thus, by Lemma <ref>, it follows that
$\pathimg{h'}{s'}{\pi} = \pathimg{h}{s}{\pi}$  and

The acyclicity of $\pi$ in $s',h'$ follows from Lemma <ref> and the fact
that $\pathimg{h'}{s'}{\pi} = \pathimg{h}{s}{\pi}$.

We show that $\pi$ is disconnected in $s',h'$.  Given
$\pathimg{h}{s}{\pi} = \pathimg{h'}{s'}{\pi}$, the second clause
of the disconnectedness definition carries over directly from
the disconnectedness of $\pi$ in $s,h$ whenever $c$ is
$x :=y $ or $x:=\pi'$ (because $h=h'$). The case where $c$ is $x:=\cmdalloc$
follows by the same reasoning, plus the fact that all newly-allocated
addresses point to the newly allocated location, and therefore by construction, cannot alias
addresses in $h$.

Given the disconnectedness of $\pi$ in $s,h$ and the fact that for all variables
$w\neq x$, $s(w)=s'(w)$, we only need to show
$s'(x)\notin\locn{\pathimg{h'}{s'}{\pi}}$ for the first clause.

If $c=(x:=y)$ then $s'(x)=s(y)$.
Given that $y\neq \rootvar{\pi}$, due to the stability assumption,
and the first clause of the definition of disconnectedness for $\pi$ in $s,h$, we have that
(since $\pathimg{h}{s}{\pi} = \pathimg{h'}{s'}{\pi}$), as required.

If $c=(x:=\pi')$ then $s'(x)=h(\lval{\pi'}{s,h})$.
Given that $h=h'$ and that $\pathimg{h}{s}{\pi} = \pathimg{h'}{s'}{\pi}$, it
suffices to show $h(\lval{\pi'}{s,h})\notin\locn{\pathimg{h}{s}{\pi}}$.
By the stability assumption,  $\pi'\not\prefix\pi$, thus applying Lemma <ref>,
we obtain $\lval{\pi'}{s,h}\notin\dom{\pathimg{h}{s}{\pi}}$.
By the contrapositive of the second clause of disconnectedness of $\pi$ in $s,h$ follows that
$h(\lval{\pi'}{s,h})\notin\dom{\pathimg{h}{s}{\pi}}$, as required.

If $c=(x:=\cmdalloc)$ then $s'(x)=\ell$ for some $\ell\notin\dom{h}$.
Clearly, $\ell\notin\locn{\pathimg{h}{s}{\pi}}$ and since
$\pathimg{h}{s}{\pi} = \pathimg{h'}{s'}{\pi}$, it follows that
$s'(x)=\ell\notin\locn{\pathimg{h'}{s'}{\pi}}$, as required.

Case: $c=(\pi':=x)$, where $\pi'\neq\pi$.

By the preservation of well-behavedness $c$ is safe to execute, thus, $\lval{\pi'}{s,h}\in\dom{h}$.
Also, $s=s'$.

The stability assumption forces $\pi'\not\prefix\pi$, thus by Lemma <ref> we have that
This, in addition with the semantics of the assignment to $\pi'$, means that
for all addresses $\alpha\in\dom{\pathimg{h}{s}{\pi}}$, $h(\alpha)=h'(\alpha)$.
Thus Lemma <ref> applies, giving
$\pathimg{h}{s}{\pi}=\pathimg{h'}{s'}{\pi}$ and
$\dom{\pathimg{h}{s}{\pi}} = \dom{\pathimg{h'}{s'}{\pi}}$.

For disconnectedness, first, $\lval{\pi}{s',h'}$ is clearly defined;
the first clause holds due to the fact that $s=s'$, $\pathimg{h}{s}{\pi}=\pathimg{h'}{s'}{\pi}$,
and the assumption of disconnectedness of $\pi$ in $s,h$.

By the semantics of the command, $\locn{h}=\locn{h'}$.
Therefore, the second clause of disconnectedness in $s',h'$ holds for all
addresses $\alpha\neq\lval{\pi'}{s,h}$,
by the second clause of the disconnectedness in $s,h$.
By the first clause and stability ($x\neq\rootvar{\pi}$) follows that
covering the case where $\alpha=\lval{\pi'}{s,h}$.

As before, acyclicity follows from the fact that $\pathimg{h}{s}{\pi}=\pathimg{h'}{s'}{\pi}$,
the acyclicity in $s,h$ and Lemma <ref>.

Case: $c=(\pi:=x)$.

By the acyclicity of $\pi$ in $s,h$ and the semantics of the store it follows that
Thus, Lemma <ref> applies and we have
$\lval{\pi}{s,h}=\lval{\pi}{s',h'}$ and
$\dom{\pathimg{h}{s}{\pi}} = \dom{\pathimg{h'}{s'}{\pi}}$.

Disconnectedness: The first clause follows from the disconnectedness in $s,h$, $s=s'$
and $\locn{\pathimg{h}{s}{\pi}} = \locn{\pathimg{h'}{s'}{\pi}}$
(because $\dom{\pathimg{h}{s}{\pi}} = \dom{\pathimg{h'}{s'}{\pi}}$).

For the second clause, we need to show that for all $\alpha\in\dom{h'}$,
if $h'(\alpha)\in\locn{\pathimg{h'}{s'}{\pi}}$, then $\alpha\in\dom{\pathimg{h'}{s'}{\pi}}$.
As $\dom{\pathimg{h}{s}{\pi}} = \dom{\pathimg{h'}{s'}{\pi}}$ and $\dom{h'}=\dom{h}$,
this is equivalent to showing that for all $\alpha\in\dom{h}$,
if $h'(\alpha)\in\locn{\pathimg{h}{s}{\pi}}$, then $\alpha\in\dom{\pathimg{h}{s}{\pi}}$.

If $\alpha\neq\lval{\pi}{s,h}=\lval{\pi}{s',h'}$, then $h(\alpha)=h'(\alpha)$. Thus,
if $h'(\alpha)\in\locn{\pathimg{h}{s}{\pi}}$
then $h(\alpha)\in\locn{\pathimg{h}{s}{\pi}}$
thus $\alpha\in\dom{\pathimg{h}{s}{\pi}}$ (by disconnectedness in $s,h$).

If $\alpha=\lval{\pi}{s,h}=\lval{\pi}{s',h'}$, then $h'(\alpha)=s(x)$.
By stability, $x\neq\rootvar{\pi}$, thus by the first clause of disconnectedness,
$s(x)\notin\locn{\pathimg{h}{s}{\pi}}$, as required.

Acyclicity: we wish to prove that for every $\pi''\prefix\pi'\prefix\pi$,
Using Lemma <ref>, this reduces to showing
(because for any path $\pi'\prefix\pi$, $\lval{\pi'}{s,h}=\lval{\pi'}{s',h'}$).
We distinguish two cases.

Let $\lval{\pi''}{s,h}\neq\lval{\pi}{s,h}$.
Then $h'(\lval{\pi''}{s,h}) = h(\lval{\pi''}{s,h})\neq\locn{\lval{\pi''}{s,h}}$,
by the acyclicity in of $\pi$ in $s,h$.

Let $\lval{\pi''}{s',h'}=\lval{\pi}{s',h'}$.
This means $h'(\lval{\pi''}{s',h'})= s(x) \neq \locn{\lval{\pi'}{s,h}}$.
by stability and the disconnectedness in $s,h$.

Let $h=h_0 \cup h_F$ (where $\cup$ is disjoint union), and let $\tr\in\sem{C}\angled{s,h,L}$ such that $\llast{\tr}=\angled{s',h',L'}$.
If $\sem{C}\angled{s,h_F,L}\neq\emptyset$ then there exists $\tr'\in\sem{C}\angled{s,h_F,L}$
such that $\llast{\tr'}=\angled{s',h'_F,L'}$, and $h'=h_0\cup h'_F$.

By the semantics of the language and standard separation-logic style reasoning.

Suppose that $\mbody(\m)$ has no procedure calls, that $\pi$ is an access path which is
disconnected and acyclic in the well-behaved state $s,h$,
that $\tr\in\sem{\cmdskip;\m(\bar{\e})}\angled{s,h,L}$ is a non-empty trace,
that $\llast{\tr}=\angled{\cmdpop{},s,h',\_}$, that $\alpha(\set{\tr})=\angled{W,\_,\_}$, and that
$\forall\e\in W.\, \e\not\properprefix\pi$.

Then, $\pi$ is preserved from $s,h$ to $s,h'$.

Let $\e=(\e_1,\ldots,\e_n)$ and let $\pi=z.\fs$.
Since $\tr$ is non-empty, for all $i\le n$, $\eval{\e_i}{s,h}=\ell_i$ (and is defined).
By the semantics for method calls and those for simple commands ($\cmdpop$) we have
$\tr = \tr_{\text{push}} \append \tr' \append \tr_{\text{pop}}$ where

* $\tr_{\text{push}} = [\angled{\cmdpush{\e_1,\ldots,\e_n}, \hat s \cons s, h, L}]$,

* $\tr' \in \sem{\body{\m}}\angled{\hat s \cons s, h, L}$, where $\llast{\tr'}=\angled{\_,\hat s' \cons s,h',L'}$, and

* $\tr_{\text{pop}}=[\angled{\cmdpop{}, s, h', L'}]$.

$\hat s \defeq s_{\nil}[\argm_1 \mapsto \ell_1]\cdots[\argm_n\mapsto\ell_n]$.

By the assumption that $\forall\e\in W.\,\e\not\properprefix\pi$ (we call this the stability assumption)
and the definition of
how $\exec$ modifies the $W$ component for
the method call case, we can assume there is at most one $\e_i$ such that
$\e_i= z$, or $\e_i=\pi_i\properprefix\pi$.

Case: there is no $i\le n$ such that $\e_i= z$, or $\e_i=\pi_i\properprefix\pi$.
Fix the stacks $\tilde s \defeq \hat s[\argm_{n+1}\mapsto s(z)]$
and $\tilde s' \defeq \hat s'[\argm_{n+1}\mapsto s(z)]$.
Recall that, by convention, $\body{\m}$ does not use $\argm_{n+1}$;
thus, since $\tr' \in \sem{\body{\m}}\angled{\hat s \cons s, h, L}$ and where $\llast{\tr'}=\angled{\_,\hat s' \cons s,h',L'}$,
it will also be the case that there is another “isomorphic” trace $\tr''$ such that
$\tr'' \in \sem{\body{\m}}\angled{\tilde s \cons s, h, L}$, where $\llast{\tr''}=\angled{\_,\tilde s' \cons s,h',L'}$.

We now show that $\rho=\argm_{n+1}.\fs$ is disconnected and acyclic in
$\tilde s, h$. Note that by Lemma <ref>,
$\pathimg{h}{s}{\pi} = \pathimg{h}{\tilde s}{\rho}$ and $\lval{\pi}{s,h}=\lval{\rho}{\tilde s, h}$,
therefore also $\dom{\pathimg{h}{s}{\pi}} = \dom{\pathimg{h}{\tilde s}{\rho}}$.

* First clause of disconnectedness: We need to show that for all variables
$x\neq \argm_{n+1}$, $\tilde s(x)\notin\locn{\pathimg{h}{\tilde s}{\rho}}$.
By the definition of $\tilde s$, this means showing that for all $i\le n$,
$\tilde s(\argm_i) = \hat s(\argm_i)\notin\locn{\pathimg{h}{\tilde s}{\rho}}$.
But $\hat s(\argm_i) = \eval{\e_i}{s,h}$.
We have two cases:

* If $\e_i= x$ for
some variable $x$, then $\eval{\e_i}{s,h} = s(x)$. By the stability assumption,
we have that $x\neq z$ and by the disconnectedness of $\pi$ in $s,h$
we have that $s(x)\notin\locn{\pathimg{h}{s}{\pi}}$. This completes the proof
through the fact that $\pathimg{h}{s}{\pi} = \pathimg{h}{\tilde s}{\rho}$.

* if $\e_i=\pi'$ then $\eval{\e_i}{s,h}=h(\lval{\pi'}{s,h})$.
By the stability assumption, we have that $\pi'\not\properprefix\pi$. If $\pi'\neq\pi$
then by Lemma <ref>, we have that
$\lval{\pi'}{s,h}\notin\dom{\pathimg{h}{s}{\pi}}$. This directly means
(through the fact that $\pathimg{h}{s}{\pi} = \pathimg{h}{\tilde s}{\rho}$) that
$\hat s(\argm_i)\notin\locn{\pathimg{h}{\tilde s}{\rho}}$, completing the proof.
If $\pi'=\pi$, then $\eval{\e_i}{s,h}=h(\lval{\pi}{s,h})$.
By acyclicity of $\pi$ in $s,h$ we have that
$\hat s(\argm_i)=h(\lval{\pi}{s,h})\notin\locn{\pathimg{h}{s}{\pi}}=\locn{\pathimg{h}{\tilde s}{\rho}}$.

* Second clause of disconnectedness: We need to show that for all $\alpha\in\dom{h}$, if
$h(\alpha)\in\locn{\pathimg{h}{\tilde s}{\rho}}$, then
$\alpha\in\dom{\pathimg{h}{\tilde s}{\rho}}$, or, that if
$h(\alpha)\in\locn{\pathimg{h}{s}{\pi}}$, then
$\alpha\in\dom{\pathimg{h}{s}{\pi}}$. This is immediate by the disconnectedness
of $\pi$ in $s,h$.

Acyclicity of $\rho$ in $\tilde s, h$ follows from Lemma <ref>,
noting that $\pathimg{h}{s}{\pi} = \pathimg{h}{\tilde s}{\rho}$ as shown above.

We now can apply Lemma <ref> to
obtain that $\rho$ is preserved from $\tilde s,h$ to $\tilde s',h'$ (note that since there are no
accesses to $\argm_{n+1}$ in $\body{\m}$, and no procedure calls, trivially
$\forall\e\in\alpha(\set{\tr''})_W.\, \e\not\properprefix\rho$).

Lemma <ref> yields
$\lval{\rho}{\tilde s',h'} = \lval{\pi}{s,h'}$,
and $\pathimg{h'}{\tilde s'}{\rho} = \pathimg{h'}{s}{\pi}$,
therefore also $\dom{\pathimg{h'}{\tilde s'}{\rho}} = \dom{\pathimg{h'}{s}{\pi}}$.
Lemma <ref> yields the acyclicity of $\pi$ in $s,h'$.

Finally, we also need to finally show that $\pi$ is disconnected in $s, h'$.
This follows from the fact that $\dom{\pathimg{h'}{s}{\pi}}=\dom{\pathimg{h'}{\tilde s'}{\rho}}$
and the disconnectedness of $\rho$ in $\tilde s',h'$.

Case: for exactly one $j\le n$, $\e_j\properprefix\pi$, and $\e_j= z$.
This case proceeds almost identically, but for the following differences:
we don't use $\argm_{n+1}$, nor “jump” to another trace $\tr''$ as above, but operate directly on
$\tr'$; the definitions $\tilde s \eqdef \hat s$ and $\tilde s' \eqdef \hat s'$ are trivial;
and $\rho \eqdef \argm_j.\fs$.

Application of Lemma <ref>: Note that $\alpha(\set{\cdot})$
is monotonic: $\tr_{\text{push}}\append\tr'\prefix\tr$, thus $\alpha(\set{\tr_{\text{push}}\append\tr'})_W \subseteq \alpha(\set{\tr})_W$.
Thus, from the top-level assumption about the $W$ component of the abstraction of $\tr$ follows that, also,
$\forall\e\in  \alpha(\set{\tr_{\text{push}}\append\tr'})_W.\,\e\not\properprefix\pi$.
By the definition of $\exec$ for $\cmdpush{}$, this means
$\forall\e\in \subst{\alpha(\set{\tr'})_W}{[\argm_i\mapsto\e_i]_{i=1}^n}.\,\e\not\properprefix\pi$,
$\forall\e\in \alpha(\set{\tr'})_W.\,\e\not\properprefix\rho$, by the properties of substitutions, the
definition of $\exec$ for load/stores, the fact that there are no procedure calls in $\m$, and
the fact that the root of $\rho$ is $\argm_j$ and thus not a local.

Case: for exactly one $j\le n$, $\e_j\properprefix\pi$, and $\e_j=z.\gs$ where $\gs\properprefix\fs$ and $\gs\neq\epsilon$.
Again, we don't use $\argm_{n+1}$; the definitions of $\tilde s \eqdef \hat s$
and $\tilde s' \eqdef \hat s'$ are again trivial.
Let $\hs$ be such that $\pi=z.\gs.\hs$; it is easy to see that $\hs\neq\epsilon$.
We set $\rho\defeq\argm_j.\hs$.

We define a heap $\tilde h$ as the subheap of $h$, reachable from any $\e_i$: that is, $\tilde h$ is
the smallest subheap of $h$ whose domain includes any address of the form $(\eval{\e_i}{s,h},\f)$
and which is closed under the image of $h$.  This is well-defined to the well-behavedness of $s,h$.
We fix as $h_0 \eqdef h\setminus\tilde h$, i.e., the remainder heap after removing $\tilde h$.
It should be easy to see that $\lval{\rho}{\tilde s, \tilde h}$ is defined.
By Lemma <ref>, there exists a trace $\tr''\in\sem{\body{\m}}\angled{\tilde s\cons s, \tilde h, l}$
such that $\llast{\tr''}=\angled{\tilde s'\cons s, \tilde h', L'}$ and $h'=h_0\cup\tilde h'$.

We now show that $\rho=\argm_j.\hs$ is disconnected and acyclic in $\tilde s, \tilde h$.

It is easily seen by the definition of $\tilde h$ that
$\pathimg{h}{s}{\pi} \supset \pathimg{h}{\tilde s}{\rho}=\pathimg{\tilde h}{\tilde s}{\rho}$
and $\lval{\pi}{s,h}=\lval{\rho}{\tilde s, h}=\lval{\rho}{\tilde s, \tilde h}$,
therefore also $\dom{\pathimg{h}{s}{\pi}} \supset \dom{\pathimg{h}{\tilde s}{\rho}}= \dom{\pathimg{\tilde h}{\tilde s}{\rho}}$.

* First clause of disconnectedness: We need to show that for all variables
$x\neq \argm_j$, $\tilde s(x)\notin\locn{\pathimg{\tilde h}{\tilde s}{\rho}}$,
or, that for all $i\le n$ and $i\neq j$,
$\tilde s(\argm_i) = \eval{\e_i}{s,h} \notin\locn{\pathimg{\tilde h}{\hat s}{\rho}}=\locn{\pathimg{h}{\hat s}{\rho}}\subset\locn{\pathimg{h}{s}{\pi}}$.

By the assumption of this case, for all $i\neq j$, $i\le n$ is it the case that $\e_i\not\properprefix\pi$.
If $\e_i$ is a variable, then by the disconnectedness of $\pi$ we have $\eval{\e_i}{s,h}\notin\locn{\pathimg{h}{s}{\pi}}$.
If $\e_i$ is a path, by Lemma <ref>, it follows that $\lval{\e_i}{s,h}\notin\dom{\pathimg{h}{s}{\pi}}$.
Thus, by the contrapositive of the 2nd clause of disconnectedness of $\pi$, $\eval{\e_i}{s,h}=h(\lval{\e_i}{s,h})\notin\locn{\pathimg{h}{s}{\pi}}$.
This suffices because $\dom{\pathimg{h}{s}{\pi}} \supset \dom{\pathimg{\tilde h}{\tilde s}{\rho}}$.

* Second clause of disconnectedness: We need to show that for all $\alpha\in\dom{\tilde h}$, if
$\tilde h(\alpha)\in\locn{\pathimg{\tilde h}{\tilde s}{\rho}}$, then
$\alpha\in\dom{\pathimg{\tilde h}{\tilde s}{\rho}}$. This is true by the construction of $\tilde h$
and the disconnectedness of $\pi$ in $s,h$.

Acyclicity of $\rho$ in $\tilde s, \tilde h$ follows the acyclicity of $\pi$ in $s,h$ and the fact
that $\tilde h \subseteq h$.

Application of Lemma <ref>: Note that $\alpha(\set{\cdot})$
is monotonic: $\tr_{\text{push}}\append\tr'\prefix\tr$, thus $\alpha(\set{\tr_{\text{push}}\append\tr'})_W \subseteq \alpha(\set{\tr})_W$.
Thus, from the top-level assumption about the $W$ component of the abstraction of $\tr$ follows that, also,
$\forall\e\in  \alpha(\set{\tr_{\text{push}}\append\tr'})_W.\,\e\not\properprefix\pi$.
By the definition of $\exec$ for $\cmdpush{}$, this means
$\forall\e\in \subst{\alpha(\set{\tr'})_W}{[\argm_i\mapsto\e_i]_{i=1}^n}.\,\e\not\properprefix\pi$,
$\alpha(\set{\tr'})_W=\alpha(\set{\tr''})_W$ (because $\ab{\tr'}=\ab{\tr''}$),
$\forall\e\in \subst{\alpha(\set{\tr''})_W}{[\argm_i\mapsto\e_i]_{i=1}^n}.\,\e\not\properprefix\pi$.
This means
$\forall\e\in \alpha(\set{\tr''})_W.\,\e\not\properprefix\rho$, by the properties of substitutions, the
definition of $\exec$ for load/stores, the fact that there are no procedure calls in $\m$, and
the fact that the root of $\rho$ is $\argm_j$ and thus not a local.

We now can apply Lemma <ref> to
obtain that $\rho$ is preserved from $\tilde s,\tilde h$ to $\tilde s',\tilde h'$.
Also, due to the stability assumption, it is clear that $\tilde{s}(\argm_j)=\tilde{s'}(\argm_j)$,
thus $\lval{\rho}{\tilde s',\tilde h'} = \lval{\rho}{\tilde s, \tilde h'}$.
Because $h'=h_0\cup \tilde h'$, it is easy to check that
$\lval{\rho}{\tilde s,\tilde h'} = \lval{\rho}{\tilde s, h'} = \lval{\pi}{s,h'}$.
By the same token, it can be shown that $\dom{\pathimg{h'}{s}{\pi}}=\dom{\pathimg{h}{s}{\pi}}$
and that $\pi$ is acyclic and disconnected in $s,h'$.

Let $\tr$ be a (prefix) trace produced by Lemma <ref>.
If additionally it is the case that $\forall\e\in\alpha(\set{\tr})_W.\,\e\not\properprefix\pi$
then $\lval{\pi}{s,h}=\lval{\pi'}{s',h'}$.

We prove by backwards induction on the prefixes $\tr'\prefix\tr$ we show the following:
Let $\tr'=\tr_1\append\tr_{\text{push}}\append\tr_0$ such that $\tr_{\text{push}}$ is the last
$\cmdpush{}$ command in $\tr'$ that is unmatched by a $\cmdpop{}$.  Let
$\tr_{\text{push}} = [\angled{\cmdpush{\e_1,\ldots,\e_n}, s_0\cons s_1\cons S, h_0}]$ and
$\llast{\tr'}=\angled{\_,s'\cons s_1\cons S, h'}$.
Let $\pi'$ be some path such that $\pi'[\bar{\theta}]\prefix\pi$, where $\bar{\theta}=\theta_0\cons\_$ is the
sequence of substitutions induced by the active calls at the end of $\tr'$.
Then, $\lval{\pi'}{s',h'}=\lval{\e_j.\fs}{s_0,h_0}$, where $\fs$ are such that $\pi'[\theta_0]=\e_j$.

By the definition of $\exec$, the assumption $\forall\e\in\alpha(\set{\tr})_W.\,\e\not\properprefix\pi$
and the fact that $\pi'[\bar{\theta}]\prefix\pi$, we can conclude that
Thus, the trace $\tr_0$ satisfies the conditions of Lemma <ref>, thus
$\pi'$ is preserved from $s_0,h_0$ to $s_0,h'$.  It is easy to see that
$s_0(\rootvar{\pi'})=s'(\rootvar{\pi'})$, thus by Lemma <ref>
$\pi'$ is preserved from $s_0,h_0$ to $s',h'$. By the observation that $s_0(\argm_j)=\eval{\e_j}{s,h_0}$,
and setting $\fs$ to be such that $\pi'[\theta_0]=\e_j$, it can be seen that

At the last step, $\pi'=\pi$ and we are done.

Given paths $\pi_1=x.\fs$ and $\pi_2=y.\fs$, there exists a well-behaved
state $(s_1,s_2,h)$ such that $\pathimg{h}{s_1}{\pi_1}=\pathimg{h}{s_2}{\pi_2}$,
$\lval{\pi_1}{s_1,h}=\lval{\pi_2}{s_2,h}$ and $\pi_i$ is acyclic and
disconnected in $s_i,h$, for $i\in\{1,2\}$.

Let $\fs = \f_1.\cdots.\f_n$ for some $n\ge 1$. Let $\ell_0,\ell_1,\ldots,\ell_n$ be distinct locations,
and set $s_1(x)=s_2(y)=\ell_1$.  Also, set $s_1(x')=s_2(y')=\ell_0$ for all $x'\in\Var\setminus\{x\}$
and $y'\in\Var\setminus\{y\}$. This defines $s_1,s_2$.

Next, set $h(\ell_0,\f)=\ell_0$ for all $\f\in\FieldName$.
For each $1\le j < n$, set $h(\ell_j, f_j)  = \ell_{j+1}$ and
$h(\ell_j, \f) = \ell_0$ for all $\f\in\FieldName\setminus\{\f_j\}$. Set $h(\ell_n,\f) = \ell_0$
for all $\f\in\FieldName$. This determines $h$.

It is easy to check that all required properties are true of $(s_1,s_2,h)$.

For any trace $\tr\in\sem{C} \astate$ and any two positions $i\le j \le |\tr|$, such that
$\tr_i=\angled{\_,\_,h,\_}$ and $\tr_j=\angled{\_,\_,h',\_}$, it is the case that

By the concrete semantics and the fact that there is no deallocation.

Let $\pi_1=x.\fs$ and $\pi_2=y.\fs$ (ie., paths with the same field sequence).
Let $s_1,s_2$ be two stacks such that $s_1(x)=s_2(y)$.
Let $h,h'$ be two heaps such that $\dom{h}\subseteq\dom{h'}$ and, in addition,
$\forall\alpha\in\dom{\pathimg{h}{s_1}{\pi_1}}\setminus\{\lval{\pi_1}{s_1,h}\}.\, h(\alpha)=h'(\alpha)$.
Assume $\pi_1$ is acyclic in $s_1,h$.
Then, $\lval{\pi_2}{s_2,h} = \lval{\pi_2}{s_2,h'}$.

We apply Lemma <ref> letting $s=s_1$, $s'=s_2$ and $h=h'$,
obtaining $\pathimg{h}{s_1}{\pi_1}=\pathimg{h}{s_2}{\pi_2}$.

As $\pi_1$ is acyclic in $s_1,h$ and $s_1(x)=s_2(y)$, we can easily
conclude that $\pi_2$ is acyclic in $s_2,h$.

we can conclude
$\dom{\pathimg{h}{s_2}{\pi_2}} \subseteq \dom{h'}$,
since $\pathimg{h}{s_1}{\pi_1}\subseteq\dom{h}\subseteq\dom{h'}$.

From the facts
$\forall\alpha\in\dom{\pathimg{h}{s_1}{\pi_1}}\setminus\{\lval{\pi_1}{s_1,h}\}.\, h(\alpha)=h'(\alpha)$
we obtain by simple substitution
$\forall\alpha\in\dom{\pathimg{h}{s_2}{\pi_2}}\setminus\{\lval{\pi_2}{s_2,h}\}.\, h(\alpha)=h'(\alpha)$.

We apply Lemma <ref> letting $\pi=\pi_2$, $s_2=s=s'$
and $h,h'$ as here to obtain $\lval{\pi_2}{s_2,h} = \lval{\pi_2}{s_2,h'}$.

Let $\pi_1=x.\fs$, $\pi_2=y.\fs$ and $s_1(x)=s_2(y)$ such that $\lval{\pi_1}{s_1,h}$ is defined.
Then for any $\gs\prefix\fs$, $\lval{x.\gs}{s_1,h}=\lval{y.\gs}{s_2,h}$, and

By induction on the length of $\gs$, the fact that $s_1(x)=s_2(y)$ and that $h$ is a function.

Let $\pi_1=x.\fs$ and $\pi_2=y.\fs$ such that $s_1(x)=s_2(y)$. Let also $h,h'$ be such that
If $\pi_1$ is disconnected in $s_1,h'$ then so is $\pi_2$ in $s_2,h'$.

We need to show that (a) for any $v\neq y$, $s_2(v)\notin\locn{\pathimg{h'}{s_2}{\pi_2}}$, and
(b) that for all $\alpha\in\dom{h'}$, if $h'(\alpha)\in\locn{\pathimg{h'}{s_2}{\pi_2}}$ then

For (a), note that $\pi_2$ is disconnected in $s_2,h$ so for any
$v\neq y$, $s_2(v)\notin\locn{\pathimg{h}{s_2}{\pi_2}}$.
But $\locn{\pathimg{h}{s_2}{\pi_2}}=\locn{\pathimg{h'}{s_2}{\pi_2}}$,
because $\dom{\pathimg{h}{s_2}{\pi_2}}=\dom{\pathimg{h'}{s_2}{\pi_2}}$, completing the proof.

For (b) observe that by Lemma <ref>,
$\pathimg{h'}{s_1}{\pi_1}=\pathimg{h'}{s_2}{\pi_2}$, thus we only need to show that
for all $\alpha\in\dom{h'}$, if $h'(\alpha)\in\locn{\pathimg{h'}{s_1}{\pi_1}}$ then
$\alpha\in\dom{\pathimg{h'}{s_1}{\pi_1}}$.  This is immediate by the disconnectedness of $\pi_1$
in $s_1,h'$.

Let $\pi_1=x.\fs$ and $\pi_2=y.\fs$ such that $s_1(x)=s_2(y)$.
If $\pi_1$ is acyclic in $s_1,h$ then so is $\pi_2$ in $s_2,h$.

Since $\pi_1$ is acyclic in $s_1,h$ we know that for any $\fs''\prefix\fs\prefix\fs$,
By Lemma <ref>, we know that
$\lval{x.\fs'}{s_1,h}=\lval{y.\fs'}{s_2,h}$ and $\lval{x.\fs''}{s_1,h}=\lval{y.\fs''}{s_2,h}$,
thus $h(\lval{y.\fs'}{s_2,h})\neq\locn{\lval{y.\fs''}{s_2,h}}$, as required.

Let $C_1,C_2$ be two programs such that $\asem{C_i} \bota= \angled{W_i,\_,A_i}$.
Let $\pi_i$ be two paths and $c_i$ two commands such that $\angled{c_i,L_i}\in A_i$, and,

* $\pi_1=v_1.\fs$ and $\pi_2=v_2.\fs$ (ie., the field sequences are the same);

* $\forall \e_i\in W_i.\, \e_i\not\properprefix \pi_i$ (for $i\in\set{1,2}$);

* $c_1=(\pi_1:=x)$ and $c_2=(\pi_2:=y)$, or, $c_1=(\pi_1:=x)$ and $c_2=(y:=\pi_2)$;

* $L_1+L_2 \le 1$.

Then there exists a well-behaved state $\astate=\angled{(s_1,s_2),h,(0,0)}$  and a
concurrent trace $\tr^{\parallel} \in \sem{C_1\parallel C_2}~\astate$ that races.

we take $L_1=0$, as at least one of $L_1,L_2$ must be zero.

Given the two paths share the same field accesses ($\fs$), Lemma <ref> provides a
a well-behaved state $\astate=\angled{(s_1,s_2),h,(0,0)}$ such that
$\lval{\pi_1}{s_1,h}=\lval{\pi_2}{s_2,h}$, and where $\pi_i$ is disconnected and acyclic in

By assumption, $\asem{C_1} \bota= \angled{W_1,\_,A_1}$, thus by Theorem <ref>,
$\angled{W_1,\_,A_1} = \alpha(\sem{C_1} \angled{s_1, h, 0})$.
Thus Lemma <ref> implies there is a trace in $\sem{C_1} \angled{s_1, h, 0}$
and a shortest prefix $\tr_1$ of that trace, whose last state is $\angled{c'_1,\_,\_,L_1}$, and:

* $c'_1$, $c_1$ are both stores or loads;

* $\exec~\ab{\tr_1}~(\bota,\epsilon)=(\angled{\_,\_,A'_1},\bar{\theta_1})$ where $\pi_1\in A'_1$;

* a path $\pi'_1$ such that $c' = (x:=\pi'_1)$ or $c'=(\pi'_1:=x)$,
and $\set{\pi_1}=\subst{\set{\pi'_1}}{\bar{\theta_1}}$.

Define the following functions. Each maps a single-thread state to a two-thread trace,
with the given parameter as the stack of the "other" thread, and zero as its lock context.

\[
\begin{array}{r@{\;}l}
\mathsf{left}(\angled{c, S, h, L}, S')  & = \angled{c\parallel\epsilon, (S, S'), (L,0)} \\
\mathsf{right}(\angled{c, S, h, L}, S')  & = \angled{\epsilon\parallel c, (S', S), (0,L)}
\end{array}
\]

We construct a concurrent trace $\tr^{\parallel}_1$ and its last state as follows:

\[\begin{array}{r@{\;}l}
\tr^{\parallel}_1 & \eqdef \left[ \mathsf{left}(\astate_i, s_2) \mid \astate_i \text{ is the $i$-th state in $\tr_1$}, 1\le i< |\tr_1| \right] \\[3pt]
\llast{\tr^{\parallel}_1} & = \angled{\_,(s'_1 \cons S'_1,s_2),h',(L'_1,0))}
\end{array}\]

that is, $C_1$ progresses along $\tr_1$, but without executing the instruction $c'_1$,
while $C_2$ never takes a step.
Observe that by construction, $L'_1=0$, as $c'_1$ is not an $\cmdunlock$ instruction (thus $L'_1=L_1$),
and by assumption $L_1=0$.

We show that the state $s_2,h'$ is well-behaved. We know that $s_1,h'$ is well-behaved, by
Lemma <ref>. Thus, for any address
$\alpha\in\dom{h'}$, $h'(\alpha)\in\locn{h'}$, taking care of the second clause of well-behavedness.
For the first clause, observe that $s_2,h$ is well-behaved thus for any variable $x$,
$s_2(x)\in\locn{h}$.  But as $\dom{h}\subseteq\dom{h'}$ by Lemma <ref>,
it follows that $\locn{h}\subseteq \locn{h'}$, thus $s_2(x)\in\locn{h'}$ too.

Thus, we have $\asem{C_2} \bota= \angled{W_2,\_,A_2} = \alpha(\sem{C_2} \angled{s_2, h', 0})$.
By identical reasoning as with $C_1$ we obtain a trace in $\sem{C_2} \angled{s_2, h', 0}$
and a shortest prefix trace $\tr_2$ of that,
whose last state is $\angled{c'_2,\_,\_,L_2}$ where

* $c'_2$, $c_2$ are both stores or loads;

* $\exec~\ab{\tr_2}~(\bota,\epsilon)=(\angled{\_,\_,A'_2},\bar{\theta_2})$ where $\pi_2\in A'_2$;

* a path $\pi'_2$ such that $c' = (x:=\pi'_2)$ or $c'=(\pi'_2:=x)$,
and $\set{\pi_2}=\subst{\set{\pi'_2}}{\bar{\theta_2}}$.

We define $\tr^{\parallel}_2$ and its last state as follows:

\[\begin{array}{r@{\;}l}
\tr^{\parallel}_2 & \eqdef \left[ \mathsf{right}(\astate_i,s'_1\cons S'_1) \mid \text{where $\astate_i$ is the $i$-th state in $\tr_2$}, 1\le i < |\tr_2| \right] \\
\llast{\tr^{\parallel}_2} & = \angled{\_,(s'_1 \cons S'_1,s'_2\cons S'_2),h'',\_)}
\end{array}\]

Finally, we define

\[
\tr^{\parallel}\eqdef\tr^{\parallel}_1\append\tr^{\parallel}_2
\]

We argue that
$\tr^{\parallel}$ is in $\sem{C_1\parallel C_2}~\angled{(s_1,s_2),h,(0,0)}$.
This is straightforward by our concurrent trace semantics (Figure <ref>)
and the observation that as $\tr^{\parallel}_1$ ends at a state with lock context $(0,0)$,
the side-condition for lock commands in $\tr^{\parallel}_2$ is satisfied by construction.
[Figure <ref> defines a prefix-closed semantics so we need not concern ourselves
with what happens after the race.]

We now show that $\lval{\pi'_1}{s'_1,h''}=\lval{\pi'_2}{s'_2,h''}$, establishing that the next
instruction following $\tr^{\parallel}$ will necessarily access the same address.
By construction we have

\begin{equation}
\begin{array}{r@{\,}l}
\lval{\pi_1}{s_1,h} & =\lval{\pi_2}{s_2,h}\\
\pathimg{h}{s_1}{\pi_1} & = \pathimg{h}{s_2}{\pi_2}\\
\multicolumn{2}{c}{\text{$\pi_i$ is acyclic and disconnected in $s_i,h$, for $i\in\set{1,2}$}}
\end{array}
\label{eq:pi1-pi2}
\end{equation}

By assumption, $\forall \e_i\in W_i.\, \e_i\not\properprefix \pi_i$.
We apply Lemma <ref> on $\tr_1$ to obtain

\begin{equation}
\pi_1 \text{ is preserved from } s_1,h \text{ to } s_1,h' \label{eq:pi-h-hprime}
\end{equation}

Lemma <ref> gives us

\begin{equation}
\lval{\pi_1}{s_1,h'}=\lval{\pi'_1}{s'_1,h'} \label{eq:pi1-endpoint}
\end{equation}

From (<ref>) we have
$\forall\alpha\in\dom{\pathimg{h}{s_1}{\pi_1}}\setminus\{\lval{\pi_1}{s_1,h}\}.\, h(\alpha)=h'(\alpha)$.
By (<ref>) we have $\pathimg{h}{s_1}{\pi_1}=\pathimg{h}{s_2}{\pi_2}$ and $\lval{\pi_1}{s_1,h}=\lval{\pi_2}{s_2,h}$,
thus $\forall\alpha\in\dom{\pathimg{h}{s_2}{\pi_2}}\setminus\{\lval{\pi_2}{s_2,h}\}.\, h(\alpha)=h'(\alpha)$.
In addition, by (<ref>), $\pi_2$ is acyclic in $s_2,h$.
Clearly, $\dom{h}\subseteq\dom{h'}$ by Lemma <ref>.
Thus Lemma <ref> gives

\begin{equation}
\lval{\pi_2}{s_2,h}=\lval{\pi_2}{s_2,h'} \label{eq:pi2-equal}
\end{equation}

The state $s_2,h'$ is well-behaved as shown above.
Since $\pi_1$ is acyclic in $s_1,h'$ (<ref>),
by Lemma <ref>, $\pi_2$ is also acyclic in $s_2,h'$.
Observe that $\dom{\pathimg{h}{s_2}{\pi_2}}=\dom{\pathimg{h}{s_1}{\pi_1}}$ (<ref>),
$\dom{\pathimg{h}{s_1}{\pi_1}}=\dom{\pathimg{h'}{s_1}{\pi_1}}$ (<ref>), and
by Lemma <ref> and the fact that $s_2(\rootvar{\pi_2})=s_1(\rootvar{\pi_1})$,

\begin{equation}
\pathimg{h'}{s_1}{\pi_1}=\pathimg{h'}{s_2}{\pi_2} \label{eq:pi1-pi2-hprime}
\end{equation}

thus $\dom{\pathimg{h}{s_2}{\pi_2}}=\dom{\pathimg{h'}{s_2}{\pi_2}}$.
Lemma <ref> applies, therefore $\pi_2$ is disconnected in $s_2,h'$.

\begin{equation}
\pi_2 \text{ is preserved from } s_2,h \text{ to } s_2,h' \label{eq:pi2-h-hprime}
\end{equation}

This means we can apply Lemma <ref> to $\tr_2$ and obtain that

\begin{equation}
\pi_2 \text{ is preserved from } s_2,h' \text{ to } s_2,h'' \label{eq:pi2-hprime-hprimeprime}
\end{equation}

and Lemma <ref> yields

\begin{equation}
\lval{\pi_2}{s_2,h'}=\lval{\pi'_2}{s'_2,h''} \label{eq:pi2-endpoint}
\end{equation}

Chaining (<ref>), (<ref>), (<ref>), (<ref>)
and (<ref>),
we conclude $\lval{\pi'_1}{s'_1,h''}=\lval{\pi'_2}{s'_2,h''}$, thus establishing the
existence of a race as per Definition <ref>.


[a]Avik Roy

# Interpretability of an Interaction Network for identifying $H\rightarrow
b\bar{b}$ jets

Mark S. Neubauer

###### Abstract

Multivariate techniques and machine learning models have found numerous
applications in High Energy Physics (HEP) research over many years. In recent
times, AI models based on deep neural networks are becoming increasingly
popular for many of these applications. However, neural networks are regarded
as black boxes- because of their high degree of complexity it is often quite
difficult to quantitatively explain the output of a neural network by
establishing a tractable input-output relationship and information propagation
through the deep network layers. As explainable AI (xAI) methods are becoming
more popular in recent years, we explore interpretability of AI models by
examining an Interaction Network (IN) model designed to identify boosted $H\to
b\bar{b}$ jets amid QCD background. We explore different quantitative methods
to demonstrate how the classifier network makes its decision based on the
inputs and how this information can be harnessed to reoptimize the model-
making it simpler yet equally effective. We additionally illustrate the
activity of hidden layers within the IN model as Neural Activation Pattern
(NAP) diagrams. Our experiments suggest NAP diagrams reveal important
information about how information is conveyed across the hidden layers of deep
model. These insights can be useful to effective model reoptimization and
hyperparameter tuning.

##  1 Introduction

Owing to their intricate internal structure, neural networks (NNs) have often
been treated as black boxes. It is difficult to understand how different input
features contribute to the network’s computational process and how the inter-
connected neural pathways convey information. In recent years, advances in
explainabale Artificial Intelligence (xAI) [1] have made it possible to build
intelligible relationship between an AI model’s inputs, architecture, and
predictions [2]. xAI has been successful in learning the underlying physics of
a number of problems in high energy detectors [3]. In this work, we apply
state-of-the art xAI techniques in interpreting an Interaction Network (IN)
model [4] developed to identify boosted $H\to b\bar{b}$ jets from QCD
background.

##  2 Evaluating Feature Importance for the IN Model

Figure 1 shows the IN model architecture and tabulates the default
hyperparameters and data dimensions. This network is built to train on graph
data structure whose nodes comprise of $N_{p}$ particle tracks, each with $P$
features, and $N_{v}$ secondary vertices, each with $S$ features, associated
with the jet. The physical description of each feature is given in Appendix C
of ref. [4]. It creates a fully connected directed graph with
$N_{pp}=N_{p}(N_{p}-1)$ edges for the particle tracks. A separate graph with
$N_{vp}=N_{v}N_{p}$ generates all possible connections between the particle
tracks and the secondary vertices.

Default IN hyperparameters
---
$(P,N_{p},S,N_{v})$ | $(30,60,14,5)$
Hidden layers | 3
Hidden layer | 60
dimension
$(D_{e},D_{o})$ | (20, 24)
Activation | ReLU

Figure 1: A schematic diagram of the network architecture and dataflow in the
IN model. This image is taken from Ref. [4]. The choice of model
hyperparameters and input data dimensions for the baseline model is given in
the accompanying table.

The node level features for the track-track (track-vertex) graph are
transformed to edge level features via a couple of interaction matrices,
identified as $R_{R[N_{p}\times N_{pp}]}$ and $R_{S[N_{p}\times N_{pp}]}$
($R_{K[N_{p}\times N_{vp}]}$ and $R_{V[N_{v}\times N_{vp}]}$). These edge-
level features are transformed via fully connected NNs, respectively called
$f_{R}^{pp}$ and $f_{R}^{vp}$, to obtain two $D_{e}$ dimensional internal
state representation of these graphs.

These internal state edge-level representations are transferred back to the
track-level representations by the interaction matrices. These particle-level
internal state representations are concatenated with the original track
features creating a $(P+2D_{e})\times N_{p}$ dimensional feature space. The
trainable dense MLP $f_{O}$ creates the post-interaction $D_{O}$ dimensional
internal representation that are summed over the tracks and then linearly
combined to produce a two-dimensional output, which is transformed to
individual class probabilities via a softmax function.

In order to realize which features play the most important role in the IN’s
decision making process, first we train the model with its default settings
(the baseline model). We mask one feature at a time for all input tracks or
secondary vertices by replacing the corresponding entries by zero values. We
obtain inference from the trained model and evaluate the Area Under the Curve
(AUC) for Region Operator Characteristic (ROC) curve (ROC AUC score) from the
model inference. The change observed in the AUC score for masking each of the
features can be seen in Figure 2. It shows that many of these input features
have a rather small impact on the model’s overall performance, reflected by
the very small change in AUC score.

Figure 2: Change in AUC score with respect to the baseline model when each of
the track and secondary vertex features is individually masked during
inference with the trained baseline model.

We can inspect the importance of these features for individual tracks and
vertices by the Layerwise Relevance Propagation (LRP) technique [5, 6]. Since
some of the input features show high degree of correlation with each other, we
use the LRP-$\gamma$ method described in ref. [6], which is designed to skew
the LRP score distributions to nodes with positive weights in the network. In
order to apply the LRP method for the IN model, we defined custom propagation
functions for (i) the aggregation of internal representation and (ii) the
transformation via interaction matrices.

We show the average scores attributed to the different features for QCD and
$H\rightarrow b\bar{b}$ jets in figure 3a. When compared with the change in
AUC score by individual features in figure 2, the features with largest
relevance scores coincide with the features that individually cause the
largest drop in AUC score. We additionally observe that vertices features play
a more important role in identifying the $H\rightarrow b\bar{b}$ jets. This
behavior is also justified from a physics stand point, since the presence of
high energy secondary vertices is an important signature for jets from $b$
quark because of its relatively longer lifetime.

(a)

(b)

Figure 3: Average relevance scores attributed to LABEL:sub@fig:LRP-feats input
track and secondary vertex features and LABEL:sub@fig:LRP-trackNverts
individual tracks and secondary vertices. The tracks and secondary vertices
are ordered according to their relative energy with respect to the jet energy.

However, the approaches also show some inconsistencies among themselves. The
secondary vertex features sv_ptrel and sv_erel are assigned very low relevance
scores while masking them independently cause very large drops in the AUC
score. These variables are highly correlated and both display very large
correlation (correlation coefficient of 0.85) with sv_pt (Figures 4a and 4b).
The LRP-$\gamma$ method skews their relevance distribution and suppresses the
LRP scores for those two variables while assigning a larger score to the
variable sv_pt.

(a)

(b)

(c)

Figure 4: Scatter plot of LABEL:sub@fig:svptrel sv_ptrel and sv_pt and
LABEL:sub@fig:sverel sv_ptrel and sv_erel. LABEL:sub@fig:trackQ shows the
distribution of the categorical variable track_quality.

We make an additional observation regarding the importance attributed to the
feature called track_quality. This feature is essentially a qualitative tag
regarding the track reconstruction status, and for most of the training data,
this has almost identical distribution for both jet categories (Figure 4c).
With such an underlying distribution, it is obvious that this variable doesn’t
contribute to the classifier’s ability to tell apart the jet categories.
However, the large relevance score associated with it indicates that the
classifier’s class-predictive output for each class somehow receives a large
contribution from its numerical embedding. The model that was trained without
these variables, along with the 11 (3) track (vertex) features that report a
change in AUC of less than $0.01\%$ converged with an AUC score of $99.00\%$,
performing as equally well.

##  3 Inspecting the Activation Layers and Model Reoptimization

As the IN processes the input, it is passed through three different MLPs that
approximate arbitrary non-linear functions identifies as
$f_{R}^{pp},f_{R}^{vp},$ and $f_{O}$. In order to explore the activity of each
neuron and compare it with the activity of neurons in the same layer, we
define the quantity Relative Neural Activity (RNA) [7] as
$\mathrm{RNA}(j,k;\mathcal{S})=\frac{\sum_{i=1}^{N}a_{j,k}(s_{i})}{\max_{j}\sum_{i=1}^{N}a_{j,k}(s_{i})}$
where $\mathcal{S}=\\{s_{i}\\}$ represents a set of samples over which the
$\mathrm{RNA}$ score is evaluated. The quantity $a_{j,k}(s_{i})$ is the
activation of $j$-th neuron in the $k$-th layer when the input to the network
is $s_{i}$. Figure 5 shows the neural activation pattern (NAP) diagram for the
baseline model, showing the $\mathrm{RNA}$ scores for the different activation
layers. The scores are separately evaluated for QCD and $H\rightarrow
b\bar{b}$ categories. To simultaneously visualize these scores, we project the
$\mathrm{RNA}$ scores of the former as negative values. The NAP diagram
clearly shows that the network’s activity level is quite sparse- while some
nodes are playing very important roles in propagating the necessary
information, other nodes don’t participate as much. We additionally observe
that the right until the very last layer of the aggregator network $f_{O}$,
the same nodes show the largest activity level for both jet categories.
However, different nodes are activated in the last layer for the two jet
categories, indicating an effective disentanglement of the jet category
information in this layer.

Figure 5:  2D map of $\mathrm{RNA}$ score for different nodes of the
activation layers. To simultaneously visualize the scores for QCD and
$H\rightarrow b\bar{b}$ jets, we project the $\mathrm{RNA}$ scores of the
former as negative values

The sparsity of NAP diagram and low feature importance for a number of input
features for the baseline IN model indicates that the model can be made
simpler, by both reducing the number of input features it relies on and the
number of trainable parameters. To demonstrate this, we train some alternate
variants of the IN models where we drop track_quality, sv_ptrel, sv_erel along
with the 11 (3) track (vertex) features that report a change in AUC of less
than $0.01\%$. The details and performance metrics of these models are given
in Table 1. They demonstrate that the baseline model can be made much simpler
without compromising the quality of its performance. As can be seen from the
results in Table 1, the ROC-AUC score of the alternate models are very close
to that of the baseline model, though the number of trainable parameters is
significantly lower.

$\Delta P$, $\Delta S$ | $h,D_{e},D_{o}$ | Parameters | AUC score | Sparsity
---|---|---|---|---
0, 0 (baseline) | 60, 20, 24 | 25554 | 99.02% | 0.56
| 32, 16, 16 | 8498 | 98.87% | 0.52
12, 5 | 32, 8, 8 | 7178 | 98.84% | 0.48
| 16, 8, 8 | 2842 | 98.62% | 0.40
Table 1: The performance of the baseline and alternate retrained models with
modified hyperparameters. Sparsity is measured by the fraction of activation
nodes with an RNA score less than 0.2

##  4 Conclusion

In this paper we have demonstrated how the application of xAI methods aided
with physical intuitions can help identify important features for the task of
identifying $H\to b\bar{b}jets$. We additionally propose a novel metric, the
RNA score, and an associated visualization tool, the NAP diagram, to
investigate information propagation through a model. These tools help
understand the sparsity of information propagation and hence optimize model
complexity without degrading the model’s performance.

## Acknowledgements

This work was supported by the U.S. Department of Energy (U.S. DOE), Office of
Science, High Energy Physics, under contract number DE-SC0023365 and the FAIR
Data program of the U.S. DOE, Office of Science, Advanced Scientific Computing
Research, under contract number DE-SC0021258.

## References

  * [1] T. Miller, _Explanation in artificial intelligence: Insights from the social sciences_ , _Artificial Intelligence_ 267 (2019) 1.
  * [2] D. Gunning, M. Stefik, J. Choi, T. Miller, S. Stumpf and G.-Z. Yang, _XAI—Explainable artificial intelligence_ , _Science Robotics_ 4 (2019) eaay7120.
  * [3] D. Turvill, L. Barnby, B. Yuan and A. Zahir, _A survey of interpretability of machine learning in accelerator-based high energy physics_ , in _2020 IEEE/ACM International Conference on Big Data Computing, Applications and Technologies (BDCAT)_ , pp. 77–86, 2020, DOI.
  * [4] E.A. Moreno, T.Q. Nguyen, J.-R. Vlimant, O. Cerri, H.B. Newman, A. Periwal et al., _Interaction networks for the identification of boosted h→ b b decays_ , _Physical Review D_ 102 (2020) 012010.
  * [5] A. Binder, S. Bach, G. Montavon, K.-R. Müller and W. Samek, _Layer-wise relevance propagation for deep neural network architectures_ , in _Information science and applications (ICISA) 2016_ , pp. 913–922, Springer (2016), DOI.
  * [6] G. Montavon, A. Binder, S. Lapuschkin, W. Samek and K.-R. Müller, _Layer-wise relevance propagation: an overview_ , _Explainable AI: interpreting, explaining and visualizing deep learning_ (2019) 193.
  * [7] A. Khot, M.S. Neubauer and A. Roy, _A detailed study of interpretability of deep neural network based top taggers_ , _arXiv preprint arXiv:2210.04371_ (2022) .