Astroparticle and particle physics at ultra-high energy: results from the
Pierre Auger Observatory

Antonella Castellina$\star$, for the Pierre Auger Collaboration$\star\star$

* INFN, Sezione di Torino and Osservatorio Astrofisico di Torino (INAF), Torino, Italy

** Full author list: http://www.auger.org/archive/authors_2022_05.html

E-mail<EMAIL_ADDRESS>
21st International Symposium on Very High Energy Cosmic Ray Interactions
(ISVHE- CRI 2022)
Online, 23-27 May 2022 10.21468/SciPostPhysProc.?

## Abstract

The scientific achievements of the Pierre Auger Collaboration cover diverse
and complementary fields of research. The search for the origin of ultra-high
energy cosmic rays (UHECRs) is based on the measurement of the energy spectrum
and mass composition of the primaries, on studies of multi-messengers, and on
extensive anisotropy searches. With the collected data it is also possible to
explore the characteristics of hadronic interactions at energies unreachable
at human-made accelerators, and to assess the existence of non-standard
physics effects. A selection of the latest results is presented and the
emerging picture is discussed.

###### Contents

  1. 1 Introduction
  2. 2 Experimental results
    1. 2.1 Arrival directions
    2. 2.2 Energy spectrum and composition
    3. 2.3 Multi-messengers
  3. 3 Constraining the sources of UHECRs
  4. 4 UHECRs and particle physics
    1. 4.1 Hadronic interactions
    2. 4.2 Beyond the standard model
  5. 5 Conclusion

##  1 Introduction

The quest for the sources of ultra-high energy cosmic rays (UHECRs) is one of
the most interesting open questions in astroparticle physics. Although we
still lack a comprehensive theory capable of explaining where the UHECRs are
produced, what is the acceleration mechanism driving nuclei to such extreme
energies and their chemical composition, a large amount of information was
gained in the recent years by exploiting the data taken with the Pierre Auger
Observatory. The largest sample of ultra-high energy events ever collected led
to major advances in our understanding of their properties, reaching a global
view which appears to be much more complex than the one available in the past.
A key role in the interpretation of the results is played by the hadronic
interaction models, which predictions, based on the most up-to-date results by
accelerators, need to be extrapolated to much higher energies. There, UHECRs
are the only probes available; besides providing new information, they can
also be exploited to explore possible effects beyond the standard model of
particle physics.

The Pierre Auger Observatory [1] has collected data for more than 15 years,
accumulating the world’s largest exposure to ultrahigh energy cosmic rays
(UHECRs). Extensive air showers induced by CR-interactions in the atmosphere
are observed using 24 fluorescence telescopes (FD) and a surface detector
(SD1500) consisting of 1600 water Cherenkov stations distributed over an area
of $3000$ km2 at mutual distance of $1500$ m. Three more FD telescopes and a
denser array of 61 detectors (SD750) are deployed in the North-Western part of
the Observatory, covering $27$ km2, to extend the measurements to lower
energies. Using the FD technique, we can measure the nearly calorimetric
energy deposited in the atmosphere by the electromagnetic component of the
EAS. The contribution of the invisible energy (carried on by muons and
neutrinos) is then added based on a data-driven method, to obtain the total
energy of the primary. Exploiting the correlation of the SD energy estimator
(the particle density at 1 km from the shower core, corrected by the
attenuation) with the FD energy in hybrid events, i.e. those detected by both
the SD and the FD, we obtain the energy scale of the Observatory and a
calibration of all SD events in an almost model independent way. The FD energy
resolution amounts to about 8%, the overall systematic uncertainty on the
energy scale to 14% [2].

##  2 Experimental results

###  2.1 Arrival directions

Searches for anisotropies in the arrival directions of UHECRs can be performed
at lower energies on large angular scales (much larger than the angular
resolution of the experiment), where the cumulative flux from different
sources could possibly be seen despite magnetic deflections. A large-scale
dipolar anisotropy in the arrival direction distribution of events recorded by
the Pierre Auger Observatory was first observed in [3], and further explored
in [4, 5], under the assumption that higher multipoles are negligible. With
more than 15 years of data, including events up to $80^{\circ}$, we could
exploit the enormous exposure of $120,000$ $\mathrm{km^{2}~{}sr~{}yr}$ to
reach a post-trial significance of $6.6\sigma$. The dipole points $\sim
125^{\circ}$ away from the Galactic centre direction, thus suggesting an
extra-galactic nature of the UHECRs above this energy. A Galactic origin would
result in a much stronger dipole, in a direction within a few tens of degrees
of the Galactic Center.

Figure 1: Left: Flux map at energies above $8$ EeV, smoothed in windows of
$45^{\circ}$. Right: Flux map at energies above $40$ EeV with a top-hat
smoothing radius of $\Psi=25^{\circ}$.

For the energies below $8\times 10^{18}$ eV, the result is consistent with
isotropy, the upper bound to the amplitude being about $1.2\%$ at $95\%$ C.L.
The dipole amplitude is found to grow with energy; this could for instance be
associated with the larger relative contribution to the flux that arises at
high energies from nearby sources. A suppression of the more isotropic
contribution from farther ones is expected to result from the strong
attenuation of the CR flux that should take place at the highest energies as a
consequence of their interactions with the background radiation. This energy
dependence can be reproduced only with an evolution of the mass composition
[6, 7].

At higher energies, intermediate scale anisotropies can be searched for, as
Galactic and extra-galactic deflections of UHECR are expected to be small
enough for point sources to be visible as warm/hot spots. The largest
available set of events at such energies was obtained thanks to an integrated
exposure of $122,000$ $\mathrm{km^{2}~{}sr~{}yr}$. Evidence of a deviation
from isotropy has been found by correlating the arrival directions of UHECRs
with the direction of starburst galaxies and active galactic nuclei (AGN)
above $\sim 40\times 10^{18}$ eV, at $4.2\sigma$ and $3.3\sigma$ significance
level respectively [8]. This result has been derived accounting for the
attenuation of the UHECR nuclei, which mixture was inferred from lower energy
observations. The observed correlation does not allow us to disentangle the
different classes of objects yet; it could furthermore be altered by the
effects of Galactic and extra-Galactic magnetic fields.
The overall fluxes in the two energy bins are shown in Fig.1.

A joint collaboration between the Pierre Auger Observatory and the Telescope
Array (TA), a hybrid experiment operating since 2008 in the Northern
hemisphere and covering an area of $700$ km2, allows us to explore possible
deviations from isotropy over the full sky. The search for large scale
anisotropies confirmed the Auger-only results, reducing the uncertainties on
the North-South components of the dipole and quadrupole components [9].
Studies on intermediate scale anisotropies found excesses lying along the
Supergalactic Plane; $11.8\%^{+5.0\%}_{-3.1\%}$ of the events show a
correlation with the position of nearby starburst galaxies on a $15^{\circ}$
angular scale, with a post-trial significance of $4.2\sigma$ [10].

###  2.2 Energy spectrum and composition

A huge exposure of about $80,000$ $\mathrm{km^{2}~{}sr~{}yr}$ allowed us to
measure the cosmic ray energy spectrum exploiting more than 800,000 events
collected by the SD1500 and the SD750 arrays with energies above $10^{17}$ eV
[11]. The ankle and the flux suppression are clearly identified at $(4.9\pm
0.1\pm 0.8)\times 10^{18}$ eV and $(4.7\pm 0.3\pm 0.6)\times 10^{19}$ eV
respectively, as shown in the left panel of Fig.2. A new feature dubbed the
instep could be identified at $1.4\times 10^{19}$ eV.
The energy at which the integral spectrum drops by a factor of two below what
would be that expected with no steepening is $E_{1/2}=(23\pm 4)\times 10^{18}$
EeV, at variance with that expected for uniformly distributed sources of
protons, as also indicated by the results from our composition measurements.

The measurement of the distributions of the depth of maximum development of
the showers in the atmosphere provides us with the most reliable information
on the UHECRs composition. Indeed, the average $\left<X_{\mathrm{max}}\right>$
is directly proportional to the logarithmic mass of the primary particle,
while its variance is a convolution of the intrinsic shower-to-shower
fluctuations and of the dispersion of masses in the primary beam.
The first two moments of the depth of shower maximum as measured at the Pierre
Auger Observatory are shown in the central and rightmost panels of Fig.2.

Figure 2: Left: SD energy spectrum after combining the individual measurements
by the SD-750 and the SD-1500 scaled by $E^{2.6}$. Right: the first two
moments of the $X_{\mathrm{max}}$ distributions measured by the FD.

The observed rate of change of $\left<X_{\mathrm{max}}\right>$ with energy is
$(77\pm 2)$ and $(26\pm 2)$ $\mathrm{g~{}cm^{-2}}$/decade below
$E_{0}=10^{18.32\pm 0.03}$ eV and at higher energies respectively. The
comparison with the predictions from models (tuned to the most recent LHC
results) assuming pure primaries clearly points to a change of composition
towards lighter elements up to $E_{0}$, while above there is an increase in
the average mass. Two data sets are plotted here: in black the direct
measurement of $\left<X_{\mathrm{max}}\right>$ by the FD [12, 13], in blue the
SD result where the rise-time of the signal in the stations is used as the
observable related to mass [14]. This latter sample allows us to extend the
composition measurements to higher energies, although at the expense of larger
systematic uncertainties. A confirmation of these results comes from the
measurement of the standard deviation. The large values of
$\sigma(X_{\mathrm{max}})$ below $E_{0}$ can be interpreted as coming from
either light or mixed primaries, whereas the subsequent decrease would
corresponds to a purer intermediate to heavy composition. A completely
independent measurement of the $X_{\mathrm{max}}$ distributions was recently
obtained exploiting the AERA radio array of the Observatory, confirming the
previous conclusion [15].
The average composition of UHECRs also appears to be heavier at low than at
high Galactic latitudes, the cut being put after a pre-scan at Galactic
latitute $|b|=30^{\circ}$ [16]. The post-trial significance of this very
interesting result is currently $4.4\sigma$. The interpretation of this result
is however not straigthforward: e.g., the Galactic magnetic field could be
responsible of the mass dependence of anisotropy, but the same effect could
arise from the different horizons probed by different primaries.

It is interesting to compare our findings with those provided by TA. The
results on the flux of UHECRs are in agreement within the systematic
uncertainties up to about $3\times 10^{19}$ eV. The difference appearing at
higher energies is still under discussion [17]. From the point of view of
composition, a detailed comparison taking into account the detector effects
shows that the TA result is compatible with the mixed composition measured by
the Pierre Auger Observatory at least up to 10 EeV [18]. The lack of
statistics does not allow TA to distinguish between a mixed, Auger-like
composition and a proton one.

###  2.3 Multi-messengers

With its high potential in the identification of photons and neutrinos, the
Pierre Auger Observatory has a key role in the field of multi-messenger
astrophysics. Astrophysical neutrinos and photons can be produced in the
sources or in their environment, and cosmogenic ones are byproducts of the
interactions of UHECRs with the radiation fields permeating the Universe.
While photons cover travel distances of the order of $\sim 4-5$ Mpc at
$10^{18}$ eV, neutrinos allow us to probe sources up to cosmological
distances. Their fluxes depend on the properties of the sources and on the
composition of the primary beam. Photon showers are selected based on the fact
that they present a more elongated profile in the atmosphere, and thus a
larger $X_{\mathrm{max}}$, a steeper lateral distribution and a reduced
production of muons with respect to hadronic showers. Neutrinos are tagged by
selecting young horizontal showers with a large electromagnetic component
close to the detector; these events can only be produced by neutrinos as in
hadronic-induced cascades the electromagnetic component would be almost
completely absorbed in the atmosphere and only muons could be detected in the
SD.

Figure 3: Upper limits on the cosmogenic neutrinos (left) and photon (right)
fluxes compared to limits from other experiments and to model predictions
([19, 20] and refs.therein).

Two main categories of events are considered: Earth-skimming, induced by
$\nu_{\tau}$ travelling from below the Earth crust in directions just below
the horizon and producing a $\tau$ lepton which can then generate a shower
above the SD, and downward-going events due to neutrinos of any flavour. The
most recent results on cosmogenic neutrinos and photons are shown in Fig.3
[19, 20].
Exploiting the very good angular resolution of the Observatory (less than
$0.5^{\circ}$ for $\theta>60^{\circ}$), point-like sources of UHE neutrinos
were searched for too [21, 22]. Upper limits on the total energy per source
have been setup following up binary black hole mergers detected via
gravitational waves (GW) [23]; preliminary results on the photon fluence from
selected GW sources were obtained in [24].

##  3 Constraining the sources of UHECRs

Any astrophysical model, however simplified or complex, aiming at describing
the UHECRs sources and their characteristics must be able to describe the
particles features as measured at Earth. A comparison with the experimentally
derived energy spectrum and composition is thus mandatory.
A combined fit of a simple astrophysical model of UHECR sources to both the
energy spectrum and mass composition data measured by the Pierre Auger
Observatory was exploited to investigate the constraining power of the
collected data on the source properties, considering the energy region
$E\gtrsim 6\times 10^{17}$ eV) [26]. Starting from a quite simple
astrophysical scenario of stationary and uniform sources in a co-moving
volume, we considered several nuclear components injected at the sources with
a power-law spectrum and with the maximal energy of the sources modeled with
an exponential cut-off.
A good description of our data is obtained by considering two extragalactic
components. Above the ankle, medium mass elements are escaping from the
sources with a very hard energy spectrum and a rather low rigidity cutoff; the
flux suppression appears to be mainly due to source exhaustion, rather than to
the anyway not negligible propagation effects. The instep reflects the
interplay between the flux contributions of the helium and the CNO group
injected at the source, shaped by photo-disintegration during the propagation.
At lower energies, an additional light or light-to-intermediate very soft
extragalactic component is needed. A similar result is obtained when this
component is taken as an extragalactic proton-only one and a medium-mass
Galactic contribution is added; this secondary Galactic component, if present,
cannot be composed by heavy elements.
In Fig.4, the contribution of the different mass groups to the energy spectrum
and the relative abundances at the top of the atmosphere are shown. The rather
large uncertainty on the predicted total fluxes (brown band) is due to
combined effects of the $\pm 14\%$ shifts in the energy scale and the
uncertainty on $X_{\mathrm{max}}$, which ranges from 6 to 9
$\mathrm{g~{}cm^{-2}}$. The experimental uncertainties are large, and mainly
dominated by those affecting the measurement of $X_{\mathrm{max}}$; they must
be taken into account in any models comparing with data. As noted in [27], our
data constrain the luminosity density that continuously emitting sources must
inject into extragalactic space in UHECRs to supply the observed energy
density to $\sim 6\times 10^{44}$ erg Mpc-3 yr-1 above the ankle (for z=0).
The unknown cosmological evolution of the sources strongly influences these
results. Testing different combinations of source evolution for the two
extragalactic components allows us to exclude very strong values (m=5) for the
high energy one, as it would produce a too large flux of secondary particles
at the ankle with respect to the measured one.

Figure 4: Left: the energy spectrum. Right: the relative abundances at the top
of atmosphere. The shaded areas indicate the combined effect of systematic
uncertainties on $E$ and $X_{\mathrm{max}}$. [A=1 (red), 2 $\leq$ A $\leq$ 4
(grey), 5 $\leq$A$\leq$ 22 (green), 23 $\leq$A$\leq$ 38 (cyan), A $\geq$ 39
(blue)].

##  4 UHECRs and particle physics

###  4.1 Hadronic interactions

The interpretation of the experimental observables in terms of primary
composition is prone to large systematic uncertainties, mainly due to the lack
of knowledge on hadronic interactions at ultra-high energies.
Measuring the attenuation length of proton-induced air showers (those mostly
contributing to the tail of the $X_{\mathrm{max}}$ distribution), we extended
the measurement of the proton-air cross section up to 57 TeV [28, 29]. The
measure of the muonic component at the ground is more sensitive to the details
of the hadronic interactions along many steps of the cascade, like the
multiplicity of the secondaries and the fraction of electromagnetic component
with respect to the total signal. On the contrary, the intrinsic muon
fluctuations mostly depend on the first interaction.
Clear evidence for a deficit in the number of muons predicted by the models
was reported by our Collaboration by exploiting inclined showers detected by
the SD1500, where the electromagnetic component has been fully absorbed by the
atmosphere [30] (see left panel of Fig.5) and hybrid events [31]. Direct
measurements at lower energies confirmed the result [32]. The problems in the
modelling of hadronic interactions at ultra-high energy have further been
proven by studying the muon production depth [33] and the time profiles of the
signals recorded with SD [34].
The measured shower-to-shower fluctuations of the muonic component are shown
in the right panel of Fig.5 [35]; the expectations for their relative values
are compatible with the experimental results, while a significant discrepancy
between models and data appears in the average muon scale. Note that
systematic errors (square brackets) are dominant for $\left<R_{\mu}\right>$,
while for the muon fluctuations the uncertainties are mainly statistical
(error bars).
A new method to explore simple ad-hoc adjustments to the predictions of
hadronic interaction models was recently developed, exploiting the correlation
of $X_{\mathrm{max}}$ with the signal at ground level as a function of the
zenith angle [36]. The uncertainties on the $X_{\mathrm{max}}$ scale among
different hadronic interaction models are of the order of 30
$\mathrm{g~{}cm^{-2}}$; on the other hand, the evaluation of the primary mass
based on the signals at ground are prone to the inability of models to
describe the muonic component. This method allows us to use both the
longitudinal and lateral development of air showers to evaluate the flawns in
the models. An overall improvement is reached by shifting the expected depths
of shower maximum deeper in atmosphere; the hadronic signal at 1000 m from the
core should in this case be increased by a factor from 15 to 25%, smaller than
the one obtained with the previous methods. This preliminary result will be
further investigated by considering less simple adjustments to predictions, by
taking into account also the fluctuations in both $X_{\mathrm{max}}$ and
signal at ground.

Figure 5: Left: average logarithmic muon content $<lnR_{\mu}>$ as a function
of the average shower depth $\left<X_{\mathrm{max}}\right>$. Right: Measured
relative fluctuations in the number of muons as a function of the energy. The
grey band represents the expectations from the measured mass composition
interpreted with the interaction models.

###  4.2 Beyond the standard model

Various quantum gravity theories suggest that the Lorentz invariance may be
violated when approaching the Planck scale, resulting in a perturbative
modification of the particle dispersion relation
$E_{i}^{2}=p_{i}^{2}+m_{i}^{2}+\sum_{n=0}^{\infty}\delta^{(n)}_{i}E_{i}^{n+2}$
(where $i$=type of particle, $n$=approximation order).
The UHECRs collected at the Pierre Auger Observatory, coming from
extragalactic sources, represent the ideal candidates to search for these
effects, as these are expected to be suppressed at low energies and for short
travel distances. In presence of Lorentz invariance violation (LIV), the
attenuation length of photo-meson production or photo-disintegration may
become extremely large and suppress particle interaction during propagation in
the extragalactic space. Introducing these modifications in either the hadron
or the photon propagation, constraints on the level of violation can be set by
comparing via a best fit procedure the obtained flux and composition to those
measured at the top of the atmosphere [37]. While in the case of the hadron
sector the sensitivity to possible violations is very low, in the photon
sector we obtained constraining limits of $\delta^{(1)}_{\gamma}\sim 10^{-40}$
eV-1 and $\delta^{(2)}_{\gamma}\sim 10^{-60}$ eV-1 for $n$=1 and $n$=2
respectively.
Depending on the strength of the violation, the high energy available in the
collision of cosmic rays with the atmosphere can futhermore lead to
modifications of the shower development with respect to the standard
invariance case. Indeed, due to the much increased probability of interaction
for neutral pions, the hadronic component of the cascade is enhanced, leading
in particular to a larger production of muons. In this case, limits on the LIV
parameter can be derived comparing the observed strong decrease of the
relative fluctuations with the muon fluctuation measurement[38].

Super-heavy relics produced in the early Universe could decay into UHE
photons; detecting a flux of astrophysical photons with energies in excess of
$\sim 10^{17}$ eV would thus be a proof of the existence of super-heavy dark
matter. The non observation of these photons in the Pierre Auger Observatory
allowed us to derive a bound on the reduced coupling constant of gauge
interactions in the dark sector: $\alpha_{eff}\lesssim 0.09$ for
$10^{10}<M_{X}/GeV<10^{16}$ [40].

##  5 Conclusion

The Pierre Auger Observatory is providing the best and most precise data
world-wide on ultra-high energy cosmic rays, and a global view emerges from
the analyses of collected data. The suppression of the UHECR flux, measured
with extremely high statistical significance, appears to be due to source
exhaustion, although a contribution from propagation effects must be also
present. This result is supported by the increasingly heavy nuclear
composition in this energy region, and by the non observation of cosmogenic
neutrinos and photons, a copious production of which should be expected for a
proton dominated composition. The transition from Galactic to extra-galactic
cosmic rays can be placed at the second knee. This is supported by the
measured composition, which becomes lighter above the 2nd knee up to $\sim
2\times 10^{18}$ eV, and by the smooth transition from isotropy to a dipolar
anisotropy above 8 EeV, where the phase of the dipole points to the anti-
center of our Galaxy. Based on the results of the combined fit extended to the
region below the ankle, two extragalactic populations, the low energy one made
by proton/light elements, the high energy one by a mixture of nuclei, can
describe well the transition region. If the extragalactic component below the
ankle is composed only by protons, a secondary Galactic component of
intermediate nuclei (CNO) must be present. The new instep feature can be
explained as due to the interplay of He and CNO nuclei. The presence of heavy
nuclei at the highest energies implies a large effect of magnetic fields, the
study of which becomes more and more mandatory. It is still to be clarified if
the onset of the different components is due to propagation effects ($\propto
A$) or to Peters cycle ($\propto Z$). Correlations have been found between the
arrival directions of UHECRs and the direction of starburst Galaxies and AGN.
About 3-4 years more data are needed to reach the $5\sigma$ significance to
claim observation. Disentangling the two types of sources remains an issue, as
some of them (e.g. NGC4945 or NGC1068) are indeed composite AGN/SBG objects.
The strong limits on cosmogenic neutrinos and photons constrain models
predicting proton sources; top-down scenarios are excluded. The deficit in the
muon component predicted by models with respect to data is confirmed by many
different analyses of independent data sets. The recent measurement of the
fluctuations in the muonic component show that models describe reasonably well
the energy partition in the first interaction up to the highest energies,
while the discrepancy in the overall muon number should be explained by
changes of different hadronic interaction characteristics along all the shower
stages. The data collected by the Pierre Auger Observatory allow us to explore
the high energy frontier even beyond the standard model of particles.
Constraining limits to a possible violation of the Lorentz invariance and on
the mass and lifetime of super-heavy dark matter particles have been obtained.
New data on composition in the suppression region and more observables, like
those that will be provided by the ongoing upgrade of the observatory,
AugerPrime [41], will allow us to face the challenges these findings open to
us for the future.

## References

  * [1] A. Aab et al. (Pierre Auger Coll.), The Pierre Auger Cosmic Ray Observatory, Nucl. Instrum. Meth. A 798, 172 (2015), 10.1016/j.nima.2015.06.058.
  * [2] B. Dawson (Pierre Auger Coll.), The energy scale of the Pierre Auger Observatory,PoS(ICRC2019) 231.
  * [3] A. Aab et al. (Pierre Auger Coll.), Observation of a large-scale anisotropy in the arrival directions of cosmic rays above $\it{8\times 10^{18}}$ eV, Science 357, 1266 (2017), 10.1126/science.aan4338
  * [4] A. Aab et al. (Pierre Auger Coll.), Large-scale cosmic-ray anisotropies above 4 EeV measured by the Pierre Auger Observatory, Astrophys.J. 868, 4 (2018), 10.3847/1538-4357/aae689
  * [5] A. Aab et al. (Pierre Auger Coll.), Cosmic-Ray Anisotropies in Right Ascension Measured by the Pierre Auger Observatory, Astrop.J. 891, 142 (2020), 10.3847/1538-4357/ab7236
  * [6] D.Harari et al., Anisotropies of UHECR nuclei diffusing from extragalactic sources, Phys.Rev.D 92, 063014 (2015), 10.1103/PhysRevD.92.063014
  * [7] C.Ding et al.,The Imprint of Large-scale Structure on the UHECR Sky, Astrop.J. 913, L13 (2021), 10.3847/2041-8213/abf11e
  * [8] P. Abreu et al. (Pierre Auger Coll.), Arrival Directions of Cosmic Rays above $32$ EeV from Phase One of the Pierre Auger Observatory, Astrop.J. 935, 170 (2022), 10.3847/1538-4357/ac7d4e
  * [9] P. Tinyakov (Pierre Auger and Telescope Array Coll.), The UHECR dipole and quadrupole in the latest data from the original Auger and TA surface detectors, PoS(ICRC2021) 375.
  * [10] A. di Matteo (Pierre Auger and Telescope Array Coll.), UHECR arrival directions in the latest data from the original Auger and TA surface detectors and nearby galaxies, PoS(ICRC2021) 308.
  * [11] P. Abreu et al. (Pierre Auger Coll.), The energy spectrum of cosmic rays beyond the turn-down around $\it{10^{17}}$ eV as measured with the surface detector of the Pierre Auger Observatory, Eur.Phys.J. C 81, 966 (2021), 10.1140/epjc/s10052-021-09700-w
  * [12] A. Aab et al. (Pierre Auger Coll.), Depth of maximum of air-shower profiles at the Pierre Auger Observatory. II. Composition implications, Phys.Rev.D 90, 122006 (2014), 10.1103/PhysRevD.90.122006
  * [13] A. Aab et al. (Pierre Auger Coll.), Inferences on mass composition and tests of hadronic interactions from 0.3 to 100 EeV using the water-Cherenkov detectors of the Pierre Auger Observatory, Phys.Rev.D 96, 122003 (2017), 10.1103/PhysRevD.96.122003
  * [14] C. J. Todeiro-Peixoto (Pierre Auger Coll.), Estimating the Depth of Shower Maximum using the Surface Detectors of the Pierre Auger Observatory, PoS(ICRC2021) 340.
  * [15] B. Pont (Pierre Auger Coll.), The depth of the shower maximum of air showers measured with AERA, PoS(ICRC2021) 387.
  * [16] E. Mayotte (Pierre Auger Coll.), Indication of a mass-dependent anisotropy above $10^{18.7}$ eV in the hybrid data of the Pierre Auger Observatory, PoS(ICRC2021) 321
  * [17] Y. Tsunesada (Pierre Auger and Telescope Array Coll.), Joint analysis of the energy spectrum of ultra-high-energy cosmic rays measured at the Pierre Auger Observatory and the Telescope Array, PoS(ICRC2021) 337.
  * [18] A. Yushkov (Pierre Auger Telescope Array Coll.), Depth of maximum of air-shower profiles: testing the compatibility of measurements performed at the Pierre Auger Observatory and the Telescope Array experiment EPJ Web of Conf. 210 010009 (2018), 0.1051/epjconf/201921001009.
  * [19] A. Aab et al. (Pierre Auger Coll.), Probing the origin of ultra-high-energy cosmic rays with neutrinos in the EeV energy range using the Pierre Auger Observatory, JCAP 10, 022 (2019), 10.1088/1475-7516/2019/10/022
  * [20] P. Abreu et al. (Pierre Auger Coll.), A Search for Photons with Energies Above $\it{2\times 10^{17}}$ eV Using Hybrid Data from the Low-Energy Extensions of the Pierre Auger Observatory, Astrop.J. 933, 125 (2022), 10.3847/1538-4357/ac7393
  * [21] A. Aab et al. (Pierre Auger Coll.), Limits on point-like sources of ultra-high-energy neutrinos with the Pierre Auger Observatory, JCAP 11 004 (2019) 10.1088/1475-7516/2019/11/004.
  * [22] A. Aab et al. (Pierre Auger Coll.), A Search for Ultra-high-energy Neutrinos from TXS 0506+056 Using the Pierre Auger Observatory, Astrop.J. 902 105 (2020), 10.3847/1538-4357/abb476.
  * [23] M. Schimp (Pierre Auger Coll.), Combined Search for UHE Neutrinos from Binary Black Hole Mergers with the Pierre Auger Observatory, PoS(ICRC2021) 968.
  * [24] P. Ruehl (Pierre Auger Coll.), Follow-up Search for UHE Photons from Gravitational Wave Sources with the Pierre Auger Observatory, PoS(ICRC2021) 973.
  * [25] A. Aab et al. (Pierre Auger Coll.), Combined fit of spectrum and composition data as measured by the Pierre Auger Observatory, JCAP 04, 038, (2017), 10.1088/1475-7516/2017/04/038
  * [26] E. Guido (Pierre Auger Coll.), Combined fit of the energy spectrum and mass composition across the ankle with the data measured at the Pierre Auger Observatory, PoS(ICRC2021) 311.
  * [27] A. Aab et al. (Pierre Auger Coll.), Features of the Energy Spectrum of Cosmic Rays above $2.5\times 10^{18}$ eV Using the Pierre Auger Observatory, Phys.Rev.Lett. 125, 121106 (2020), 10.1103/PhysRevLett.125.121106
  * [28] P. Abreu et al. (Pierre Auger Coll.), Measurement of the Proton-Air Cross Section at $\sqrt{s}=57$ TeV with the Pierre Auger Observatory, Phys.rev.Lett. 109, 062002 (2012), 10.1103/PhysRevLett.109.062002 .
  * [29] R.Ulrich (Pierre Auger Coll.), Extension of the measurement of the proton-air cross section with the Pierre Auger Observatory, PoS(ICRC2015) 401.
  * [30] A. Aab et al. (Pierre Auger Coll.), Muons in air showers at the Pierre Auger Observatory: Mean number in highly inclined events, Phys.Rev.D 91, 032003 (2015), 10.1103/PhysRevD.91.032003
  * [31] A. Aab et al. (Pierre Auger Coll.), Testing Hadronic Interactions at Ultrahigh Energies with Air Showers Measured by the Pierre Auger Observatory, Phys. Rev. Lett. 117 192001 (2016) 10.1103/PhysRevLett.117.192001.
  * [32] A. Aab et al. (Pierre Auger Coll.), Direct measurement of the muonic content of extensive air showers between $\it{2\times 10^{17}}$ eV and $\it{2\times 10^{18}}$ eV at the Pierre Auger Observatory, Eur.Phys.J.C 80, 751 (2020), 10.1140/epjc/s10052-020-8055-y
  * [33] A. Aab et al. (Pierre Auger Coll.), Muons in air showers at the Pierre Auger Observatory: Measurement of atmospheric production depth, Phys. Rev. D 90 (2014) 012012; D92 (2015) 019903(E) 10.1103/PhysRevD.92.019903.
  * [34] A. Aab et al. (Pierre Auger Coll.), Inferences on mass composition and tests of hadronic interactions from 0.3 to 100 EeV using the water-Cherenkov detectors of the Pierre Auger Observatory, Phys. Rev. D 96 122003 (2017) 10.1103/PhysRevD.96.122003.
  * [35] A. Aab et al. (Pierre Auger Coll.), Measurement of the Fluctuations in the Number of Muons in Extensive Air Showers with the Pierre Auger Observatory, Phys.Rev.Lett. 126, 152002 (2021), 10.1103/PhysRevLett.126.152002
  * [36] J. Vicha (Pierre Auger Coll.), Adjustments to Model Predictions of Depth of Shower Maximum and Signals at Ground Level using Hybrid Events of the Pierre Auger Observatory, PoS(ICRC2021) 310.
  * [37] P. Abreu et al. (Pierre Auger Coll.), Testing effects of Lorentz invariance violation in the propagation of astroparticles with the Pierre Auger Observatory, JCAP 01, 023 (2022), 10.1088/1475-7516/2022/01/023
  * [38] C. Trimarelli (Pierre Auger Coll.), Constraining Lorentz Invariance Violation using the muon content of extensive air showers measured at the Pierre Auger Observatory, PoS(ICRC2021) 340.
  * [39] P. Abreu et al. (Pierre Auger Coll.), Limits to gauge coupling in the dark sector set by the non-observation of instanton-induced decay of Super-Heavy Dark Matter in the Pierre Auger Observatory data, arXiV:2203.08854, 10.48550/arXiv.2203.08854
  * [40] O.Deligny (Pierre Auger Coll.),Limits to gauge coupling in the dark sector of super-heavy dark matter particles from the Pierre Auger Observatory data, this conference.
  * [41] A. Castellina (Pierre Auger Coll.), AugerPrime: the Pierre Auger Observatory Upgrade, EPJ Web of Conf. 210, 06002 (2019), 10.1051/epjconf/201921006002


# Geometrical bounds on irreversibility in squeezed thermal bath

Chen-Juan Zou1 Yue Li2 Jia-Bin You3 Qiong Chen1 Wan-Li Yang2 Mang Feng2 1Key
Laboratory of Low-Dimensional Quantum Structures and Quantum Control of
Ministry of Education, Key Laboratory for Matter Microstructure and Function
of Hunan Province, Department of Physics, Hunan Normal University, Changsha
410081, China  2State Key Laboratory of Magnetic Resonance and Atomic and
Molecular Physics, Innovation Academy for Precision Measurement Science and
Technology, Chinese Academy of Sciences, Wuhan 430071, China  3Institute of
High Performance Computing, Agency for Science, Technology, and Research
(A*STAR), 1 Fusionopolis Way, #16-16 Connexis, Singapore 138632.

###### Abstract

Irreversible entropy production (IEP) plays an important role in quantum
thermodynamic processes. Here we investigate the geometrical bounds of IEP in
nonequilibrium thermodynamics by exemplifying a system coupled to a squeezed
thermal bath subject to dissipation and dephasing, respectively. We find that
the geometrical bounds of the IEP always shift in contrary way under
dissipation and dephasing, where the lower and upper bounds turning to be
tighter occurs in the situation of dephasing and dissipation, respectively.
However, either under dissipation or under dephasing, we may reduce both the
critical time of the IEP itself and the critical time of the bounds for
reaching an equilibrium by harvesting the benefits of squeezing effects, in
which the values of the IEP, quantifying the degree of thermodynamic
irreversibility, also becomes smaller. Therefore, due to the nonequilibrium
nature of the squeezed thermal bath, the system-bath interaction energy brings
prominent impact on the IEP, leading to tightness of its bounds. Our results
are not contradictory with the second law of thermodynamics by involving
squeezing of the bath as an available resource, which can improve the
performance of quantum thermodynamic devices.

Keywords: squeezed thermal bath, irreversible entropy production, geometrical
bounds

##  I Introduction

Over the past two decades, the nonequilibrium phenomena and __ thermodynamic
irreversibility quantified by irreversible entropy production (IEP) have drawn
much attention, since this fundamental concept is one of the cornerstones of
classical and quantum thermodynamics TM1 ; TM2 ; TM3 ; TM4 ; TM5 ; TM6 ; TM7 ;
TM8 ; TM9 ; TM10 . As is well known, the positivity of entropy production has
been universally captured by the conventional second law of thermodynamics
(SLT) TM9 , which quantitatively characterize the interplay between the
exchange of energy and the irreversibility by introducing the state function
entropy fun1 ; fun2 . Also, through specifying a lower bound for the
irreversible entropy change, the related Clausius inequality provides a
fundamental feature of irreversible phenomena. Noted that this lower bound
(zero) is trivially independent of how far from equilibrium a process operates
TM6 .

Recently, with restriction to a specific class of nonequilibrium phenomena,
such as thermal relaxation process, rich features of thermodynamic
irreversibility have been found successively TM2 . Such a lower bound for
classical, near-equilibrium transformation processes has been derived by mean
of a geometric approach rmp1 . S. Deffner and E. Lutz treated the system’s
Hilbert space as a Riemannian manifold and extended the classical case rmp1
to the nonequilibrium closed quantum system. The obtained generalized Clausius
inequality states that the thermodynamic irreversibility is bounded in terms
of the Bures length between the final state and the corresponding equilibrium
state, by using information geometry TM6 _._ Soon afterwards, they broadened
the closed quantum system further __ to the weakly coupled open quantum
system, then obtained the exact microscopic expressions for the nonequilibrium
entropy production TM7 . Along this direction, reference fun1  theoretically
and experimentally determined a sharper geometrical bound for a qubit
thermalization process and obtained a tighter version of the Clausius
inequality following a similar approach.

In many cases of interest, however, having a sharper lower bound is essential.
Therefore, theses above-mentioned progress hence stimulated successive studies
on the tightness of the geometrical bounds on irreversibility in open quantum
system TM6 ; TM7 ; TM8 ; TM9 ; TM10 ; fun1 . These recent publications tried
to develop theories to further understand the thermodynamic irreversibility
inherent to nonequilibrium processes. In particular, based on the variational
principle and time-reversed map, the authors in reference TM10  obtained a
information-theoretical bound for entropy production in a relaxation process
by a geometrical distance on the Riemannian manifold Rie , which was
experimentally validated by a single ultracold trapped ion 40Ca Ca .

T. V. Vu and Y. Hasegawa strengthened the Clausius inequality and proved that
IEP is bounded from below by a modified Wasserstein distance (quantum
generalization of the Wasserstein metric) between the initial and final states
TM8 .__ Thereafter, they extended this single-bath case to the case of
multiple-bath, and refined the SLT in a quantum regime, through deriving the
fundamental __ bound on irreversibility for thermal relaxation processes of
Markovian open quantum systems TM9 . __

On the other hand, quantum bath engineering techniques are powerful tools that
enable the realization of arbitrary thermal and nonthermal environments no1 ;
no2 ; no3 ; no4 ; no5 ; no6 ; no7 ; no8 . Additionally, due to the unique
vantage of the quantum control, bath engineering has aroused widespread
interest in the context of quantum thermodynamical processes. Various
strategies have attempted to improve the the performance of thermodynamic
devices dev1 ; dev2 ; dev3 ; dev5 ; dev6 ; dev7 ; dev8 ; dev9 ; dev10 , whose
efficiency is usually reduced by the presence of IEP __TM7 . For instance, the
use of a squeezed thermal bath allows us to operate thermodynamic devices
beyond the classical bound no8 ; dev1 . In particular, the experiment in
reference dev8  demonstrated that the efficiency of the quantum heat engine
may go beyond the standard Carnot efficiency by employing a squeezed thermal
bath.

Figure 1: Schematic diagram for two-level system (excited state
$\left|e\right\rangle$ and ground state $\left|g\right\rangle$) interacting
with a squeezed thermal bath at temperature $T$ with squeezing parameter $s$.
The transition frequency between the two levels is $\omega_{0}$. The evolution
induced by the system-bath interaction produces irreversible entropy.

However, the influence of squeezed thermal bath on the bounds of IEP in the
quantum thermodynamics has been largely unexplored. A precise characterization
of the IEP in such an unchartered domain, and a general framework providing a
deeper understanding of the associated quantum thermodynamic phenomena
therefore appear necessary. Therefore, it is instructive to look into the role
of squeezed thermal bath playing during the process of thermodynamic
irreversibility.

In this work, we study and quantify the geometrical bounds on irreversibility
of a quantum system contacting with a squeezed thermal bath. Here we consider
two distinct situations including the dissipation model and dephasing model,
respectively. Starting with Born-Markovian quantum master equation, we derived
analytical expressions for time-dependent reduced density matrix of system, in
which the squeezing parameter of bath is involved. Then through quantifying
the degree of irreversibility by the IEP, we find that the geometrical bounds
of IEP decrease (increase) with the growth of degree of squeezing in the case
of dissipation (dephasing) model, respectively. Additionally, the common
feature is that the critical times of the IEP itself and that of IEP’s bounds
reaching equilibrium, as well as the values of IEP quantifying the degree of
thermodynamic irreversibility are reduced due to the presence of squeezing
effect, for these two models. Therefore, due to the nonequilibrium nature of
squeezed thermal bath, the interaction energy between the system and bath
brought important impact on the irreversibility, as well as the tightness of
its bounds. As expected, our finding obeys the principles of thermodynamics
and reveals richer features of thermodynamics relaxation process, and the
presence of a quantum property, such as the squeezing effects included in the
bath, could be used to serve as an available resource to improve the
performance of the quantum thermodynamic devices.

The paper is organized as follows. We give a brief account of the method about
geometrical bounds on irreversibility in open quantum system (Section II). We
evaluate the geometrical bounds on irreversibility for squeezed thermal bath
in the dissipation model (Section III) and dephasing model (Section IV),
respectively, and we conclude and give prospects for future developments in
Section V.

##  II Geometrical bounds of irreversible entropy production

Considering an arbitrary quantum system with the Hamiltonian coupled to a
thermal bath. Usually, the quantum system is initialized in a given state
$\rho_{0}$, then interacts with a bath at temperature $T$. The evolution
induced by the interactions brings the system in a state $\rho(t)$ and
produces irreversible entropy. Furthermore, the system will thermalize with
the bath and then asymptotically reaching the unique canonical equilibrium
state $\rho_{th}$ if the system Hamiltonian $H$ maintains constant. The total
entropy variation of the system is defined as

$\Delta S_{tot}=\Delta S_{ir}+\Delta
S_{re}=S\left[\rho(t)\right]-S\left[\rho_{0}\right],$ (1)

where $S\left[\rho\right]=-tr\left(\rho\ln\rho\right)$ is the von Neumann
entropy. The IEP (irreversible part of the total entropy variation $\Delta
S_{tot}$) is denoted by

$\Delta
S_{ir}(t)=S\left(\rho_{0}\|\rho_{th}\right)-S\left(\rho(t)\|\rho_{th}\right),$
(2)

which is the thermodynamic irreversibility under consideration in the present
work TM2 ; fun1 .
$S\left(\rho_{1}\|\rho_{2}\right)=tr\left(\rho_{1}\ln\rho_{1}\right)-tr\left(\rho_{1}\ln\rho_{2}\right)$
represents the quantum relative entropy of $\rho_{1}$ to $\rho_{2}$. The
entropy flow between the system and the environment (reversible part of the
total entropy variation $\Delta S_{tot}$) is $\Delta S_{re}=\Delta Q/T$ ,
where $\Delta Q=tr\left(H\rho(t)\right)-tr\left(H\rho_{0}\right)$ is the heat
absorbed by the system fun1 . On the other hand, the Clausius inequality
$\Delta S_{ir}\geq 0$ putting forward the lower limit of IEP is always non-
negative. In order to deepen the understanding of how much energy in the
irreversible process is consumed, it is essential to search a sharper or
tighter bound on irreversibility.__ By treating the Hilbert space of the
system as a Riemannian manifold, the relationship between IEP and the geodesic
distance $D$ corresponding to the metric that is contractive under complete
positive and trace preserving maps can be directly established, and a
generalized form of Clausius inequality can be obtained by deriving the
Wigner-Yanase length between the initial and final states of the system. Note
that the only cases in which an analytical expression for the geodesic
distance is known are the Wigner-Yanase metric
$D_{WY}\left(\rho_{1},\rho_{2}\right)=\arccos\left[tr\left\\{\sqrt{\rho_{1}}\sqrt{\rho_{2}}\right\\}\right]$
and the quantum Fisher information metric
$D_{QF}\left(\rho_{1},\rho_{2}\right)=\arccos\left[tr\left\\{\sqrt{\sqrt{\rho_{1}}\rho_{2}\sqrt{\rho_{1}}}\right\\}\right]$
fun1 . Based on these analytical expressions, one can obtain the __
geometrical LB of IEP as

$\Delta
S_{ir}\left(t\right)\geq\frac{8}{\pi^{2}}\max_{\left\\{X=QF,WY\right\\}}D_{X}^{2}\left(\rho_{0},\rho(t)\right),$
(3)

and the __ geometrical UB is

$\Delta S_{ir}\left(t\right)\leq
S\left(\rho_{0}\|\rho_{th}\right)-\frac{8}{\pi^{2}}\max_{\left\\{X=QF,WY\right\\}}D_{X}^{2}\left(\rho(t),\rho_{th}\right).$
(4)

Figure 2: (a) Time dependence of LB and UB of IEP under different degree of
squeezing $s$ in the dissipation model. Time dependence of IEP (blue solid
line) and its LB (red dotted line) and UB (green dashed line) in the case of
(b) $s=0$ and (c) $s=2$. The initial state of system is
$\left(\left|e\right\rangle+\left|g\right\rangle\right)/\sqrt{2}$. Hereafter,
we choose $\omega_{0}$ as the frequency unit, $T=0.34$, and $\phi=0.$

From the above two relations one can define the related bound gap as

$\Delta U=UB-LB,$ (5)

and the deviation of the IEP from the LB or UB is given by

$\Delta\delta_{L}=\Delta S_{ir}\left(t\right)-LB,$ (6)

or

$\Delta\delta_{U}=UB-\Delta S_{ir}\left(t\right),$ (7)

respectively. In this paper, we say a LB (UB) is relative tighter if the LB
(UB) takes a larger (smaller) value compared with the case of conventional
thermal bath.

##  III Geometrical bounds on irreversibility in the dissipation model

Here we consider the dissipation model taking into account the effect of the
squeezed thermal bath with temperature $T$, in the case of single-excitation.
The total Hamiltonian is (in units of $\hbar=1$)

$\displaystyle H$ $\displaystyle=$ $\displaystyle H_{S}+H_{B}+H_{SB}$ (8)
$\displaystyle=$
$\displaystyle\omega_{0}\hat{\sigma}_{+}\hat{\sigma}_{-}+\sum\limits_{k}\omega_{k}\hat{b}_{k}^{\dagger}\hat{b}_{k}+\sum\limits_{k}(g_{k}\hat{\sigma}_{+}\hat{b}_{k}+H.c.),$

where $H_{S}$, $H_{B}$, and $H_{SB}$ stand for the Hamiltonians of the system,
bath, and system-bath interaction, respectively.
$\hat{\sigma}_{+}(\hat{\sigma}_{-})=\left|e\right\rangle\left\langle
g\right|(\left|g\right\rangle\left\langle e\right|)$ and $\omega_{0}$ are the
inversion operator and transition frequency of the system with
$\left|e\right\rangle$ and $\left|g\right\rangle$ being the excited and ground
states. $\hat{b}_{k}^{\dagger}(\hat{b}_{k})$ are the creation (annihilation)
operators of the k-th mode of the bath. The coupling strength between the
system and the bath is denoted by $g_{k}$.

The master equation, in the interaction picture, is given by the following
Lindblad form master1 ; master2

$\displaystyle\dot{\rho}_{s}(t)$ $\displaystyle=$ $\displaystyle\frac{\gamma
N}{2}D[\sigma_{+}]\rho_{s}(t)+\frac{\gamma\left(1+N\right)}{2}[\sigma_{-}]\rho_{s}(t)$
(9) $\displaystyle-\gamma M\sigma_{+}\rho_{s}(t)\sigma_{+}-\gamma
M^{\ast}\sigma_{-}\rho_{s}(t)\sigma_{-},$

where $\dot{\rho}_{s}(t)=d\rho_{s}(t)/dt$, and $D[A]\rho=2A\rho
A^{+}-A^{+}A\rho-\rho A^{+}A$. The spontaneous emission __ rate of system is
$\gamma$ and
$N=N_{th}\left[\cosh^{2}\left(s\right)+\sinh^{2}\left(s\right)\right]+\sinh^{2}\left(s\right)$,
and $M=-\sinh\left(2s\right)e^{i\phi}\left(2N_{th}+1\right)/2$ with $s$ and
$\phi$ the bath squeezing strength and phase, respectively.
$N_{th}=1/\left(e^{\omega_{0}\beta}-1\right)$ is the Plank distribution, where
$\beta=1/k_{B}T$ with $k_{B}=1$ the Boltzmann constant.

Figure 3: (a) The values of the IEP $\triangle S_{irr}$ in the parameter plane
of $\\{s,t\\}$ in the dissipation model. (b) Time dependence of IEP $\triangle
S_{irr}$ under the different degrees of squeezing $s$.

Rewriting the density matrix as $\rho_{s}(t)$
$=\left(\mathbf{I}+\vec{r}\left(t\right)\cdot\vec{\sigma}\right)/2$ by mean of
the Bloch vector $\vec{r}\left(t\right)=Tr[\vec{\sigma}\rho_{s}(t)]$ with the
Identity matrix $\mathbf{I}$, we can transform the master equation (9) into
the Bloch equation

$\frac{d}{dt}\vec{r}\left(t\right)=\xi\vec{r}\left(t\right)+\vec{m},$ (10)

with

$\xi=\left(\begin{array}[]{ccc}-\frac{\tilde{\gamma}+2\gamma M}{2}&0&0\\\
0&-\frac{\tilde{\gamma}-2\gamma M}{2}&0\\\
0&0&-\tilde{\gamma}\end{array}\right)$ (11)

and $\vec{m}=(0,0,,-\gamma_{0})^{T}$. Here
$\tilde{\gamma}=\gamma\left(2N+1\right)$ is the total transition rate.
Assuming the system to be initially in the ground state
$\left|g\right\rangle$, a straightforward calculation yields the analytical
solution

$\rho_{s}(t)=\left(\begin{array}[]{cc}\frac{1-\nu+(1-\nu)\left\langle\sigma_{z}\right\rangle_{ss}}{2}&\mu\left\langle\sigma_{-}\right\rangle_{ss}\\\
(1-e^{\frac{-\gamma_{s}t}{2}})\left\langle\sigma_{+}\right\rangle_{ss}&\frac{1+\nu-(1-\nu)\left\langle\sigma_{z}\right\rangle_{ss}}{2}\end{array}\right),$
(12)

Here $\left\langle\sigma_{\pm}\right\rangle_{ss}$ and
$\left\langle\sigma_{z}\right\rangle_{ss}\ $are the stationary solutions of
differential equation (10).
$\mu=\\{\gamma_{s}-e^{\left(-4\tilde{\gamma}+\gamma_{s}\right)t/4}[\gamma_{s}\cos(\gamma_{s}t/4)+(\gamma_{s}+\gamma
M)\sin(\gamma_{s}t/4)]\\}/\gamma_{s}$ and
$\nu=\tilde{\gamma}e^{\left(-4\tilde{\gamma}+\gamma_{s}\right)t/4}[\cos(\gamma_{s}t/4)-\sin(\gamma_{s}t/4)]/\gamma_{s}$,
where $\gamma_{s}=\tilde{\gamma}+2\gamma M$.

The time dependence of the IEP, and its geometrical bounds are explored by
numerically calculating the quantities (equations (2) and (3), (4)) and
plotted in figure 2. It is clear from figure 2 that IEP and its bounds (UB and
LB) increase monotonically with time toward the corresponding equilibrium
values for a squeezed thermal bath. For the non-squeezing case ($s=0$), our
results based on quantum master equation are fully consistent with that for
the thermal bath in reference fun1  using the method of Kraus operators. In
the case of squeezed thermal bath ($s>0$), with increase in degree of
squeezing the UB becomes higher in the early stage of evolution $t\lesssim
0.15$, and then becomes reduced with the growth of $s$ when $t\gtrsim 0.45$,
compared with the thermal bath ($s=0$), as shown in figure 2(a). Also, the
dynamical behavior of LB shares similar features with that of UB. Another
common feature between UB and LB is that their equilibrium times decrease with
the growing degree of squeezing, for instance, the equilibrium times for UB
(LB) are 4.5 (7.5), 2 (1.4) and 0.15 (0.2) when $s=0$, $1$, and $2$,
respectively. We conclude that the summarized overall trends of geometrical
bounds (LB and UB) found in reference fun1  also hold for squeezed thermal
bath, while only the UB becomes tighter in the long-time limit, and the LB
exhibits subtle tightness in the early stage of evolution, compared with the
traditional thermal bath fun1 .

Let us next examine the time dependence of the IEP, the numerical simulation
plotted in figure 2(b), and 2(c) suggests that the values of IEP are well
limited in the region between the LB and UB, and the squeezed thermal bath has
prominent influences in both the concrete values of IEP and the critical times
$T_{c}$ of IEP reaching equilibrium, where the values of $T_{c}$ become less
for growing values of $s$. Another observation in figure 2(b) is that the
deviation $\Delta\delta_{L}$ (equation (6))increases monotonically and
gradually coincides with $\Delta U$ (equation (5)), whereas the deviation
$\Delta\delta_{U}$ (equation (7)) decreases monotonically and gradually
disappears to zero in the case of $s=0$. With the increase of squeezing
strength, i.e., $s=2$, as shown in figure 2(c), the deviation
$\Delta\delta_{L}$ ($\Delta\delta_{U}$) increases (decreases) gradually to a
fixed value 0.51 (0.542), and the rates of change of deviations
$\Delta\delta_{L}$ and $\Delta\delta_{U}$ become slower. From the above
results we deduce that the actual amount of IEP becomes depart from its UB
gradually and approaches its LB with the growth of $s$.

To get a clear picture of how the IEP evolves in the parameter space of
$\\{s,t\\}$, we plot the evolution of IEP in figure 3. It shows that the
values of IEP are obviously dependent on the squeezing parameters of bath.
Here we provide remarks on the parameter dependence. Although the values of
$s$ are not directly related to the system-bath interaction, whereas they
depend on both correlation time and occupation number of the bath and then
affect immensely the relaxation dynamics and steady state of the relevant
system during the thermodynamic process. As a result, any change of this key
parameter will bring significant influences on the irreversibility, and the
inherent squeezing effect stemming from the bath plays a crucial role in
understanding the relative tightness of the bounds. As we have derived
analytically in the previous paragraphs, the above summarized dependence of
bounds and IEP on the squeezing parameter is reflected in the equations (3, 4,
12).

Traditionally, the IEP could be used to evaluate the performance of
thermodynamic devices, such as the availability (or exergy) that can be
extracted from a given system, and the maximal useful work, which are usually
reduced by the presence of irreversibility TM7 ; book1 . It implies that the
irreversibility can be restrained by controlling the amount of IEP through
adjusting the degree of squeezing effect, as shown in figure 3(b), where the
values of IEP become less for growing values of $s$. From the perspective of
quantum bath engineering, employing a squeezed thermal bath is a promising
avenue of using squeezing effect as a quantum resource to improve the
performance thermodynamic devices no8 ; dev1 ; dev8 .

As we known, squeezing effect that is rooted in Heisenberg’s uncertainty
principle can be defined as the reduction in the uncertainty of some
observable, at the cost of the build-up in the conjugate one dev3 ; book2 ;
Natt . Physically, the squeezing involved in the bath thereby inevitably
modifies the nonunitary relaxation dynamics of system and the relevant
irreversibility during thermodynamic process. Compared with the thermal bath,
the squeezed thermal bath is taken out of thermodynamics equilibrium through
the squeezing operating, with a consequence that its excitation number
changing from $N_{th}$ to $N=N_{th}(\cosh^{2}r+\sinh^{2}r)+\sinh^{2}r$ dev3 ;
enn , which can be seen as an increase in its effective temperature
$T_{eff}=\omega_{h}/k_{B}\ln[1/(N_{th}^{-1}+1)]$ with a higher frequency
$\omega_{h}>\omega_{0}$. Therefore, being purely quantum mechanical fuel in
nature, squeezed thermal bath is beneficial in its own way by providing us
with more compact energy-storage and higher effective high-temperature bath
without being actually too hot TM5 . That is to say, the squeezed thermal
state has the same entropy as the Gibbs state, but increased mean energy,
which is instrumental in the suppression of irreversibility.

##  IV Geometrical bounds on irreversibility in the dephasing model

Next we focus on the dephasing model with respect to the squeezed thermal
bath, where the bath operator is simply a sum of linear couplings to the
coordinates of a continuum of harmonic oscillators described by a spectral
density function $J(\omega)$ open1 ; open2 ; open3 ; open4 ; open5 ; open6 ,
and the decay of the coherence occurs without a decay of the corresponding
populations. Now the total Hamiltonian is

$\displaystyle H$ $\displaystyle=$ $\displaystyle H_{S}+H_{B}+H_{SB}$ (13)
$\displaystyle=$
$\displaystyle\omega_{0}\hat{\sigma}_{+}\hat{\sigma}_{-}+\sum\limits_{k}\omega_{k}\hat{b}_{k}^{\dagger}\hat{b}_{k}+\hat{\sigma}_{z}\sum\limits_{k}(g_{k}\hat{b}_{k}+H.c.).$

The dynamics of the system can be characterized by the reduced density matrix
which is obtained by tracing out the degrees of freedom of the bath. In the
interaction picture, the reduced density matrix of system can be written as
open4 ; open5 ; open6 ; dep1 ; dep2

$\rho_{s}(t)=\left[\begin{array}[]{cc}\rho_{ee}&\rho_{eg}\Gamma\left(t\right)\\\
\rho_{ge}\Gamma^{\ast}\left(t\right)&\rho_{gg}\end{array}\right],$ (14)

where the phase decay behavior of the qubit under the influence of the bath is
denoted by the factor
$\Gamma\left(t\right)=Tr_{B}\rho_{B}\prod\nolimits_{k}\exp[\alpha_{k}(t)\hat{b}_{k}^{\dagger}-\alpha_{k}^{\ast}(t)\hat{b}_{k}]$,
where $\alpha_{k}(t)=2\frac{g_{k}}{\omega_{k}}\left(1-e^{i\omega_{k}t}\right)$
open4 . The associated master equation is given by

$\dot{\rho}_{s}(t)=\frac{-i\tilde{\epsilon}(t)}{2}[\sigma_{z},\rho_{s}(t)]+\frac{\tilde{D}(t)}{2}[\sigma_{z}\rho_{s}(t)\sigma_{z}-\rho_{s}(t)],$
(15)

where $\tilde{D}(t)=-\frac{d\ln(\left|\Gamma(t)\right|)}{dt}$and
$\tilde{\epsilon}(t)=-$Im$[\frac{d\Gamma(t)/dt}{\Gamma(t)}].$

In the following we consider that the bath starts from a squeezed thermal
state dep1 ; dep2 ; therm1 ; therm2

$\rho_{B}(0)=\hat{\varsigma}\rho_{th}\hat{\varsigma}^{\dagger},$ (16)

where $\rho_{th}=e^{-\beta H_{B}}/Z_{\beta}$ is thermal state with $Z_{\beta}$
the partition function. $\hat{\varsigma}$ $=\sum\nolimits_{k}\hat{s}_{k}$,
where
$\hat{s}_{k}=e^{[(s_{k}^{\ast}e^{-i\phi_{k}}\hat{b}_{k}^{2}-s_{k}e^{i\phi_{k}}(\hat{b}_{k}^{\dagger})^{2})/2]}$
is the squeezing operator for the boson bath mode $\hat{b}_{k}$ with $s_{k}$
and $\phi_{k}$ being the bath squeezing strength and phase, respectively. In
this situation, the function $\Gamma\left(t\right)$ could be evaluated under
the summation of the modes of the squeezed thermal bath as open4 ; open5 ;
open6 ; dep1 ; dep2

$\Gamma\left(t\right)=e^{[-\sum_{k}\frac{4\left|g_{k}\right|^{2}}{\omega_{k}^{2}}\left(1-\cos\omega_{k}t\right)\gamma_{k}\left(t\right)\coth(\frac{\omega_{k}}{2T})]}$
(17)

with $\gamma_{k}\left(t\right)=\cosh 2s_{k}-\sinh
2s_{k}\cos\left(\omega_{k}t-\Delta\phi_{k}\right)$ and $\Delta\phi_{k}$ is the
phase difference between the squeezing phase $\phi_{k}$ relative to the phase
of the coupling strength $g_{k}$.

Figure 4: (a) Time dependence of LB and UB of IEP under different degrees of
squeezing $s$ in the dephasing model. Time dependence of IEP (blue solid line)
and its LB (red dotted line) and UB (green dashed line) in the case of (b)
$s=0$ and (c) $s=1$. The initial state of system is
$\left(\left|e\right\rangle+\left|g\right\rangle\right)/\sqrt{2}$, and the
parameter $\Delta\phi=\pi/4$.

Substituting the coupling spectral density
$J(\omega)=2\pi\sum_{k}\left|g_{k}\right|^{2}\delta\left(\omega-\omega_{k}\right)$
into the equation (17), we can transform the above summation in
$\Gamma\left(t\right)$ into an integral for continuous bath modes as

$\displaystyle\Gamma\left(t\right)$ $\displaystyle=$
$\displaystyle\exp\\{-\int_{0}^{\infty}\frac{d\omega}{\pi\omega^{2}}2J(\omega)(1-\cos\omega
t)\coth(\omega/2T)$ (18)
$\displaystyle\times\left[\cosh(2s)-\sinh(2s)\cos\left(\omega
t-\Delta\phi\right)\right]\\}.$

Figure 5: (a) The values of the IEP $\triangle S_{irr}$ in the parameter plane
of $\\{s,t\\}$ in the dephasing model. (b) Time dependence of IEP $\triangle
S_{irr}$ under the different degrees of squeezing $s$. Here the parameter
$\Delta\phi=\pi/4$.

In the present work we adopt the Ohmic coupling spectral density
$J(\omega)=\eta\omega e^{-\omega/\varpi}$ with $\varpi$ cutoff frequency and
$\eta$ is the coupling strength noise . Note that such engineering of the
spectrum’s Ohmicity seems possible when simulating the dephasing model using
trapped ultracold atom, as demonstrated in reference atom1 ; atom2 . In the
high-temperature regime, the dephasing process is Markovian, and after
straightforward algebra, one finds

$\Gamma\left(t\right)=e^{\\{-\frac{2\eta Tt}{\pi}\left[\pi\cosh(2s)-\ln
4\sinh(2s)\sin\Delta\phi\right]\\}},$ (19)

where the approximation $\coth(\omega/2T)\approx 2T/\omega$ has been taken.
Expression (equation (19)) is the exact analytic for the time-dependent
dephasing rate $\Gamma\left(t\right)$ in the present model.

Special attention was paid to the time dependence of the IEP as well as its
geometrical bounds, as shown in figure 4, where one can find that all the
quantities exhibit asymptotic behaviors approaching their stationary values.
They correspond to relaxations of the system through the dephasing channel due
to the system-bath coupling. Regarding the time dependence of the bounds, we
encounter another common feature is that the equilibrium times of bounds
decrease with the growing degree of squeezing, for instance, the equilibrium
times for UB (LB) are 3 (7.3), 2.5 (3.3) and 0.4 (0.6) when $s=0$, $1$, and
$2$, respectively. But unlike the dissipation model, both UB and LB increase
monotonically as the squeezing character of the bath grows, and only the LB
becomes tighter (compared with the thermal bath) during the whole dynamic
process in the dephasing model. As a contrast, the tightness of LB only
appears in the early stage of evolution for the dissipation model. By
comparing the time evolutions of IEP with the two blue solid-lines in panel
(b,c) of figure 4, we find that the IEP is well located inside the region
between the LB and UB, and reaches its stationary value faster with the
increase of the squeezing parameter $s$, and the values of the IEP is reduced
due to the existence of squeezing effect.

In figure 5, we provide numerical estimates of the IEP in the parameter space
of $\\{s,t\\}$. The figure 5 tells us that one can precisely control the
thermodynamic irreversibility through adjusting the parameters of the bath. As
shown in figure 5(b), under the dephasing model, the IPE reaches equilibrium
faster as the squeezed parameter increases, and the value of the IEP in long-
time limit $\bigtriangleup S_{irr}(\infty)$ will eventually converge together,
irrespective of the values of parameter $s$. It means that, in the dephasing
model, the squeezing effect could not make too much impact on the
thermodynamic irreversibility in the long-time limit, although the existence
of squeezing drives the system into equilibrium faster. Physically, on a
fundamental level, quantum coherence and the related dephasing process could
also alter the possible state transitions in thermodynamic processes Los , and
may even modify the fluctuation–dissipation relation FD1 ; FD2  and quantum
nonequilibrium work relation smith . Additionally, when a system relaxes to
equilibrium through contacting with a thermal bath, quantum coherences are
known to contribute an additional term to the IEP Iep1 ; Iep2 ; Iep3 .
Different from the dissipation case where the system can exchange energy with
its bath, in the dephasing model, this open system can never exchange energy
with its bath. But the information and correlation exchange between the system
and the bath are dominant during the dynamics and this exchange also influence
the IEP. As a result, any alteration in the von Neumann entropy (basis of IEP
and its bound) resulting from the relaxation process (dissipation or
dephasing) has contributions not only from the change in population but also
from decoherence. In this regard, it was pointed out that the entropy
production can be split in two contributions, an incoherent one (stemming from
populations) and a coherent one (stemming from quantum coherences) ep1 ; ep2 ;
ep3 .

##  V Discussion and conclusions

The study of IEP is of importance due to its intimate relation with the arrow
of time in classical and quantum systems arr1 ; arr2 , the SLT law1 ; law2 ;
law3 , thermodynamic operations and thermal machines mac1 ; mac2 ; mac3 ; mac4
, and quantum and classical speed limits TM8 ; speed1 ; speed2 . Therfore,
tightening the bounds of IEP not only deepens our understanding of how much
entropy production changes during the thermodynamic process but also provides
insights into how to improve the performance of quantum thermodynamic devices.

Interaction with a squeezed thermal bath is not the only generalized process
that goes beyond the typical settings in classical thermodynamics. Our
findings demonstrate how to utilize the squeezing effect of bath as a resource
to control the irreversibility, where the use of nonthermal bath offers more
degrees of control and manipulation, such as the amount of squeezing. Note
that quantum bath engineering techniques has become powerful tools that enable
the realization of arbitrary thermal and nonthermal bath, for instance,__
experimental realizations of squeezed thermal states range from
superconducting circuit QED mw1 ; mw2 ; mw3  to optomechanical mechanical
oscillators osc1 ; osc2 . The key parameters considered in our numerical
simulation, such as the inverse temperature $\beta$ and the degree of
squeezing $s$, could be experimentally controlled using the current
technologies demonstrated in the above-mentioned experiments. Additionally,
there have been many experiments focus on the assessment of nonequilibrium
thermodynamics irreversibility using the technology of quantum trajectories of
stochastic dynamics in nuclear magnetic resonance setups TM3 , superconducting
qubit sc , and mechanical resonator MR , respectively. Our results reveal more
detailed properties of thermodynamic irreversibility that are stronger than
the conventional SLT, for given a restricted class of irreversible processes.
Along with other studies addressing squeezing effects in quantum
thermodynamics, we hope that our analyses help to unveil the role of squeezing
effects in quantum thermodynamics devices.

In summary, we studied the influence of squeezed characteristics of bath on
the IEP of open quantum systems. The results show that the equilibrium rates
of IEP and its bounds become faster, and the values of IEP are reduced through
harvesting the benefits of squeezing effects in the case of both dissipation
model and dephasing model. In the dissipation model, the summarized overall
trends of geometrical bounds (LB and UB) found in reference fun1  also hold
for squeezed thermal bath, while only the UB becomes tighter in the long-time
limit, and the LB exhibits subtle tightness in the early stage of evolution,
compared with the thermal bath. Unlike the dissipation model, both UB and LB
increase monotonically as the squeezing character of the bath grows, and only
the LB becomes tighter (compared with the thermal bath) during the whole
dynamic process in the dephasing model.

Moreover, the concrete amount of thermodynamic bounds greatly depends on the
explicit form of the system-bath coupling, whereas the trends of them is
independent of the details of the models. Also, the above-summarized trends
for the bounds are independent of system size and hold for systems having more
degrees of freedom. Our results do not contradict the SLT, which is modified
by the inclusion of squeezing as an available resource in the bath. It is
worth noting that a general evolution and the associated geometrical bounds of
irreversibility of two-level system in thermal bath were theoretically
analyzed and experimentally demonstrated in fun1 . Here, we further highlight
the role of adjustable parameters in bath, such as temperture and squeezing
degree, on the reduction of thermodynamic irreversibility. It is expected that
the present work helps in developing a better understanding of the
irreversibility under ambient conditions.

## ACKNOWLEDGMENT

The authors thank Jun-Hong An and Zhang-Qi Yin enlightening discussions. This
work is supported by the Hubei Province Science Fund for Distinguished Young
Scholars under Grant No. 2020CFA078, and by the Special Project for Research
and Development in Key areas of Guangdong Province under Grant No.
2020B0303300001, National Natural Science Foundation of China under Grants No.
U21A20434, No. 12074390, No. 11835011, No. 11734018. JBY acknowledges supports
from A*STAR Career Development Award (C210112010), A*STAR (#21709), and
National Research Foundation Singapore via Grant No. NRF2021-QEP2-02-P01.

## References

  * (1) Seifert U 2005 Phys. Rev. L 95 040602
  * (2) Landi G T and Paternostro M 2021 Rev. Mod. Phys. 93 035008
  * (3) Batalhão T B, Souza A M, Sarthour R S, Oliveira I S, Paternostro M, Lutz E and Serra R M 2015 Phys. Rev. L 115 190601
  * (4) Campbell S, Guarnieri G, Paternostro M and Vacchini B 2017 Phys. Rev. A 96 042109
  * (5) Singh V and Müstecaplıoğlu Ö 202 Phys. Rev. E 102 062123
  * (6) Deffner S and Lutz E 2010 Phys. Rev. L 105 170402
  * (7) Deffner S and Lutz E 2011 Phys. Rev. L 107 140404
  * (8) Vu T V and Hasegawa Y 2021 Phys. Rev. L 126 010601
  * (9) Vu T V and Hasegawa Y 2021 Phys. Rev. L 127 190601
  * (10) Shiraishi N and Saito K 2019 Phys. Rev. L 123 110603
  * (11) Mancino L, Cavina V, Pasquale A D, Sbroscia M, Booth R I, Roccia E, Gianani I, Giovannetti V and Barbieri M 2018 Phys. Rev. L 121 160602
  * (12) Latune C L, Sinayskiy I and Petruccione F 2020 Phys. Rev. A  102 042220
  * (13) Ruppeiner G 1995 Rev. Mod. Phys. 67 605
  * (14) Alhambra A M and Woods M P 2017 Phys. Rev. A 96 022118
  * (15) Zhang J W, Rehan K, Li M, Li J C, Chen L, Su S L, Yan L L, Zhou F and Feng M 2020 Phys. Rev. R 2 033082
  * (16) Poyatos J F, Cirac J I and Zoller P 1996 Phys. Rev. L 77 4728
  * (17) Hou Q Z, Yang C J, Chen C Y, An J H, Yang W L and Feng M 2019 Phys. Rev. A 100 032302
  * (18) Myatt C J, King B E, Turchette Q A, Sackett C A, Kielpinski D, Itano W M, Monroe C and Wineland D J 2000 Nature (London) 403 269
  * (19) Murch K W, Vool U, Zhou D, Weber S J, Girvin S M and Siddiqi I 2012 Phys. Rev. L 109 183602
  * (20) Krauter H, Muschik C A, Jensen K, Wasilewski W, Petersen J M, Cirac J I and Polzik E S 2011 Phys. Rev. L 107 080503
  * (21) Shankar S, Hatridge M, Leghtas Z, Sliwa K, Narla A, Vool U, Girvin S, Frunzio L, Mirrahimi M and Devoret M 2013 Nature(London) 504 7480
  * (22) Lin Y, Gaebler J P, Reiter F, Tan T R, Bowler R, S ørensen A S, Leibfried D and Wineland D J 2013 Nature (London) 504 415
  * (23) Ro$\beta$nagel J, Abah O, Schmidt-Kaler F, Singer K and Lutz E 2014 Phys. Rev. L 112 030602
  * (24) Long R and Liu W 2015 Phys. Rev. E 91 062137
  * (25) Manzano G, Galve F, Zambrini R and Parrondo J M R 2016 Phys. Rev. E 93 052120
  * (26) Manzano G 2018 Phys. Rev. E 98 042123
  * (27) Huang X L, Wang T and Yi X X 2012 Phys. Rev. E 86 051105
  * (28) Abah O and Lutz E 2014 Europhys. Lett. 106 20001
  * (29) Elouard C, Bernardes N K, Carvalho A R R, Santos M Fand Aufféves A 2017 New J. Phys. 19 103011
  * (30) Klaers J, Faelt S, Imamoglu A and Togan E 2017 Phys. Rev. X 7 031044
  * (31) de Assis R J, Sales J S, da Cunha J A R and de Almeida N G 2020 Phys. Rev. E 102 052131
  * (32) Wang J, He J and Ma Y 2019 Phys. Rev. E 100 052126
  * (33) Orszag M 2007 Quantum Optics 2nd ed (Springer-Verlag, Berlin)
  * (34) Hernandez M and Orszag M 2018 Phys. Rev. A 78 042114
  * (35) Cengel Y A and Boles M A 2001 Thermodynamics: An Engineering Approach (McGraw-Hill, New York)
  * (36) Drummond P D and Ficek Z (eds.) 2004 Quantum Squeezing (Springer-Verlag, Berlin/Heidelberg)
  * (37) Polzik E S 2008 Nature (London) 453 45 (2008)
  * (38) Srikanth R and Banerjee S 2008 Phys. Rev. A 77 012318
  * (39) Chin A W, Huelga S F and Plenio M B 2012 Phys. Rev. L 109 233601
  * (40) Weiss U 2008 Quantum Dissipative Systems (World Scientific, Singapore)
  * (41) Leggett A J, Chakravarty S, Dorsey A, Fisher M, Garg A and Zwerger W 1987 Rev. Mod. Phys. 59 1
  * (42) Breuer H-P and Petruccione F 2002 The Theory of Open Quantum Systems (Oxford University Press, New York)
  * (43) Agarwal G S 2012 Quantum Optics 1st ed. (Cambridge University Press, Cambridge)
  * (44) Schaller G 2014 Open Quantum Systems Far from Equilibrium 1st ed (Springer, New York)
  * (45) Li S-W, Cai C Y and Sun C P 2015 Ann. Phys. 360 19
  * (46) You Y-N and Li S-W Phys. Rev. A 97 012114
  * (47) Banerjee S and Ghosh R 2007 J. Phys. A: Math. Theor 40 13735
  * (48) Chruściński D and Kossakowski A 2012 J. Phys. B 45 154002
  * (49) Gardiner C and Zoller P 2004 Quantum Noise (Springer New York)
  * (50) Haikka P, McEndoo S, De Chiara G, Palma G M and Maniscalco S 2011 Phys. Rev. A 84 031602
  * (51) Haikka P, Johnson T H and Maniscalco S 2013 Phys. Rev. A 87 010103(R)
  * (52) Lostaglio M, Jennings D and Rudolph T 2015 Nat. Commun. 6 6383 (2015)
  * (53) Miller H J D, Scandi M, Anders J and Perarnau-Llobet M 2019 Phys. Rev. L 123 230603
  * (54) Scandi M, Miller H J D, Anders J and Perarnau-Llobet M 2020 Phys. Rev. R 2 023377
  * (55) Smith A, Lu Y, An S, Zhang X, Zhang J-N, Gong Z, Quan H T, Jarzynski C and Kim K 2018 New J. Phys. 20 013008
  * (56) Hekking F W J and Pekola J P 2013 Phys. Rev. L 111 093602
  * (57) Plastina F, Alecce A, Apollaro T J G, Falcone G, Francica G, Galve F, Gullo N L and Zambrini R 2014 Phys. Rev. L 113 260601
  * (58) Varizi A D, Cipolla M A, Perarnau-Llobet M, Drumond R C and Landi G T 2021 New J. Phys. 23 063027
  * (59) Santos J P, Céleri L C, Landi G T and Paternostro M 2019 npj Quantum Inf. 5 23
  * (60) Francica G, Goold J and Plastina F 2019 Phys. Rev. E 99 042105
  * (61) Camati P A, Santos J F G and Serra R M 2019 Phys. Rev. A 99 062103
  * (62) Santos J P, Landi G T and Paternostro M 2017 Phys. Rev. L 118 220601
  * (63) Brunelli M, Fusco L, Landig R, Wieczorek W, Hoelscher-Obermaier J, Landi G,Semião F L,Ferraro A, Kiesel N, Donner T, De Chiara G and Paternostro M 2018 Phys. Rev. L 121 160604
  * (64) Spohn H 1987 J. Math. Phys. 19 1227
  * (65) Spohn H and Lebowitz J L 1978 Adv. Chem. Phys. 38 109
  * (66) Alicki R 1979 J Phys. A: Math. Gen 12 L103
  * (67) Pietzonka P, Barato A C and Seifert U 2016 Phys. Rev. E 93 052145
  * (68) Pietzonka P and Seifert U 2018 Phys. Rev. L 120 190602
  * (69) Guarnieri G, Landi G T, Clark S R and Goold J 2019 Phys. Rev. R 1 033021
  * (70) Timpanaro A M, Guarnieri G, Goold J and Landi G T 2019 Phys. Rev. L 123 090604
  * (71) Okuyama M and Ohzeki M 2018 Phys. Rev. L 120 070402
  * (72) Shiraishi N, Funo K and Saito K 2018 Phys. Rev. L 121 070601
  * (73) Murch K W, Weber S J, Beck K M, Ginossar E and Siddiqi I 2013 Nature (London) 499 62
  * (74) Eddins A, Kreikebaum J M, Toyli D M, Levenson-Falk E M, Dove A, Livingston W P, Levitan B A, Govia L C G, Clerk A A and Siddiqi I 2019 Phys. Rev. X 9 011004
  * (75) Toyli D M, Eddins A W, Boutin S, Puri S, Hover D, Bolkhovsky V, Oliver W D, Blais A and Siddiqi I 2016 Phys. Rev. X 6 031004
  * (76) Wollman E E, Lei C U, Weinstein A J, Suh J, Kronwald A, Marquardt F, Clerk A A and Schwab K C 2015 Science 349 952
  * (77) Clark J B, Lecocq F, Simmonds R W, Aumentado J and Teufel J D 2017 Nature (London) 541 191
  * (78) Poyatos J F, Cirac J I and Zoller P 1996 Phys. Rev. L 77 4728
  * (79) Myatt C J, King B E, Turchette Q A, Sackett C A, Kielpinski D, Itano W M, Monroe C and Wineland D J 2000 Nature (London) 403 269
  * (80) Lin Y, Gaebler J P, Reiter F, Tan T R, Bowler R, Sørensen A S, Leibfried D and Wineland D J 2013 Nature 504 415
  * (81) Harrington P M, Tan D, Naghiloo M and Murch K W Phys. Rev. L 123 020502
  * (82) Rossi M, Mancino L, Landi G T, Paternostro M, Schliesser A and Belenchia A 2020 Phys. Rev. L 125 080601