# Scalar Leptoquark Contribution to Two Photon Elastic Scattering

Muhamad Omar Nadeem

Abstract

We study the contribution of Scalar Leptoquark (SLQ) loops to the two photon
scattering cross section at high energies. The leading order helicity
amplitudes are discussed at coupling $\mathcal{O}(e^{4})$ for Standard Model
(SM) fermions, weak bosons and SLQs. Helicity amplitudes for the different
types of loops are calculated for photon centre of mass energies being much
higher than weak boson mass scale but still much lower than the predicted
Leptoquark mass scales. Finally, the cross section due to the SLQ loops is
presented in conjunction with the known SM contribution.

##  1 INTRODUCTION

Among the several candidates of new physics interactions, the notion of direct
couplings of quarks and leptons has remained persistent throughout the years
and is gaining renewed traction in recent times. While such phenomena are
forbidden in the known SM, a certain class of new physics models do allow for
such interactions by postulating the existence of bosonic particles called
Leptoquarks which couple directly with quarks and leptons. This can be
implemented either through an entirely different gauge group or an extension
of the particle content of the known SM, or even a combination of the two.
The idea of Leptoquarks has been around since before 1974 when they appeared
in the Pati-Salam Model [1]. In this paper we shall be considering an
extension of the SM to include SLQs, as presented in [2]. It postulates the
existence of four types of massive scalar particles which carry both colour
and electric charges. These Leptoquarks can couple up and down type quarks
both among themselves and with leptons and their corresponding neutrinos. The
bosonic sector of the Lagrangian also receives new interactions as Leptoquarks
can couple to all of the SM bosons. This allows for a myriad of new topologies
which have been used to explain anomalous phenomenon such as Lepton Flavour
Universality violation [3, 4, 5] and the $W$ boson mass anomaly [6].
The purpose of this paper is to explore the SLQ-Photon interactions in this
model [2] by analysing the SLQ loops’ contribution to the cross section of
elastic two-photon scattering at high energies. The
$\gamma\gamma\rightarrow\gamma\gamma$ channel has been touted as useful tool
to reveal new physics [7] due to the relative ease of producing and
controlling photons, as compared to $W$, $Z$ or Higg’s bosons. The
$\gamma\gamma\rightarrow\gamma\gamma$ channel also has a much cleaner
background and fewer new parameters as compared to decay processes of heavier
particles [7]. As such this channel has been applied to the searches for
Unparticles [8] and investigating low scale quantum gravity [9, 10].
Recent studies [11] place the mass of Leptoquarks well into the TeV scale so
to appreciably detect any Leptoquark loop effects in the
$\gamma\gamma\rightarrow\gamma\gamma$ cross section we require higher photon
energies. This is becoming an increasingly attractive possibility due to the
Photon Linear Collider (PLC) [12] and LHC [13] being capable of achieving
photon centre of mass energies well in excess of 200 GeV.
The layout of this paper is as follows: In Section 2 we present the kinematics
of the process and in Section 3 we use it to calculate the leading helicity
amplitudes in both SM and SLQ model. Section 4 presents the scattering cross
section and in Section 5 we compare, via plots, the contributions of the
various helicity amplitudes to the cross section. Finally in Section 6 we give
our conclusion and discuss the possible extensions to this work.
All Feynman diagrams in this paper have been made using the software
JaxoDraw [14].

##  2 Kinematics

The scattering process is considered in the centre of mass frame with the two
initial photons colliding head on and the scattered photons going nearly
collinear to the initial incident direction. Hence the cross section is
analyzed at small scattering angle.
We shall be borrowing from the kinematics of [15] which presents the momenta
and polarization vectors in the photons’ centre of mass frame. The photon
momenta are

$\begin{split}k_{1}^{\mu}&=(E_{\gamma},0,0,E_{\gamma})\;,\\\
k_{2}^{\nu}&=(E_{\gamma},0,0,-E_{\gamma})\;,\\\
k_{3}^{\sigma}&=(E_{\gamma},0,E_{\gamma}\sin{\theta},E_{\gamma}\cos{\theta})\;,\\\
k_{4}^{\lambda}&=(E_{\gamma},0,-E_{\gamma}\sin{\theta},-E_{\gamma}\cos{\theta})\;.\end{split}$
(1)

Here we keep with the convention of associating a particular Lorentz index
with each momentum four vector. The circularly polarized vectors for each
photon are :

$\begin{split}\varepsilon_{\mu}^{1,\pm}&=\frac{1}{\sqrt{2}}(0,1,\pm i,0)\;,\\\
\varepsilon_{\nu}^{2,\pm}&=\frac{1}{\sqrt{2}}(0,-1,\pm i,0)\;,\\\
\varepsilon_{\sigma}^{3,\pm}&=\frac{1}{\sqrt{2}}(0,1,\pm i\cos{\theta},\mp
i\sin{\theta})\;,\\\ \varepsilon_{\lambda}^{4,\pm}&=\frac{1}{\sqrt{2}}(0,1,\mp
i\cos{\theta},\pm i\sin{\theta})\;.\end{split}$ (2)

Considering the energies at hand and the possible scale of Leptoquark masses,
we do the calculations with a scale hierarchy. Leptoquark mass scale is much
larger than photon energies hence the ratio
$k_{i}.k_{j}/\hat{\mathbb{M}_{Q}^{2}}$ is considered a small parameter and its
higher powers will result in the suppression of certain amplitudes and terms.

##  3 Feynman Amplitudes

The general form of all two photon scattering amplitudes is thus

$\mathcal{M}=\varepsilon^{*}_{\lambda}\varepsilon^{*}_{\sigma}\varepsilon_{\nu}\varepsilon_{\mu}\mathcal{M}^{\mu\nu\sigma\lambda}\;,$
(3)

with the Lorentz indices switched around for the different topologies. The
rank four tensor structure $\mathcal{M}^{\mu\nu\sigma\lambda}$ shows up for
all SM and SLQ contributions to the $\gamma\gamma\rightarrow\gamma\gamma$
amplitude. The leading order contribution $(e^{4})$ includes all six flavours
of the SM quarks and the three charged leptons which all form one loop box
diagrams. Further contributions are from the $W^{\pm}$ bosons which give one
loop bubble, triangle and box diagrams due to the $W^{\pm}\gamma$ triplet and
quartic vertices. The four SLQs also give bubble, triangle and box one loop
diagrams, as in Figure 1, with exactly the same coupling order as the SM
particles’ loops. Each amplitude will be written in helicity basis, determined
by the helicity eigen state of each of the four participating photons.

$\mathcal{M}^{r_{1},r_{2},r_{3},r_{4}}=\mathcal{M}^{\mu\nu\sigma\lambda}\varepsilon_{\lambda}^{r_{4}*}\varepsilon_{\sigma}^{r_{3}*}\varepsilon_{\nu}^{r_{2}}\varepsilon_{\mu}^{r_{1}}\;,$
(4)

where $r_{i}=\pm$ denotes the helicity state of each photon. Due to parity
invariance, time reversal symmetry and Bose symmetry , only some of the
sixteen amplitudes are distinct [12, 15]. Those being $\mathcal{M}^{++++}$,
$\mathcal{M}^{++--}$, $\mathcal{M}^{+-+-}$, $\mathcal{M}^{+--+}$ and
$\mathcal{M}^{+++-}$.

Figure 1: An example each of the (a) Box, (b) Triangle and (c) Bubble scalar
Leptoquark loops in two photon scattering. Wiggly lines denote the photon and
straight lines the Scalar Leptoquarks.

###  3.1 Standard Model Amplitudes

Even though our aim here is to investigate the effect of SLQ loops on the
total cross section, SM amplitudes are essential to it because the observable
cross section is proportional to $|\mathcal{M}|^{2}$ where
$\mathcal{M}=\mathcal{M}_{SM}+\mathcal{M}_{SLQ}$.
The thesis [15] gives an excellent treatment of low energy fermionic
contributions to the $\gamma\gamma\rightarrow\gamma\gamma$ cross section but
it is only valid for photon energies much smaller than the light fermionic
mass scales and hence unsuitable for our purposes where we consider energies
much higher than $m_{W}$. As such the scale dependant fine structure is taken
to be $\alpha=1/128$ and $m_{w}=80.22\text{GeV}$ [12]. Above this scale
fermionic contributions are dominated by $W^{\pm}$ loop contributions [12]. We
shall instead use the high energy approximation to both the fermionic and
$W^{\pm}$ bosonic loops as given in [12]. For photon energies much higher than
the weak boson mass scale and at a fixed small scattering angle, taken to be
$\theta=0.1^{\circ}$ the fermionic and bosonic leading helicity amplitudes are

$\mathcal{M}_{W}^{++++}=\mathcal{M}_{W}^{+-+-}\;,$ (5)

$\mathcal{M}_{f}^{++++}=\mathcal{M}_{f}^{+-+-}\;,$ (6)

$\mathcal{M}_{W}^{++++}=16i\alpha^{2}\pi\frac{k_{1}.k_{2}}{k_{1}.k_{3}}\ln\left(\frac{k_{1}.k_{2}}{k_{1}.k_{3}}\right)\;,$
(7)

$\mathcal{M}_{f}^{++++}=-4\alpha^{2}\ln^{2}\left(\frac{k_{1}.k_{2}}{k_{1}.k_{3}}\right)\;.$
(8)

###  3.2 Scalar Leptoquark Amplitudes

Figure 2: Each of these bubble diagrams has an equivalent accompanying
diagram. This is because the loop consists of a positive SLQ and an equally
negatively charged SLQ. In each case the accompanying diagram is due to the
momenta swapped between the postive and negative SLQs.

In keeping with the condition that the energy of the photons is much smaller
than TeV scale ie $E_{\gamma}<<\hat{\mathbb{M}}_{Q}$, SLQ box loop
contributions can be dropped as all the terms in those integrals scale as
$1/\hat{\mathbb{M}}_{Q}^{4}$. SLQ triangle loops have a larger contribution
but still smaller than SLQ bubble loops due to the significantly smaller pre-
factor coefficients of the two triplet vertices [2] in the triangle loop.
Bubble loops are the ones giving the leading order contribution to the
possible new physics in this process. Even so, the leading order term of the
bubble loop is not interesting from a dynamical point of view as it has no
dependance on photon energies and scattering angles. Those terms only show up
in the subleading terms, scaling as $k_{i}.k_{j}/\hat{\mathbb{M}}_{Q}^{2}$.
Using the Feynman rules of [2] one of the six bubble loops (Figure 2) is

$\mathcal{M}_{B1}^{\mu\nu\sigma\lambda}=\sum_{X,Q}\mu^{2\epsilon}\int\frac{d^{D}l}{(2\pi)^{D}}\frac{\text{Tr}[(ie^{2}\frac{X}{9})g^{\mu\nu}(ie^{2}\frac{X}{9})g^{\sigma\lambda}]}{[l^{2}-\hat{\mathbb{M}_{Q}^{2}}][(k_{1}+k_{2}-l)^{2}-\hat{\mathbb{M}_{Q}^{2}}]}\;.$
(9)

The sum is implied over the four different types of SLQs in this model with
$Q$ being their electric charge and $X$ being a number in the vertex factor
for each type of Leptoquark. Here $Q=\pm 1/3,\pm 2/3,\pm 4/3,\pm 5/3$ and
$X=2,8,32,50$ respectively. Expanding the fully integrated result of (9) in
powers of the regulator $\epsilon$, using the package [16], gives

$\begin{split}\mathcal{M}_{B1}^{\mu\nu\sigma\lambda}&=\sum_{X,Q}-4i\alpha^{2}\left(\frac{X}{9}\right)^{2}g^{\mu\nu}g^{\sigma\lambda}\times\Biggl{[}\frac{1}{\epsilon}+1+\ln\left(\frac{\mu^{2}}{2\hat{\mathbb{M}_{Q}^{2}}}\right)\\\
&+\frac{1}{2}\ln\left(1+\frac{k_{1}.k_{2}}{2\hat{\mathbb{M}_{Q}^{2}}}\right)-\sqrt{-\frac{2\hat{\mathbb{M}_{Q}^{2}}}{k_{1}.k_{2}}}\text{ArcTanh}\left(\sqrt{-\frac{k_{1}.k_{2}}{2\hat{\mathbb{M}_{Q}^{2}}}}\right)\Biggr{]}\;.\end{split}$
(10)

This particular package properly works only for Mathematica versions 12.0 and
higher [17]. Another exactly equivalent loop, with the same tensor structure,
is possible by swapping the loop momenta between the positive and negative
SLQs in each loop. Hence

$\mathcal{M}_{B1}^{\mu\nu\sigma\lambda}=\mathcal{M}_{B2}^{\mu\nu\sigma\lambda}\;.$
(11)

Other loops are

$\begin{split}\mathcal{M}_{B3}^{\mu\sigma\nu\lambda}&=\sum_{X,Q}-4i\alpha^{2}\left(\frac{X}{9}\right)^{2}g^{\mu\nu}g^{\sigma\lambda}\times\Biggl{[}\frac{1}{\epsilon}+1+\ln\left(\frac{\mu^{2}}{2\hat{\mathbb{M}_{Q}^{2}}}\right)\\\
&+\frac{1}{2}\ln\left(1-\frac{k_{1}.k_{3}}{2\hat{\mathbb{M}_{Q}^{2}}}\right)-\sqrt{\frac{2\hat{\mathbb{M}_{Q}^{2}}}{k_{1}.k_{3}}}\text{ArcTanh}\left(\sqrt{\frac{k_{1}.k_{3}}{2\hat{\mathbb{M}_{Q}^{2}}}}\right)\Biggr{]}\;,\end{split}$
(12)
$\mathcal{M}_{B3}^{\mu\sigma\nu\lambda}=\mathcal{M}_{B4}^{\mu\sigma\nu\lambda}\;,$
(13)

and

$\mathcal{M}_{B5}^{\mu\lambda\nu\sigma}=\mathcal{M}_{B6}^{\mu\lambda\nu\sigma}\;.$
(14)

where (13) is exactly the same as (14) except for $k_{1}.k_{3}$ replaced with
$k_{1}.k_{4}$. Hence, from the kinematics of the previous section, it is clear
that (9) only has photon energy dependance while (10) and (11) have dependance
on the scattering angle as well.
The SLQ bubbles can be renormalized by invoking a renormalization condition
that SLQ loop amplitudes are negligible at energies much smaller than electron
mass scale and by expanding up to leading power in
$x_{ij}=k_{i}.k_{j}/\hat{2\mathbb{M}_{Q}^{2}}$, the renormalized bubble loops
become

$\mathcal{M}_{B1,R}^{\mu\nu\sigma\lambda}=\sum_{X,Q}-4i\alpha^{2}\left(\frac{X}{9}\right)^{2}g^{\mu\nu}g^{\sigma\lambda}\left[\frac{x_{12}}{6}\right]\;,$
(15)
$\mathcal{M}_{B3,R}^{\mu\sigma\nu\lambda}=\sum_{X,Q}-4i\alpha^{2}\left(\frac{X}{9}\right)^{2}g^{\mu\sigma}g^{\nu\lambda}\left[-\frac{5}{6}x_{13}\right]\;,$
(16)
$\mathcal{M}_{B5,R}^{\mu\lambda\nu\sigma}=\sum_{X,Q}-4i\alpha^{2}\left(\frac{X}{9}\right)^{2}g^{\mu\lambda}g^{\nu\sigma}\left[-\frac{5}{6}x_{14}\right]\;.$
(17)

Adding all the above bubble amplitudes gives

$\begin{split}\mathcal{M}_{B,R}^{r_{1},r_{2},r_{3},r_{4}}&=2(4i\alpha^{2})\left(\sum_{X}\left(\frac{X}{9}\right)^{2}\right)\\\
&\times\Biggl{[}\left(-\frac{x_{12}}{6}\right)[(\varepsilon_{1}^{r_{1}}.\varepsilon_{2}^{r_{2}})(\varepsilon_{3}^{r_{3}*}.\varepsilon_{4}^{r_{4}*})]\\\
&+\left(\frac{5}{6}x_{13}\right)[(\varepsilon_{1}^{r_{1}}.\varepsilon_{3}^{r_{3}*})(\varepsilon_{2}^{r_{2}}.\varepsilon_{4}^{r_{4}*})]\\\
&+\left(\frac{5}{6}x_{14}\right)[(\varepsilon_{1}^{r_{1}}.\varepsilon_{4}^{r_{4}*})(\varepsilon_{2}^{r_{2}}.\varepsilon_{3}^{r_{3}*})]\Bigg{]}\quad.\end{split}$
(18)

##  4 Scattering Cross Section

The unpolarized differential scattering cross section for
$\gamma\gamma\rightarrow\gamma\gamma$ is given by [15]

$\frac{d\sigma}{d\Omega}=\frac{|\mathcal{M}|^{2}}{256\pi^{2}E_{\gamma}^{2}}\;,$
(19)

where sum is implied over the two helicity states of each photon. To get the
most significant term in the cross section, we combine the leading SM and SLQ
contributions in each amplitude.

$\begin{split}\mathcal{M}^{++++}&=\mathcal{M}_{W}^{++++}+\mathcal{M}_{SLQ,B}^{++++}\;,\\\
\mathcal{M}^{+-+-}&=\mathcal{M}_{W}^{+-+-}\;,\\\
\mathcal{M}^{++--}&=\mathcal{M}_{SLQ,B}^{++--}\;.\end{split}$ (20)

Hence

$|\mathcal{M}|^{2}=\frac{1}{4}(2|\mathcal{M}^{++++}|^{2}+2|\mathcal{M}^{++--}|^{2}+2|\mathcal{M}^{+-+-}|^{2})\;.$
(21)

It can also be seen that the leading order new physics effects come from the
cross terms of the total amplitude mod square, in the following form

$|\mathcal{M}|^{2}_{W,SLQ}=\mathcal{M}_{W}^{++++}(\mathcal{M}_{SLQ,B}^{++++})^{*}+(\mathcal{M}_{W}^{++++})^{*}\mathcal{M}_{SLQ,B}^{++++}+...\;.$
(22)

The $W^{\pm}$ loop amplitude acts as the large coefficient for the much
smaller SLQ amplitude.

##  5 Numerical Results

Figures 3 to 5 show the variation of the positive real values of the various
helicity amplitudes contributing to this process. The general trend is that
amplitudes increase upon increasing photon energies but the effect is tempered
in the total differential cross section by the quadratic photon energy term in
the denominator . All the plots in Figures 3 to 6 have been made under the
assumption that all four types of SLQs have equal masses of order 3 TeV. Only
small corrections are noted by the addition of SLQ loops or by varying the
predicted mass of SLQs. Nevertheless, such things do have an effect on the
cross section, which is most evident by analyzing the cross section
contribution from new physics effects as given in Figure 7. They provide an
increase of the order $10^{-2}$ to the SM cross section. Figure 7 also shows
the leading corrections to the cross section for three different cases of SLQ
masses. It shows that higher values for the SLQ masses will result in smaller
contributions to the cross section. By giving each SLQ a different mass, lying
in the 2-5 TeV range, the new plot will lie between the 2 TeV and the 5 TeV
plots.

Figure 3: The dominant contribution to the two photon unpolarized scattering
cross section comes from the $W^{\pm}$ amplitudes. Figure 4: Leading new
physics effects are due to two helicities of SLQ bubbles. Figure 5: Suppressed
helicities of bubbles and fermionic loops. Figure 6: The profile of the SLQ
model improved cross section closely follows the one due to the SM, as the
overall magnitude is set by the $W^{\pm}$ contribution. Figure 7: Cross
section contributions from different values of SLQ mass.

##  6 Conclusion

The contribution of SLQ one loop amplitudes to the two photon scattering cross
section is smaller than the $W^{\pm}$ contribution but still significantly
larger than that of SM fermionic loops. The SLQ improved cross section closely
follows the trend of the SM predictions. At energies in the high GeV scale,
the $W^{\pm}$ loop contributions dominate, followed by
$\mathcal{M}_{SLQ,B}^{++--}$ and $\mathcal{M}_{SLQ,B}^{++++}$. Other SLQ
bubble helicities are comparable to the SLQ triangle helicities and fermionic
loops, being several orders of magnitude smaller than the leading terms.
Higher predicted values of the SLQ mass will result in progressively smaller
cross section contributions.
These results can be significantly improved by QCD corrections to the SLQ
loops. While QCD corrections have been applied to SM fermionic loops [18], at
these energies even the leading order one loop fermionic amplitude is
significantly suppressed. The $W$ boson loops on the other hand, only allow
for the much smaller QED corrections. As such, the QCD improved SLQ loops
could come very close to the QED improved $W$ boson loops.

## ACKNOWLEDGEMENTS

I extend my gratitude to my colleagues I.Hussain, B.Altaf and F.Ahmad for
vital tech support. And to Dr F.M.Bhutta and Dr S.Ishaq for several helpful
discussions and feedback.

## References

  * [1] C. Murgui and M. B. Wise, “Scalar leptoquarks, baryon number violation, and pati-salam symmetry,” Physical Review D, vol. 104, Aug. 2021.
  * [2] A. Crivellin and L. Schnell, “Complete lagrangian and set of feynman rules for scalar leptoquarks,” Computer Physics Communications, vol. 271, p. 108188, Feb. 2022.
  * [3] D. Bečirević, O. Sumensari, and R. Z. Funchal, “Lepton flavor violation in exclusive $b\rightarrow s$ decays,” The European Physical Journal C, vol. 76, Mar. 2016.
  * [4] L. Calibbi, A. Crivellin, and T. Li, “Model of vector leptoquarks in view of the $b$-physics anomalies,” Phys. Rev. D, vol. 98, p. 115002, Dec 2018.
  * [5] D. Becirevic, M. Fedele, I. Nisandzic, and A. Tayduganov, “Lepton flavor universality tests through angular observables of $\overline{B}\to d^{(\ast)}\ell\overline{\nu}$ decay modes,” 2019.
  * [6] S.-P. He, “Leptoquark and vector-like quark extended model for simultaneousexplanation of w boson mass and muon g–2 anomalies*,” Chinese Physics C, vol. 47, p. 043102, Apr. 2023.
  * [7] G. Gounaris, P. Porfyriadis, and F. Renard, “Light by light scattering at high energy: a tool to reveal new particles,” Physics Letters B, vol. 452, p. 76–82, Apr. 1999.
  * [8] O. Çakır and K. Ozansoy, “Unparticle searches through $\gamma\>\gamma$ scattering,” The European Physical Journal C, vol. 56, p. 279–285, July 2008.
  * [9] M. S. Berger and B. Zerbe, “Signals for low scale gravity in the process $\gamma\>\gamma\rightarrow Z\>Z$,” Phys. Rev. D, vol. 72, p. 095007, 2005.
  * [10] H. Davoudiasl, “$\gamma\gamma\rightarrow\gamma\gamma$ as a test of weak scale quantum gravity at the nlc,” Physical Review D, vol. 60, Sept. 1999.
  * [11] N. Desai and A. Sengupta, “Status of leptoquark models after lhc run-2 and discovery prospects at future colliders,” 2023.
  * [12] G. Jikia and A. Tkabladze, “Photon-photon scattering at the photon linear collider,” Physics Letters B, vol. 323, p. 453–458, Mar. 1994.
  * [13] D. d’Enterria and G. G. da Silveira, “Observing light-by-light scattering at the large hadron collider,” Physical Review Letters, vol. 111, Aug. 2013.
  * [14] D. Binosi, J. Collins, C. Kaufhold, and L. Theussl, “Jaxodraw: A graphical user interface for drawing feynman diagrams. version 2.0 release notes,” Computer Physics Communications, vol. 180, p. 1709–1715, Sept. 2009.
  * [15] Hussain, I., “Light by Light Scattering in the Low Energy Regime,” 2016.
  * [16] T. Huber and D. Maître, “Hypexp 2, expanding hypergeometric functions about half-integer parameters,” Computer Physics Communications, vol. 178, p. 755–776, May 2008.
  * [17] W. R. Inc., “Mathematica, Version 12.0.”  Champaign, IL, 2019.
  * [18] Z. Bern, A. D. Freitas, A. Ghinculov, H. Wong, and L. Dixon, “Qcd and qed corrections to light-by-light scattering,” Journal of High Energy Physics, vol. 2001, p. 031–031, Nov. 2001.


# Work statistics in slow thermodynamic processes

Jie Gu<EMAIL_ADDRESS>Chengdu Academy of Education Sciences, Chengdu 610036,
China

###### Abstract

We apply the adiabatic approximation to slow but finite-time thermodynamic
processes and obtain the full counting statistics of work. The average work
consists of change in free energy and the dissipated work, and we identify
each term as a dynamical- and geometric-phase-like quantity. An expression for
the friction tensor, the key quantity in thermodynamic geometry, is explicitly
given. The dynamical and geometric phases are proved to be related to each
other via the fluctuation-dissipation relation.

††preprint: AIP/123-QED

##  I Introduction

Reversibility is valid only for infinitely long quasi-static transformations,
in which both dissipation and fluctuations are absent. Despite perfect
efficiency and constancy, these processes are not realistic. Slow but finite-
time processes add non-trivial ingredients in quasi-static processes, while
remaining analytically tractable. In slow thermodynamic processes, the driving
forces are varying slowly compared with relaxation time scales. It implies
that after an exponentially short transient depending on the initial states,
the state is all the time near the instantaneous steady state [1].

There are at least five frameworks for investigating slow thermodynamic
processes. The first framework is a perturbation theory developed in Refs. [2,
1]. The slowness of driving allows a perturbative expansion of the probability
distribution or density matrix in powers of $\tau^{-1}$, where $\tau$ is the
time duration of the process. The second framework is the adiabatic response
theory [3, 4], which is an extension of the well-established dynamic linear
response theory [5]. It exploits the property that states are close to
equilibrium. Both frameworks reveal a geometric structure of thermodynamic
processes in small systems [6, 4, 7, 8, 9, 10, 11, 12, 13, 14], extending the
concept of thermodynamic length for macroscopic endoreversible thermodynamics
[15, 16, 17, 18, 19, 20, 21]. The key notion is the so-called _friction
tensor_ [15, 3, 4, 2] (also called _thermodynamic metric tensor_ [11, 12,
14]). The dissipated work during $[0,\tau]$ is expressed in terms of the
friction tensor $\bm{g}$ as

$W_{\text{diss}}=\int_{0}^{\tau}\dot{\bm{\Lambda}}\cdot\bm{g}\cdot\dot{\bm{\Lambda}}^{\text{T}}\text{d}t,$
(1)

where $W_{\text{diss}}=\langle W\rangle-\Delta F$, i.e., the dissipated work
$W_{\text{diss}}$ is the difference between average work $\langle W\rangle$
done on the system and the change in equilibrium free energy $\Delta F$, and
$\bm{\Lambda}$ denotes time-dependent control parameters. By the Cauchy-
Schwarz inequality, the dissipated work satisfies

$\tau W_{\text{diss}}\geq\mathcal{L}^{2},$ (2)

where
$\mathcal{L}=\int_{0}^{\tau}\sqrt{\dot{\bm{\Lambda}}\cdot\bm{g}\cdot\dot{\bm{\Lambda}}^{\text{T}}}\text{d}t$
is the thermodynamic length as defined by the metric.

The third framework treats a slow process as a sequence of quench-and-
relaxation processes where the number of steps is large but finite. Its
equivalence with the previous two frameworks has been proved in Ref. [10]. The
fourth is stochastic thermodynamics [22, 23]. As an nonperturbative formalism,
its applicability is not limited to slow processes. For a system in contact
with a heat bath with inverse temperature $\beta$ and undergoing slow
processes, the work distribution is Gaussian, and the Jarzynski equality [24,
25], one of the most notable results in this framework, gives the fluctuation-
dissipation relation (FDR)

$W_{\text{diss}}=\frac{1}{2}\beta\sigma_{W}^{2},$ (3)

where the dissipated work is related to the variance of work $\sigma_{W}^{2}$.
The FDR in this form has also been derived by the third framework [10] and a
variant of the perturbation theory (the first framework) [26].

The fifth framework applies the adiabatic approximation in quantum mechanics
to classical stochastic processes for the full counting statistics of
thermodynamic quantities [27, 28]. It was introduced for classical chemical
kinetics [27], but has been soon extended to a geometrical expression of
excess entropy production [29], particle statistics in quantum transport
problems[30, 31, 32, 33] and quantum heat engines [34], in which Berry-phase-
induced effects are the focus [35]. Nonetheless, this simple but powerful
framework has not been applied to work statistics to the best of our
knowledge. In addition, the connections between the last framework and the
other four have not been explored. Both the Berry’s geometric phase and the
thermodynamic metric tensor are geometric quantities. Is there any connection
between them? Since the fifth framework enables one to calculate the full
counting statistics of work, the variance of work can in principle be
obtained, so what role does the fluctuation-dissipation relation play in this
framework? In this study, we will answer these two questions and unveil the
connections among different frameworks.

##  II Full Counting Statistics of Work

Consider a system with $N$ energy levels $\\{\epsilon_{m}(t)\\}$ interacting
with a heat bath with fixed inverse temperature $\beta$. The energy levels
$\\{\epsilon_{m}(t)\\}$ are varying slowly by external manipulations, serving
as control parameters of the protocol. Assume the system dynamics follows a
stochastic process modeled by a continuous-time Markov chain. The time
evolution of the probability distribution
${{|p(t)\rangle}}=[p_{1},p_{2},\ldots,p_{N}]^{\text{T}}$ is described by a
Pauli master equation [36, 37]

$\frac{\text{d}|{p}(t)\rangle}{\text{d}t}=\bm{R}(t)|p(t)\rangle,$ (4)

where $p_{m}$, the $m$-th entry of $|p\rangle$, is the probability of state
$m$; $R_{mn}\ (n\neq m)$, the $mn$-th entry of the matrix $\bm{R}$, is the
time-dependent transition rate from state $n$ to $m$, and
$R_{mm}=-\sum_{n,n\neq m}R_{nm}$. We assume that the Markov chain is ergodic
and the detailed balance condition

$\frac{R_{mn}}{R_{nm}}=e^{\beta(\epsilon_{n}-\epsilon_{m})}$ (5)

is satisfied, so at any instant there exists an instantaneous Gibbs state
$|\lambda(t)\rangle$ satisfying $\bm{R}(t)|\lambda(t)\rangle=\bm{0}$.

The full counting statistics of work can be obtained as follows. Let
$P(W,\tau)$ be the probability of having done work $W$ on the system within
time $[0,\tau]$. The full counting statistics of work is characterized by the
cumulant generating function $\mathcal{K}_{\chi}(\tau)$ and moment generating
function $\mathcal{M}_{\chi}(\tau)$ [36],

$\mathcal{K}_{\chi}(\tau)=\ln\mathcal{M}_{\chi}(\tau)\equiv\ln\int_{-\infty}^{\infty}e^{i\chi
W}P(W,\tau)\text{d}W,$ (6)

where $\chi$ is the counting field. Then the average and variance of work are
given by, respectively,

$\langle
W\rangle=\mathcal{K}^{\prime}\equiv\left.\frac{\partial\mathcal{K}_{\chi}(\tau)}{\partial(i\chi)}\right|_{\chi=0}$
(7)

and

$\sigma^{2}_{W}=\mathcal{K}^{\prime\prime}\equiv\left.\frac{\partial^{2}\mathcal{K}_{\chi}(\tau)}{\partial(i\chi)^{2}}\right|_{\chi=0}.$
(8)

In order to find the cumulant generating function, let $P_{m}\left(w,t\right)$
denote the probability that within time $t$ the work done to system is $w$ and
the system is dwelling on the $m$-th state. Then $P_{m}\left(w+dw,t+dt\right)$
consists of two contributions given by

$\displaystyle P_{m}\left(w+\text{d}w,t+\text{d}t\right)$ (9)
$\displaystyle=P_{m}\left(w,t\right)\left(1+R_{mm}\text{d}t\right)+\sum_{n(\neq
m)}P_{n}\left(w+\text{d}w,t\right)R_{mn}\text{d}t.$

The first term represents that the system dwell on the $m$-th state during
$[t,t+\text{d}t)$ , and the energy level $m$-th state is manipulated so the
work $\text{d}w=\text{d}\epsilon_{m}$ is performed on the system on top of
$w$. The second term represents jumps from other states, during which no work
is done. Up to the first order of $\text{d}t$, the equation becomes

$\frac{\partial P_{m}}{\partial w}\dot{\epsilon}_{m}+\frac{\partial
P_{m}}{\partial t}=P_{m}R_{mm}+\sum_{n(\neq m)}P_{n}R_{mn}$ (10)

By multiplying $e^{i\chi w}$ on both sides and integrating with respect to
$w$, the time evolution of the state-resolved moment generating function is
given by

$\frac{\partial\mathcal{M}_{\chi,m}(t)}{\partial
t}=\sum_{n}{R}_{\chi,mn}(t)\mathcal{M}_{\chi,n}(t)$ (11)

where

${R}_{\chi,mn}(t)=R_{mn}(t)+i\chi\dot{\epsilon}_{m}(t)\delta_{mn},$ (12)

or in matrix form,

$\bm{R}_{\chi}(t)=\bm{R}(t)+i\chi\dot{\bm{H}}(t),$ (13)

where the Hamiltonian $\bm{H}(t)=\text{diag}\\{\epsilon_{m}(t)\\}.$ Similar
results for a two-state system has been obtained in Ref. [38].

##  III Main results

For later use, we define the eigenvalue of $\bm{R}_{\chi}(t)$ with the largest
real part as $\lambda_{\chi}(t)$, with the corresponding left and right
eigenvectors being $\langle\lambda_{\chi}(t)|$ and $|\lambda_{\chi}(t)\rangle$
by a slight abuse of notation. The left and right eigenvectors associated with
same eigenvalue are not Hermitian-conjugate of each other due to the non-
Hermiticity of $\bm{R}_{\chi}$. Evaluating at $\chi=0$ gives
$\lambda_{\chi=0}=0$,
$\langle\lambda|\equiv\langle\lambda_{\chi=0}(t)|=[1,1,\ldots,1]$ and
$|\lambda(t)\rangle\equiv|\lambda_{\chi=0}(t)\rangle=|p^{\text{G}}(t)\rangle$,
the instantaneous Gibbs state, i.e.,

$|\lambda(t)\rangle=\frac{1}{Z}\big{[}e^{-\beta\epsilon_{1}(t)},e^{-\beta\epsilon_{2}(t)},\ldots,e^{-\beta\epsilon_{N}(t)}\big{]}^{\text{T}},$
(14)

where $Z=\text{Tr}[\exp(-\beta\bm{H})]$ is the partition function.

The cumulant generating function is then

$\mathcal{K}_{\chi}(\tau)=\ln\sum_{m}\mathcal{M}_{\chi,m}(\tau)=\ln\langle\lambda|\mathcal{M}_{\chi}(\tau)\rangle,$
(15)

where we have defined
$|\mathcal{M}_{\chi}(\tau)\rangle=[\mathcal{M}_{\chi,1},\mathcal{M}_{\chi,2},\ldots,\mathcal{M}_{\chi,N}]^{\text{T}}$.

Given the differential equations Eq. (11), the cumulant generating function
can be formally solved by
$\mathcal{K}_{\chi}(\tau)=\ln\langle\lambda|\mathbb{T}\big{[}\exp({\int_{0}^{\tau}\bm{R}_{\chi}(t)dt}\big{]}|p(0)\rangle$,
where $\mathbb{T}$ is the time-ordering operator, and $|p(0)\rangle$ is the
initial state. We consider slow processes and assume that the initial state is
the Gibbs state $|p(0)\rangle=|\lambda(0)\rangle$. Given the similarity
between Eq. (11) and the Schrödinger equation, by invoking the adiabatic
approximation analogous to that in quantum mechanics and following Refs. [27,
28, 30, 29, 33], we find

$\mathcal{M}_{\chi}(\tau)=e^{\int_{0}^{\tau}\lambda_{\chi}(t)\text{d}t}e^{-\int_{C}\left\langle\lambda_{\chi}|\dot{\lambda}_{\chi}\right\rangle\text{d}t}\left\langle\lambda_{\chi}(0)|\lambda(0)\right\rangle\left\langle\lambda|\lambda_{\chi}(\tau)\right\rangle$
(16)

or equivalently

$\displaystyle\mathcal{K}_{\chi}=\mathcal{K}_{\text{dyn}}+\mathcal{K}_{\text{geo}}+\big{[}\ln\left\langle\lambda_{\chi}(0)|\lambda(0)\right\rangle+\ln\left\langle\lambda|\lambda_{\chi}(\tau)\right\rangle\big{]},$
(17a)
$\displaystyle\mathcal{K}_{\text{dyn}}=\int_{0}^{\tau}\lambda_{\chi}(t)\text{d}t,$
(17b)
$\displaystyle\mathcal{K}_{\text{geo}}=-\int_{0}^{\tau}\langle\lambda_{\chi}|\dot{\lambda}_{\chi}\rangle\text{d}t.$
(17c)

We have omitted the argument $\tau$ for convenience. It is clear that
$\mathcal{K}_{\text{dyn}}$ and $\mathcal{K}_{\text{geo}}$ are analogous to the
dynamical and Berry’s geometric phase in quantum mechanics, respectively [35].
The dynamical-phase part corresponds to a fictitious process during which each
state at an arbitrary instant is the instantaneous Gibbs state, and the
geometric-phase part corresponds to effects on top of the former due to the
slow but finite-time driving. As will be seen shortly,
$\mathcal{K}_{\text{geo}}$ here is not a truly geometric quantity, but we will
still keep the subscript “geo” due to its apparent resemblance to the
geometric phase in quantum mechanics. The boundary terms (the two terms in the
square bracket) are reminiscent of the noncyclic geometric phase [39].

Now we proceed to find the average and variance of work by taking the first
and second derivatives of $\mathcal{K}_{\chi}$. Differentiating both sides of
$\langle\lambda_{\chi}|\bm{R}_{\chi}=\lambda_{\chi}\langle\lambda_{\chi}|$
gives

$\langle\lambda_{\chi}^{\prime}|\bm{R}_{\chi}+\langle\lambda_{\chi}|\bm{R}_{\chi}^{\prime}=\lambda_{\chi}^{\prime}\langle\lambda_{\chi}|+\lambda_{\chi}\langle\lambda_{\chi}^{\prime}|,$
(18)

where we have defined $A_{\chi}^{\prime}\equiv\partial
A_{\chi}/\partial(i\chi)$ and $A^{\prime}\equiv\partial
A_{\chi}/\partial(i\chi)|_{\chi=0}$.

By evaluating Eq. (18) at $\chi=0$ and multiplying both sides by
$|\lambda\rangle$ we obtain the derivative of the eigenvalue,

$\lambda^{\prime}=\langle\lambda|\dot{\bm{H}}|\lambda\rangle=\text{Tr}(\dot{\bm{H}}e^{-\beta\bm{H}})/Z,$
(19)

This result obviously coincides with the time-independent first-order
perturbation theory.

By evaluating Eq. (18) at $\chi=0$ and multiplying by the Drazin inverse
$\bm{R}^{+}$ we obtain

$\langle\lambda^{\prime}|=-\langle\lambda|\dot{\bm{H}}\bm{R}^{+},$ (20)

where we have used the property of Drazin inverse that
$\langle\lambda|\bm{R}^{+}=\bm{0}$ [2, 40], the definition of $\bm{R}_{\chi}$
and the fact that $\lambda_{\chi=0}=0$. See Appendix A for definition and
properties of Drazin inverse.

The derivative of the right eigenvector can be obtained similarly,

$|\lambda^{\prime}\rangle=-\bm{R}^{+}\dot{\bm{H}}|\lambda\rangle.$ (21)

According to Eq. (19), the derivative of the dynamical phase is, up to the
first order,

$\mathcal{K}^{\prime}_{\text{dyn}}=\int_{0}^{\tau}\lambda^{\prime}dt=-\beta^{-1}\ln\frac{Z(\tau)}{Z(0)}.$
(22)

Therefore, we obtain the identity

$\Delta F=\mathcal{K}^{\prime}_{\text{dyn}},$ (23)

where $F(t)=-\beta^{-1}\ln Z(t)$ is the equilibrium free energy. This makes
sense because as mentioned above, the dynamical-phase part corresponds to a
fictitious process during which the system is in equilibrium at any instant,
and for this kind of process, the average work is equal to the change in
equilibrium free energy.

By using the properties of Drazin inverse again, we prove in Appendix A that
the boundary terms do not contribute to $\mathcal{K}^{\prime}$, so $\langle
W\rangle=\mathcal{K}^{\prime}=\mathcal{K}^{\prime}_{\text{dyn}}+\mathcal{K}^{\prime}_{\text{geo}}$.
Since $\langle W\rangle=\Delta F+W_{\text{{diss}}}$, Eq. (23) implies that the
dissipated work $W_{\text{diss}}$ can be identified as the derivative of the
geometric phase,

$W_{\text{diss}}=\mathcal{K}^{\prime}_{\text{geo}}.$ (24)

Plugging Eqs. (20) and (21) into Eq. (17c) gives the dissipated work:

$\displaystyle W_{\text{diss}}=\mathcal{K}^{\prime}_{\text{geo}}$
$\displaystyle=\int_{0}^{\tau}\langle\lambda|\dot{\bm{H}}\bm{R}^{+}|\dot{\lambda}\rangle\text{d}t$
(25)
$\displaystyle=-\beta\int_{0}^{\tau}\langle\lambda|\dot{\bm{H}}\bm{R}^{+}\dot{\bm{H}}|\lambda\rangle\text{d}t,$

where in order to obtain the second line we have used Eq. (14) and the
property of Drazin inverse that $\bm{R}^{+}|\lambda\rangle=\mathbf{0}$. One
can see the symmetric structure in Eq. (25). The counterpart of Eq. (25) in
open quantum systems has been derived in Ref. [9] by the perturbation theory
[1], and reduces to our result for classical states, i.e., when the initial
state is a Gibbs state, and $[H(t),H(t^{\prime})]=0$ such that there is no
drive-induced coherence throughout.

As derived in Appendix A, the dissipated work can be explicitly written as

$W_{\text{diss}}=\mathcal{K}_{\text{geo}}^{\prime}=\int_{0}^{\tau}\sum_{m,n}\dot{\epsilon}_{m}g_{mn}\dot{\epsilon}_{n}\text{d}t,$
(26)

where the friction tensor (thermodynamic metric tensor) is

$g_{mn}=-\frac{\beta\big{[}R_{mn}^{+}e^{-\beta\epsilon_{n}}+R_{nm}^{+}e^{-\beta\epsilon_{m}}\big{]}}{2Z}.$
(27)

Here, $R_{mn}^{+}$ is the $mn$-th entry of $\bm{R}^{+}$. This tensor is
symmetric and positive semi-definite. The explicit expression of the friction
tensor for slow classical thermodynamic processes has not been reported to our
knowledge. From Eqs. (26) and (27), it can be seen that the dissipated power
$\dot{W}_{\text{diss}}$ and entropy production rate in slow processes is time-
local: they only depends on the instantaneous values and derivatives of the
control parameters [3].

Now let us prove the fluctuation-dissipation theorem for this setup.
Differentiating both sides of
$\langle\lambda_{\chi}|\bm{R}_{\chi}=\lambda_{\chi}\langle\lambda_{\chi}|$
twice gives

$\langle\lambda_{\chi}^{\prime\prime}|\bm{R}_{\chi}+2\langle\lambda_{\chi}^{\prime}|\bm{R}_{\chi}^{\prime}=\lambda_{\chi}^{\prime\prime}\langle\lambda_{\chi}|+2\lambda_{\chi}^{\prime}\langle\lambda_{\chi}^{\prime}|+\lambda_{\chi}\langle\lambda_{\chi}^{\prime\prime}|.$
(28)

Evaluating at $\chi=0$ and multiplying by $|\lambda\rangle$ lead to

$\lambda^{\prime\prime}=-2\langle\lambda|\dot{\bm{H}}\bm{R}^{+}\dot{\bm{H}}|\lambda\rangle.$
(29)

In Appendix A we show that, up to the first order of $\tau^{-1}$, the second
derivatives of terms other than the dynamical phase do not contribute. Thus,

$\sigma_{W}^{2}=\mathcal{K}^{\prime\prime}=\mathcal{K}_{\text{dyn}}^{\prime\prime}=\int_{0}^{\tau}\lambda^{\prime\prime}\text{d}t=-2\int_{0}^{\tau}\langle\lambda|\dot{\bm{H}}\bm{R}^{+}\dot{\bm{H}}|\lambda\rangle\text{d}t.$
(30)

Comparing this equation with Eq. (25), we immediately obtain the fluctuation-
dissipation relation in the form of Eq. (3) [26, 8, 9, 10]. Note that the
geometric phase does not vanish in general and still contributes, but it is a
higher order term of $\tau^{-1}$ compared with the dynamical phase.

Next we discuss the condition under which the dissipated work and also the
fluctuation of work vanish. Obviously the dissipated work is vanishingly small
for a quasi-static (infinitely slow) process. Other than the previous trivial
case, it is not hard to see that there is no fluctuation of work if all energy
levels are shifted globally; the populations are unchanged, so there is no
entropy production associated with this process. Mathematically, this
condition implies that the rates of change in energy levels are identical,
i.e., $\dot{\bm{H}}\propto\bm{I}$, where $\bm{I}$ is the identity operator.
Plugging it into Eq. (25), we indeed have $\sigma_{W}^{2}\propto
W_{\text{diss}}=0$, due to the property that
$\langle\lambda|\bm{R}^{+}=\mathbf{0}$ [2, 40]. To prove the necessity of the
condition, we exploit the strong result by Shiraishi _et al._ [41],

$\Sigma\geq\frac{L(|p(0)\rangle,|p(\tau)\rangle)^{2}}{2A\tau},$ (31)

where $\Sigma$ is the entropy production, $\tau$ the time duration of the
process, $L(|p(0)\rangle,|p(\tau)\rangle)$ the trace distance between
$|p(0)\rangle$ and $|p(\tau)\rangle$, and $A$ the dynamical activity. This
speed limit implies that, in order for zero entropy production during any
finite-time processes, the probability distribution does not change with time.
Combining with the detailed balance condition, one can see that the entropy
production vanishes only if the energy level spacing between an arbitrary pair
of energy levels remains constant. Furthermore, it remains an open question
whether vanishing entropy is a sufficient or necessary condition (or neither)
for the fluctuation relation for currents to hold [42], an interesting
connection we will not pursue any further here.

On the other hand, one might recall that, in quantum mechanics, if there is
only one control parameter in the Hamiltonian that is changing with time
periodically, i.e., $\bm{H}(\Lambda(t))=\bm{H}(\Lambda(t+\mathcal{T}))$ where
$\mathcal{T}$ is the period, the geometric phase vanishes for this case [43].
It is also well known that no net charges are transported in a single-
parameter adiabatic pump [44, 30, 31, 32, 33, 34]; during the second half
period, any charge that has flowed during the first will return. Therefore,
one might expect $W_{\text{diss}}=\mathcal{K}_{\text{geo}}^{\prime}=0$ if
there is only one control parameter. Nonetheless, similar reasoning does not
apply to $W_{\text{diss}}$ because the dissipated work and entropy production
are monotonically increasing with time; the integrand in Eq. (26) is not
guaranteed to vanish and in general positive even if there is only one control
parameter (see the next section for a concrete example), so the dissipated
work during the second half period obviously does not cancel the first half.
In other words, $W_{\text{diss}}$ is generally positive even for single-
parameter periodic systems, contradicting the expectation above. This
contradiction can be resolved by the observation that
$\mathcal{K}_{\text{geo}}$ is not truly geometric even if Eq. (17c) shows a
structure similar to the geometric phase. A geometric quantity should be
independent of how the path in the parameter space is traversed, especially of
the changing rate of control parameters (of course as long as the adiabatic
approximation holds). However, $\langle\lambda_{\chi}|$ and
$|\lambda_{\chi}\rangle$ depend on not only the control parameters
$\bm{\Lambda}(t)$, but also the time derivative ${\dot{\bm{\Lambda}}}(t)$ as
can be seen from Eq. (13). Therefore, $\mathcal{K}_{\text{geo}}$ and
$W_{\text{diss}}$ are not geometric quantities. By contrast, the geometric
phase in quantum mechanics [35], geometric-phase-induced transport [44, 30,
31, 32, 33, 34] and the thermodynamic length defined below Eq. (2) are
genuinely geometric.

##  IV Example

As a simple illustration, let us consider a two-state system. The transition
rate matrix is given by

$\bm{R}=\begin{bmatrix}-k_{1}&k_{2}\\\ k_{1}&-k_{2}\end{bmatrix}.$ (32)

The eigensystem associated to the nonzero eigenvalue is

$\displaystyle\lambda_{1}=-\left(k_{1}+k_{2}\right),$ (33)
$\displaystyle|\lambda_{1}\rangle=\begin{bmatrix}1,&-1\end{bmatrix}^{\text{T}},$
$\displaystyle\langle\lambda_{1}|=\frac{1}{k_{1}+k_{2}}\begin{bmatrix}k_{1},&-k_{2}\end{bmatrix},$

so the Drazin inverse is simply given by

$\bm{R}^{+}=\frac{1}{(k_{1}+k_{2})^{2}}\bm{R}.$ (34)

By using the first-order perturbation theory, we find, up to the first order
of $\chi$, the eigensystem of $\bm{R}_{\chi}$ associated to $\lambda_{\chi}$
(the eigenvalue with the largest real part) is:

$\displaystyle\lambda_{\chi}$
$\displaystyle=\frac{i\chi}{2}\left[(\dot{\epsilon}_{1}+\dot{\epsilon}_{2})+\frac{(k_{1}-k_{2})(\dot{\epsilon}_{2}-\dot{\epsilon}_{1})}{k_{1}+k_{2}}\right],$
(35a) $\displaystyle|\lambda_{\chi}\rangle$
$\displaystyle=\frac{1}{k_{1}+k_{2}}\begin{bmatrix}k_{2},&k_{1}\end{bmatrix}^{\text{T}}+\frac{i\chi
k_{1}k_{2}(\dot{\epsilon}_{1}\dot{\epsilon}_{2})}{(k_{1}+k_{2})^{3}}\begin{bmatrix}1,&-1\end{bmatrix}^{\text{T}},$
(35b) $\displaystyle\langle\lambda_{\chi}|$
$\displaystyle=\begin{bmatrix}1,&1\end{bmatrix}+\frac{i\chi(\dot{\epsilon}_{1}-\dot{\epsilon}_{2})}{(k_{1}+k_{2})^{2}}\begin{bmatrix}k_{1},&-k_{2}\end{bmatrix}.$
(35c)

Exploiting the detailed balance condition Eq. (5), we have

$\lambda^{\prime}=\left.\frac{\partial\lambda_{\chi}}{\partial(i\chi)}\right|_{\chi=0}=\frac{e^{-\beta\epsilon_{1}}\epsilon_{1}+e^{-\beta\epsilon_{2}}\epsilon_{2}}{Z}$
(36)

The R.H.S. coincides with $\dot{F}$, the time derivative of equilibrium free
energy, so we have verified Eq. (23), i.e., the dynamical-phase part of
$\mathcal{K}^{\prime}$ is equal to the change in equilibrium free energy.

Substituting Eqs. (35b) and (35c) into Eq. (17c) and using the detailed
balance condition, one can obtain the dissipated power

$\dot{W}_{\text{diss}}=\frac{\beta(\dot{\epsilon}_{1}-\dot{\epsilon}_{2})^{2}k_{1}k_{2}}{(k_{1}+k_{2})^{3}}.$
(37)

It is straightforward to check that plugging the Drazin inverse into Eq. (25)
results in exactly the same result. It verifies the discussion in the previous
section that the dissipated power vanishes if and only if
$\dot{\epsilon}_{1}=\dot{\epsilon}_{2}$. Accordingly, the friction tensor is

$\bm{g}=\frac{\beta k_{1}k_{2}}{(k_{1}+k_{2})^{3}}\begin{bmatrix}1&-1\\\
-1&1\end{bmatrix}.$ (38)

This tensor is positive semi-definite as expected.

Even more concretely, we consider a special two-state system—a model for the
bit erasure [45], whose Hamiltonian is given by

$\bm{H}=\frac{\Lambda(t)}{2}\sigma_{z}.$ (39)

During the population reduction step, the energy level splitting $\Lambda(t)$
is increased from zero to an energy much greater than $\beta^{-1}$. The
transition rates are related to $\Lambda(t)$ as

$k_{1}=\gamma n(\Lambda),\quad k_{2}=\gamma[n(\Lambda)+1],$ (40)

where $\gamma$ is the coupling between the system and bath, and
$n(\Lambda)=(e^{\beta\Lambda}-1)^{{-1}}$ is the Bose-Einstein distribution
function. Plugging the transition rates into Eq. (37) leads to

$\dot{W}_{\mathrm{diss}}=\frac{\beta\dot{\Lambda}^{2}\left(1-e^{-\beta\Lambda}\right)e^{-\beta\Lambda}}{\gamma\left(1+e^{-\beta\Lambda}\right)^{3}}.$
(41)

Even though there is only one control parameter, $W_{\text{diss}}$ remains
positive. We note that exactly the same expression has been derived recently
by the perturbation theory [46].

##  V Conclusion

The adiabatic approximation is a powerful tool not only in quantum mechanics,
but also in slow stochastic processes. With this tool, the following
connections have been unveiled. From Eqs. (23) and (24), one can see that
$\langle W\rangle$ consists of $\Delta F$ and $W_{\text{diss}}$, with each
contribution identified as a dynamical- and geometric-phase-like quantity,
respectively. Furthermore, in Eqs. (26) and (27), we calculated the friction
tensor $\bm{g}$ from the geometric phase, which is consistent with the results
obtained by the first two frameworks described in Sect. I. The fluctuation-
dissipation relation was derived in an elementary and transparent manner,
complementary to the third and fourth frameworks. These results can be
summarized in a relation between the dynamical phase and geometric phase,

$\mathcal{K}^{\prime}_{\text{geo}}=\frac{\beta}{2}\mathcal{K}^{\prime\prime}_{\text{dyn}}.$
(42)

Whether similar relations hold for higher cumulants and/or for other
thermodynamic quantities is an interesting open question. Results in slow
quantum processes can be obtained by applying this framework to quantum master
equations [31, 32, 34], and can be compared with existent results such as Ref.
[8, 9, 10, 11, 12].

## Data Availability Statement

The data that supports the findings of this study are available within the
article.

## Appendix A

The Drazin inverse of a transition rate matrix $\bm{R}$ can be found as
follows [2, 40]. Since $\bm{R}$ is assumed to be diagonalizable, it can be
expanded as

$\bm{R}=\sum_{\xi_{\alpha}\neq
0}\xi_{\alpha}|\xi_{\alpha}\rangle\langle\xi_{\alpha}|,$ (43)

where $|\xi_{\alpha}\rangle$ and $\langle\xi_{\alpha}|$ are, respectively, the
left and right eigenvector of $\bm{R}$ associated with the eigenvalue
$\xi_{\alpha}$. Together with $|\lambda\rangle$ and $\langle\lambda|$ that are
associated with the zero eigenvalue, they satisfy [47]:

  1. 1.

$\langle\lambda|\xi_{\alpha}\rangle=0$,
$\langle\xi_{\alpha}|\lambda\rangle=0$,
$\langle\xi_{\alpha}|\xi_{\beta}\rangle=0$.

  2. 2.

$\langle\xi_{\alpha}|\xi_{\alpha}\rangle=1$,
$\langle\lambda|\lambda\rangle=1$.

  3. 3.

$\sum_{\xi_{\alpha}\neq
0}|\xi_{\alpha}\rangle\langle\xi_{\alpha}|+|\lambda\rangle\langle\lambda|=\bm{I}$,
where $\bm{I}$ is the identity matrix.

Its Drazin inverse is then given by

$\bm{R}^{+}=\sum_{\xi_{\alpha}\neq
0}\frac{1}{\xi_{\alpha}}|\xi_{\alpha}\rangle\langle\xi_{\alpha}|.$ (44)

The Drazin inverse has the following properties [2, 40]:

  1. 1.

$\langle\lambda|\bm{R}^{+}=\bm{0}$ because
$\langle\lambda|\xi_{\alpha}\rangle=0$.

  2. 2.

Similarly, $\bm{R}^{+}|\lambda\rangle=\mathbf{0}$.

  3. 3.

The biorthogonality and completeness relations between eigenvectors imply

$\displaystyle\bm{R}^{+}\bm{R}$ $\displaystyle=\sum_{\xi_{\alpha}\neq
0,\xi_{\beta}\neq
0}\frac{\xi_{\beta}}{\xi_{\alpha}}|\xi_{\alpha}\rangle\langle\xi_{\alpha}|\xi_{\beta}\rangle\langle\xi_{\beta}|$
(45) $\displaystyle=\sum_{\xi_{\alpha}\neq
0}|\xi_{\alpha}\rangle\langle\xi_{\alpha}|$
$\displaystyle=\bm{I}-|\lambda\rangle\langle\lambda|.$

Now we prove that the boundary terms do not contribute to
$\mathcal{K}^{\prime}$. The derivative of the boundary terms with respect to
$\chi$ is

$\displaystyle\left.\frac{\partial\big{[}\ln\left\langle\lambda_{\chi}(0)|\lambda(0)\right\rangle+\ln\left\langle\lambda|\lambda_{\chi}(\tau)\right\rangle\big{]}}{\partial\chi}\right|_{\chi=0}$
(46) $\displaystyle=$
$\displaystyle\left\langle\lambda^{\prime}(0)|\lambda(0)\right\rangle+\left\langle\lambda|\lambda^{\prime}(\tau)\right\rangle$

By substituting Eq. (20) into
$\left\langle\lambda^{\prime}(0)|\lambda(0)\right\rangle$ and making use of
$\bm{R}^{+}|\lambda\rangle=\mathbf{0}$, we see that the first term vanishes.
Similarly, the second term also vanishes.

The explicit expression for the dissipated work, Eq. (26), can be derived by
expressing the vectors and matrices in terms of entries,

$\langle\lambda|\dot{\bm{H}}\bm{R}^{+}|\dot{\lambda}\rangle=\sum_{m,n}\dot{\epsilon}_{m}R^{+}_{mn}\big{(}-\frac{\beta
e^{-\beta\epsilon_{n}}\dot{\epsilon}_{n}}{Z}+\frac{\beta
e^{-\beta\epsilon_{n}}}{Z^{2}}\sum\limits_{l}e^{-\beta\epsilon_{l}}\dot{\epsilon}_{l}\big{)}.$
(47)

The second term vanishes due to the second property of Drazin inverse. By
symmetrizing the first term and plugging it into Eq. (25), we obtain Eqs. (26)
and (27).

Finally, let us prove that, up to the first order of $\tau^{-1}$, the second
derivatives of terms other than the dynamical phase do not contribute to
$\mathcal{K}^{\prime\prime}$. Evaluating Eq. (28) at $\chi=0$ gives

$\langle\lambda^{\prime\prime}|\bm{R}+2\langle\lambda^{\prime}|\dot{\bm{H}}=\lambda^{\prime\prime}\langle\lambda|+2\lambda^{\prime}\langle\lambda^{\prime}|.$
(48)

Multiplying both sides by $\bm{R}^{+}$ we have

$\langle\lambda^{\prime\prime}|(\bm{I}-|\lambda\rangle\langle\lambda|)+2\langle\lambda^{\prime}|\dot{\bm{H}}\bm{R}^{+}=2\lambda^{\prime}\langle\lambda^{\prime}|\bm{R}^{+},$
(49)

and by intermediate normalization
$\langle\lambda^{\prime\prime}|\lambda\rangle=0$, then
$\langle\lambda^{\prime\prime}|$ is given by

$\langle\lambda^{\prime\prime}|=2\lambda^{\prime}\langle\lambda^{\prime}|\bm{R}^{+}-2\langle\lambda^{\prime}|\dot{\bm{H}}\bm{R}^{+}.$
(50)

Recall that $\lambda^{\prime},|\lambda^{\prime}\rangle$ an
$\langle\lambda^{\prime}|$ are all on the order of
$\dot{\bm{H}}={O}(\tau^{-1})$, then
$\langle\lambda^{\prime\prime}|=O(\tau^{-2})$. Similar argument applies to
$|\lambda^{\prime\prime}\rangle$, too. Taking second derivatives of the last
two terms (geometric phase and boundary terms) in Eq. (17a) with respect to
$i\chi$ and evaluating at $\chi=0$, one can see that all terms are of order
$O(\tau^{-2})$.

## References

  * Cavina, Mari, and Giovannetti [2017] V. Cavina, A. Mari, and V. Giovannetti, “Slow Dynamics and Thermodynamics of Open Quantum Systems,” Phys. Rev. Lett. 119, 050601 (2017).
  * Mandal and Jarzynski [2016] D. Mandal and C. Jarzynski, “Analysis of slow transitions between nonequilibrium steady states,” J. Stat. Mech. Theory Exp. 2016, 063204 (2016).
  * Sivak and Crooks [2012] D. A. Sivak and G. E. Crooks, “Thermodynamic Metrics and Optimal Paths,” Phys. Rev. Lett. 108, 190602 (2012).
  * Zulkowski  _et al._ [2012] P. R. Zulkowski, D. A. Sivak, G. E. Crooks, and M. R. DeWeese, “Geometry of thermodynamic control,” Phys. Rev. E 86, 041148 (2012).
  * Zwanzig [2001] R. Zwanzig,  _Nonequilibrium statistical mechanics_ (Oxford university press, 2001).
  * Crooks [2007] G. E. Crooks, “Measuring Thermodynamic Length,” Phys. Rev. Lett. 99, 100602 (2007).
  * Zulkowski and DeWeese [2015] P. R. Zulkowski and M. R. DeWeese, “Optimal protocols for slowly driven quantum systems,” Phys. Rev. E 92, 032113 (2015).
  * Miller  _et al._ [2019] H. J. Miller, M. Scandi, J. Anders, and M. Perarnau-Llobet, “Work Fluctuations in Slow Processes: Quantum Signatures and Optimal Control,” Phys. Rev. Lett. 123, 230603 (2019).
  * Scandi and Perarnau-Llobet [2019] M. Scandi and M. Perarnau-Llobet, “Thermodynamic length in open quantum systems,” Quantum 3, 197 (2019).
  * Scandi  _et al._ [2020] M. Scandi, H. J. D. Miller, J. Anders, and M. Perarnau-Llobet, “Quantum work statistics close to equilibrium,” Phys. Rev. Res. 2, 023377 (2020).
  * Miller and Mehboudi [2020] H. J. Miller and M. Mehboudi, “Geometry of Work Fluctuations versus Efficiency in Microscopic Thermal Machines,” Phys. Rev. Lett. 125, 260602 (2020).
  * Abiuso  _et al._ [2020] P. Abiuso, H. J. D. Miller, M. Perarnau-Llobet, and M. Scandi, “Geometric Optimisation of Quantum Thermodynamic Processes,” Entropy 22, 1076 (2020).
  * Brandner and Saito [2020] K. Brandner and K. Saito, “Thermodynamic Geometry of Microscopic Heat Engines,” Phys. Rev. Lett. 124, 040602 (2020).
  * Frim and DeWeese [2022] A. G. Frim and M. R. DeWeese, “Geometric Bound on the Efficiency of Irreversible Thermodynamic Cycles,” Phys. Rev. Lett. 128, 230601 (2022).
  * Kirkwood [1946] J. G. Kirkwood, “The Statistical Mechanical Theory of Transport Processes I. General Theory,” J. Chem. Phys. 14, 180–201 (1946).
  * Weinhold [1975] F. Weinhold, “Metric geometry of equilibrium thermodynamics,” J. Chem. Phys. 63, 2479–2483 (1975).
  * Ruppeiner [1979] G. Ruppeiner, “Thermodynamics: A Riemannian geometric model,” Phys. Rev. A 20, 1608–1613 (1979).
  * Salamon and Berry [1983] P. Salamon and R. S. Berry, “Thermodynamic length and dissipated availability,” Phys. Rev. Lett. 51, 1127–1130 (1983).
  * Janyszek and Mrugal/a [1989] H. Janyszek and R. Mrugal/a, “Riemannian geometry and the thermodynamics of model magnetic systems,” Phys. Rev. A 39, 6515–6523 (1989).
  * Brody and Rivier [1995] D. Brody and N. Rivier, “Geometrical aspects of statistical mechanics,” Phys. Rev. E 51, 1006–1011 (1995).
  * Ruppeiner [1995] G. Ruppeiner, “Riemannian geometry in thermodynamic fluctuation theory,” Rev. Mod. Phys. 67, 605–659 (1995).
  * Seifert [2008] U. Seifert, “Stochastic thermodynamics: Principles and perspectives,” Eur. Phys. J. B 64, 423–431 (2008).
  * Seifert [2012] U. Seifert, “Stochastic thermodynamics, fluctuation theorems and molecular machines.” Rep. Prog. Phys. 75, 126001 (2012).
  * Jarzynski [1997a] C. Jarzynski, “Nonequilibrium Equality for Free Energy Differences,” Phys. Rev. Lett. 78, 2690–2693 (1997a).
  * Jarzynski [1997b] C. Jarzynski, “Equilibrium free-energy differences from nonequilibrium measurements: A master-equation approach,” Phys. Rev. E 56, 5018–5035 (1997b).
  * Speck and Seifert [2004] T. Speck and U. Seifert, “Distribution of work in isothermal nonequilibrium processes,” Phys. Rev. E 70, 066112 (2004).
  * Sinitsyn and Nemenman [2007] N. A. Sinitsyn and I. Nemenman, “The Berry phase and the pump flux in stochastic chemical kinetics,” EPL 77, 58001 (2007).
  * Sinitsyn [2009] N. A. Sinitsyn, “The stochastic pump effect and geometric phases in dissipative and stochastic systems,” J. Phys. A 42, 193001 (2009).
  * Sagawa and Hayakawa [2011] T. Sagawa and H. Hayakawa, “Geometrical expression of excess entropy production,” Phys. Rev. E 84, 051110 (2011).
  * Ren, Hänggi, and Li [2010] J. Ren, P. Hänggi, and B. Li, “Berry-Phase-Induced Heat Pumping and Its Impact on the Fluctuation Theorem,” Phys. Rev. Lett. 104, 170601 (2010).
  * Yuge  _et al._ [2012] T. Yuge, T. Sagawa, A. Sugita, and H. Hayakawa, “Geometrical pumping in quantum transport: Quantum master equation approach,” Phys. Rev. B 86, 235308 (2012).
  * Goswami, Agarwalla, and Harbola [2016] H. P. Goswami, B. K. Agarwalla, and U. Harbola, “Geometric effects in nonequilibrium electron transfer statistics in adiabatically driven quantum junctions,” Phys. Rev. B 93, 195441 (2016).
  * Gu  _et al._ [2018] J. Gu, X.-G. Li, H.-P. Cheng, and X.-G. Zhang, “Adiabatic Spin Pump through a Molecular Antiferromagnet Ce 3 Mn 8 III,” J. Phys. Chem. C 122, 1422–1429 (2018).
  * Giri and Goswami [2017] S. K. Giri and H. P. Goswami, “Geometric phaselike effects in a quantum heat engine,” Phys. Rev. E 96, 052129 (2017).
  * Berry [1984] M. V. Berry, “Quantal Phase Factors Accompanying Adiabatic Changes,” Proc. R. Soc. A 392, 45–57 (1984).
  * Van Kampen [1992] N. G. Van Kampen,  _Stochastic processes in physics and chemistry_ , Vol. 1 (Elsevier, 1992).
  * Breuer and Petruccione [2002] H.-P. Breuer and F. Petruccione,  _The theory of open quantum systems_ (Oxford University Press, 2002).
  * Verley, Van den Broeck, and Esposito [2013] G. Verley, C. Van den Broeck, and M. Esposito, “Modulated two-level system: Exact work statistics,” Phys. Rev. E 88, 032137 (2013).
  * Ohkubo and Eggel [2010] J. Ohkubo and T. Eggel, “Noncyclic and nonadiabatic geometric phase for counting statistics,” J. Phys. A 43, 425001 (2010).
  * Crooks [2018] G. E. Crooks, “On the Drazin inverse of the rate matrix,”  (2018).
  * Shiraishi, Funo, and Saito [2018] N. Shiraishi, K. Funo, and K. Saito, “Speed Limit for Classical Stochastic Processes,” Phys. Rev. Lett. 121, 070601 (2018).
  * Sinitsyn, Akimov, and Chernyak [2011] N. A. Sinitsyn, A. Akimov, and V. Y. Chernyak, “Supersymmetry and fluctuation relations for currents in closed networks,” Phys. Rev. E 83, 021107 (2011).
  * Griffiths [2004] D. J. Griffiths,  _Introduction to Quantum Mechanics_ , 2nd ed. (Pearson Prentice Hall, 2004).
  * Switkes  _et al._ [1999] M. Switkes, C. M. Marcus, K. Campman, and A. C. Gossard, “An Adiabatic Quantum Electron Pump,” Science 283, 1905–1908 (1999).
  * Landauer [1961] R. Landauer, “Irreversibility and Heat Generation in the Computing Process,” IBM J. Res. Dev 5, 183–191 (1961).
  * Ma  _et al._ [2022] Y.-H. Ma, J.-F. Chen, C. P. Sun, and H. Dong, “Minimal energy cost to initialize a bit with tolerable error,” Phys. Rev. E 106, 034112 (2022).
  * Brody [2014] D. C. Brody, “Biorthogonal Quantum Mechanics,” J. Phys. A 47, 035305 (2014).


∎

11institutetext: Ya-Nan Wang22institutetext: School of Mathematical Sciences,
Nanjing Normal University, Nanjing, 210097, China
22email<EMAIL_ADDRESS>33institutetext: Jun Yan 44institutetext:
School of Mathematical Sciences, Fudan University & Shanghai Key Laboratory
for Contemporary Applied Mathematics, Shanghai 200433, China
44email<EMAIL_ADDRESS>55institutetext: Jianlu Zhang 66institutetext:
Hua Loo-Keng Key Laboratory of Mathematics & Mathematics Institute, Academy of
Mathematics and systems science, Chinese Academy of Sciences, Beijing 100190,
China
66email<EMAIL_ADDRESS>
# Parameterized viscosity solutions of convex Hamiltonian systems with time
periodic damping

Ya-Nan Wang  Jun Yan  Jianlu Zhang

(Received: date / Accepted: date)

###### Abstract

In this article we develop an analogue of Aubry Mather theory for time
periodic dissipative equation

$\left\\{\begin{aligned} \dot{x}&=\partial_{p}H(x,p,t),\\\
\dot{p}&=-\partial_{x}H(x,p,t)-f(t)p\end{aligned}\right.$

with $(x,p,t)\in T^{*}M\times\mathbb{T}$ (compact manifold $M$ without
boundary). We discuss the asymptotic behaviors of viscosity solutions of
associated Hamilton-Jacobi equation

$\partial_{t}u+f(t)u+H(x,\partial_{x}u,t)=0,\quad(x,t)\in M\times\mathbb{T}$

w.r.t. certain parameters, and analyze the meanings in controlling the global
dynamics. We also discuss the prospect of applying our conclusions to many
physical models.

###### Keywords:

viscosity solution, weak KAM solution, convex Hamiltonian, Aubry Mather
theory, global attractor, rotation number

###### MSC:

37J50, 37J55, 35B40,49L25

##  1 Introduction

For a smooth compact Riemannian manifold $M$ without boundary, the Hamiltonian
$H$ is usually characterized as a $C^{r\geq 2}-$smooth function on the
cotangent bundle $T^{*}M$, with the associated Hamilton equation defined by

$\displaystyle{\sf(Conservative)\quad}\left\\{\begin{aligned}
\dot{x}&=\partial_{p}H(x,p)\\\
\dot{p}&=-\partial_{x}H(x,p)\end{aligned}\right.$ (1)

for each initial point $(x,p)\in T^{*}M$. From the physical aspect, the
Hamiltonian equation describes the movement of particles with conservative
energy, since the Hamiltonian $H(x,p)$ verifies to be a First Integral of (1).
In particular, if the potential periodically depends on the time $t$ (for
systems with periodic propulsion or procession), we can introduce an augmented
Hamiltonian

$\displaystyle\widetilde{H}(x,p,t,I)=I+H(x,p,t),\quad\quad(x,p,t,I)\in
T^{*}M\times T^{*}\mathbb{T}$ (2)

such that the associated Hamiltonian equation

$\displaystyle{\sf(Conservative)\quad}\left\\{\begin{aligned}
\dot{x}&=\partial_{p}H(x,p,t)\\\ \dot{p}&=-\partial_{x}H(x,p,t)\\\
\dot{t}&=1\\\ \dot{I}&=-\partial_{t}H(x,p,t)\end{aligned}\right.$ (3)

still preserves $\widetilde{H}$.

However, the realistic motion of the masses inevitably sustains a dissipation
of energy, due to the friction from the environment, e.g. the wind, the fluid,
interface etc. That urges us to make rational modification of previous
equations. In the current paper, the damping is assumed to be time-
periodically proportional to the momentum. Precisely, we modify (3) into

$\displaystyle{\sf(Dissipative)\quad}\left\\{\begin{aligned}
\dot{x}&=\partial_{p}H(x,p,t)\\\ \dot{p}&=-\partial_{x}H(x,p,t)-f(t)p\\\
\dot{t}&=1\\\ \dot{I}&=-\partial_{t}H(x,p,t)-f^{\prime}(t)u-f(t)I\\\
\dot{u}&=\langle H_{p},p\rangle-H+\alpha-f(t)u\end{aligned}\right.$ (4)

with $\alpha\in\mathbb{R}$ being a constant of initial energy and $f\in
C^{r\geq 2}(\mathbb{T}:=\mathbb{R}/[0,1],\mathbb{R})$. Notice that the former
three equations of (19) is decoupled with the latter two, so we can denote the
flow of the former three equations in (4) by $\varphi_{H}^{t}$ and by
$\widehat{\varphi}_{H}^{t}$ the flow of the whole (4). The following
individual cases of $f(t)$ will be considered:

  * •

(H0-) $[f]:=\int_{0}^{1}f(t)dt>0$

  * •

(H0+) $[f]<0$

  * •

(H00) $[f]=0$

Besides, we propose the following standing assumptions for the Hamiltonian:

  * •

(H1) [Smoothness] $H:TM\times\mathbb{T}\rightarrow\mathbb{R}$ is $C^{r\geq 2}$
smooth;

  * •

(H2) [Convexity] For any $(x,t)\in M\times\mathbb{T}$, $H(x,\cdot,t)$ is
strictly convex on $T_{x}^{*}M$;

  * •

(H3) [Superlinearity] For any $(x,t)\in M\times\mathbb{T}$,
$\lim_{|p|_{x}\rightarrow+\infty}H(x,p,t)/|p|_{x}=+\infty$ where $|\cdot|_{x}$
is the norm deduced from the Riemannian metric.

  * •

(H4) [Completeness] For any $(x,p,\theta)\in T^{*}M\times\mathbb{T}$, the flow
$\varphi_{H}^{t}(x,p,\theta)$ exists for all $t\in\mathbb{R}$.

###### Remark 1.1

  * i)

As we can see, the three different cases of (H0) respectively leads to a
dissipation, acceleration and periodic conservation of energy along
$\widehat{\varphi}_{H}^{t}$ in the forward time, if we take

$\displaystyle\widehat{H}(x,p,t,I,u)=\widetilde{H}(x,p,t,I)+f(t)u-\alpha.$ (5)

This is because $\dfrac{d}{dt}\widehat{H}=-f(t)\widehat{H}$.

  * •

(H1-H3) are usually called Tonelli conditions. As for (H4), the completeness
of $\varphi_{H}^{t}$ is actually equivalent to the completeness of
$\widehat{\varphi}_{H}^{t}$. A sufficient condition to (H4) is the following:

$|H_{x}|\leq\kappa(1+|p|_{x})\text{ for all }(x,p,t)\in
T^{*}M\times\mathbb{T}$

for some constant $\kappa$.

  * ii)

Observe that the time-1 map
$\varphi_{H}^{1}:\\{(x,p,t=0)\\}\rightarrow\\{(x,p,t=0)\\}$ is conformally
symplectic, i.e.

$(\varphi_{H}^{1})^{*}dp\wedge dx=e^{[f]}dp\wedge dx.$

Such maps have wide applications in astronomy CC, optimal transport WL,
biological physics Ca2 and economics B etc (see Sec. 2 for more details).

In the following, we will explore the dynamics of (4) by using an analogue of
Aubry Mather theory Mat or weak KAM theory Fa. Similar research was exploited
in CCJW; WWY2; WWY3, for generalized $1^{st}$ order PDEs. The current paper
have a lot of similarities in the methodology with these works.

###  1.1 Variational Principle and Hamilton Jacobi equation

As the dual of the Hamiltonian, the Lagrangian can be defined as the following

$\displaystyle L(x,v,t):=\max_{p\in T^{*}_{x}M}\langle
p,v\rangle-H(x,p,t),\quad(x,v,t)\in TM\times\mathbb{T}.$ (6)

of which the maximum is achieved for $v=H_{p}(x,p,t)\in T_{x}M$, once (H1-H3)
are assumed. Therefore, the Legendre transformation

$\mathcal{L}:T^{*}M\times\mathbb{T}\rightarrow
TM\times\mathbb{T},\quad\text{via }(x,p,t)\rightarrow(x,H_{p}(x,p,t),t)$ (7)

is a diffeomorphism. Notice that the Lagrangian
$L:TM\times\mathbb{T}\rightarrow\mathbb{R}$ is also $C^{r}-$smooth and convex,
superlinear in $p\in T_{x}M$, so by a slight abusion of notions we say it
satisfies (H1-H3) as well. As a conjugation of $\varphi^{t}_{H}$, the Euler-
Lagrangian flow $\varphi_{L}^{t}$ is defined by

$\left\\{\begin{aligned} &\dot{x}=v,\\\
&\frac{d}{dt}L_{v}(x,v,t)=L_{x}(x,v,t)-f(t)L_{v}(x,v,t),\\\
&\dot{t}=1.\end{aligned}\right.$ (E-L)

It has equivalently effective in exploring the dynamics of (4), and the
completeness of $\varphi_{L}^{t}$ is equivalent to the completeness of
$\varphi_{L}^{t}$. To present (E-L) a variational characterization, we
introduce a minimal variation on the absolutely continuous curves with fixed
endpoints

$h_{\alpha}^{s,t}(x,y)=\inf_{\begin{subarray}{c}\gamma\in C^{ac}([s,t],M)\\\
\gamma(s)=x,\gamma(t)=y\end{subarray}}\int^{t}_{s}e^{F(\tau)}(L(\gamma,\dot{\gamma},\tau)+\alpha)\mbox{d}\tau,$

where $F(t)=\int^{t}_{0}f(\tau)\mbox{d}\tau$ and $\alpha\in\mathbb{R}$. It is
a classical result in the calculus of variations, the infimum is always
available for all $s<t\in\mathbb{R}$, which is actually $C^{r}-$ smooth and
satisfies (E-L), once (H4) is assumed (due to the Weierstrass Theorem in Mat
or Theorem 3.7.1 of Fa).

######  Theorem 1.2 (main 1)

For $f(t)$ satisfying (H0-), $H(x,p,t)$ satisfying (H1-H4) and any
$\alpha\in\mathbb{R}$, the following

$u_{\alpha}^{-}(x,t):=\inf_{\begin{subarray}{c}\gamma\in
C^{ac}((-\infty,t],M)\\\
\gamma(t)=x\end{subarray}}\int^{t}_{-\infty}e^{F(s)-F(t)}(L(\gamma(s),\dot{\gamma}(s),s)+\alpha)\mbox{d}s$

is well defined for $(x,t)\in M\times\mathbb{R}$ and satisfies

  1. 1.

(Periodicity) $u_{\alpha}^{-}(x,t+1)=u_{\alpha}^{-}(x,t)$ for any $x\in M$ and
$t\in\mathbb{R}$. By taking $\bar{t}\in[0,1)$ with $t\equiv\bar{t}\ (mod\ 1)$
for any $t\in\mathbb{R}$, we can interpreted $u_{\alpha}^{-}$ as a function on
$M\times\mathbb{T}$.

  2. 2.

(Lipschitzness) $u_{\alpha}^{-}:M\times\mathbb{T}\rightarrow\mathbb{R}$ is
Lipschitz, with the Lipschitz constant depending on $L$ and $f$;

  3. 3.

(Domination111Any function $\omega\in C(M\times\mathbb{T},\mathbb{R})$
satisfying (8) is called a (viscosity) subsolution of (HJ+) and denoted by
$\omega\prec_{f}L+\alpha$. ) For any absolutely continuous curve
$\gamma:[s,t]\rightarrow M$ connecting $(x,\bar{s})\in M\times\mathbb{T}$ and
$(y,\bar{t})\in M\times\mathbb{T}$, we have

$\displaystyle
e^{F(t)}u_{\alpha}^{-}(y,\bar{t})-e^{F(s)}u_{\alpha}^{-}(x,\bar{s})\leq\int_{s}^{t}e^{F(\tau)}\Big{(}L(\gamma,\dot{\gamma},\tau)+\alpha\Big{)}d\tau.$
(8)

  4. 4.

(Calibration) For any $(x,\theta)\in M\times\mathbb{T}$, there exists a
backward calibrated curve curve
$\gamma_{x,\theta}^{-}:(-\infty,\theta]\rightarrow M$, $C^{r}-$smooth and
ending with $\gamma_{x,\theta}^{-}(\theta)=x$, such that for all $s\leq
t\leq\theta$, we have

$\displaystyle
e^{F(t)}u_{\alpha}^{-}(\gamma_{x,\theta}^{-}(t),\bar{t})-e^{F(s)}u_{\alpha}^{-}(\gamma_{x,\theta}^{-}(s),\bar{s})$
(9) $\displaystyle=$
$\displaystyle\int_{s}^{t}e^{F(\tau)}\Big{(}L(\gamma_{x,\theta}^{-},\dot{\gamma}_{x,\theta}^{-},\tau)+\alpha\Big{)}d\tau.$

  5. 5.

(Viscosity) $u_{\alpha}^{-}:M\times\mathbb{T}\rightarrow\mathbb{R}$ is a
viscosity solution of the following Stationary Hamilton-Jacobi equation (with
time periodic damping):

$\partial_{t}u+f(t)u+H(x,\partial_{x}u,t)=\alpha,\quad(x,t)\in
M\times\mathbb{T},\ \alpha\in\mathbb{R}.$ (HJ+)

######  Theorem 1.3 (main 1’)

For $f(t)$ satisfying (H00) and $H(x,p,t)$ satisfying (H1-H4), there exists a
unique $c(H)\in\mathbb{R}$ (Mañé Critical Value) such that

$u^{-}_{z,\bar{\varsigma}}(x,\bar{t}):=\varliminf_{\begin{subarray}{c}\bar{\varsigma}\equiv\varsigma,\bar{t}\equiv
t(mod\;1)\\\
t-\varsigma\to+\infty\end{subarray}}\bigg{(}\inf_{\begin{subarray}{c}\gamma\in
C^{ac}([\varsigma,t],M)\\\
\gamma(\varsigma)=z,\gamma(t)=x\end{subarray}}\int^{t}_{\varsigma}e^{F(\tau)-F(t)}\big{(}L(\gamma,\dot{\gamma},\tau)+c(H)\big{)}\mbox{d}\tau\bigg{)}$
(10)

is well defined on $M\times\mathbb{T}$ (for any fixed $(z,\bar{\varsigma})\in
M\times\mathbb{T}$) and satisfies

  1. 1.

(Lipschitzness)
$u_{z,\bar{\varsigma}}^{-}:M\times\mathbb{T}\rightarrow\mathbb{R}$ is
Lipschitz.

  2. 2.

(Domination) For any Lipschitz continuous curve $\gamma:[s,t]\rightarrow M$
connecting $(x,\bar{s})\in M\times\mathbb{T}$ and $(y,\bar{t})\in
M\times\mathbb{T}$, we have

$\displaystyle
e^{F(t)}u_{z,\bar{\varsigma}}^{-}(y,\bar{t})-e^{F(s)}u_{z,\bar{\varsigma}}^{-}(x,\bar{s})\leq\int_{s}^{t}e^{F(\tau)}\Big{(}L(\gamma,\dot{\gamma},\tau)+c(H)\Big{)}d\tau.$
(11)

Namely, $u_{z,\bar{\varsigma}}^{-}\prec_{f}L+c(H)$.

  3. 3.

(Calibration) For any $(x,\theta)\in M\times\mathbb{T}$, there exists a
$C^{r}$ curve $\gamma_{x,\theta}^{-}:(-\infty,\theta]\rightarrow M$ with
$\gamma_{x,\theta}^{-}(\theta)=x$, such that for all $s\leq t\leq\theta$, we
have

$\displaystyle
e^{F(t)}u_{z,\bar{\varsigma}}^{-}(\gamma_{x,\theta}^{-}(t),\bar{t})-e^{F(s)}u_{z,\bar{\varsigma}}^{-}(\gamma_{x,\theta}^{-}(s),\bar{s})$
(12) $\displaystyle=$
$\displaystyle\int_{s}^{t}e^{F(\tau)}\Big{(}L(\gamma_{x,\theta}^{-},\dot{\gamma}_{x,\theta}^{-},\tau)+c(H)\Big{)}d\tau.$

  4. 4.

(Viscosity) $u_{z,\bar{\varsigma}}^{-}$ is a viscosity solution of

$\partial_{t}u+f(t)u+H(x,\partial_{x}u,t)=c(H),\quad(x,t)\in
M\times\mathbb{T}.$ (HJ0)

Following the terminologies in Fa; MS, it’s appropriate to call the function
given in Theorem 1.2 (resp. Theorem 1.3) a weak KAM solution. Such a solution
can be used to pick up different types of invariant sets with variational
meanings of (4):

######  Theorem 1.4 (main 2)

For $f(t)$ satisfying (H0-), $H(x,p,t)$ satisfying (H1-H4) and any
$\alpha\in\mathbb{R}$, we can get the following sets:

  * •

(Aubry Set) $\gamma:\mathbb{R}\rightarrow M$ is called globally calibrated, if
for any $s<t\in\mathbb{R}$, (9) holds on $[s,t]$. There exists a
$\varphi_{L}^{t}-$invariant set defined by

$\widetilde{\mathcal{A}}:=\\{(\gamma(t),\dot{\gamma}(t),\bar{t})\in
TM\times\mathbb{T}|\gamma\text{ is globally calibrated}\\}$

with the following properties:

    * –

$\widetilde{\mathcal{A}}$ is a Lipschitz graph over the projected Aubry set
$\mathcal{A}:=\pi\widetilde{\mathcal{A}}\subset M\times\mathbb{T}$, where
$\pi:T^{*}M\times\mathbb{T}\rightarrow M\times\mathbb{T}$ is the standard
projection.

    * –

$\widetilde{\mathcal{A}}$ is upper semicontinuous w.r.t.
$L:TM\times\mathbb{T}\rightarrow\mathbb{R}$

    * –

$u_{\alpha}^{-}$ is differentiable on $\mathcal{A}$.

  * •

(Mather Set) Suppose $\mathfrak{M}_{L}$ is the set of all
$\varphi_{L}^{t}-$invariant probability measure, then
$\tilde{\mu}\in\mathfrak{M}_{L}$ is called a Mather measure if it minimizes

$\min_{\tilde{\nu}\in\mathfrak{M}_{L}}\int_{TM\times\mathbb{T}}L+\alpha-f(t)u_{\alpha}^{-}\mbox{d}\tilde{\nu}.$

Let’s denote by $\mathfrak{M}_{m}$ the set of all Mather measures.
Accordingly, the Mather set is defined by

$\widetilde{\mathcal{M}}:=\overline{\bigcup\\{supp\
\tilde{\mu}|\tilde{\mu}\in\mathfrak{M}_{m}\\}}$

which satisfies

    1. 1.

$\widetilde{\mathcal{M}}\neq\emptyset$ and
$\widetilde{\mathcal{M}}\subset\widetilde{\mathcal{A}}$.

    2. 2.

$\widetilde{\mathcal{M}}$ is a Lipschitz graph over the projected Mather set
$\mathcal{M}:=\pi\widetilde{\mathcal{M}}\subset M\times\mathbb{T}$.

  * •

(Maximal Global Attractor) Define

$\displaystyle\widehat{\Sigma}_{H}^{-}$ $\displaystyle:=$
$\displaystyle\big{\\{}(x,p,\bar{s},\alpha-f(s)u-H(x,p,s),u)\in T^{*}M\times
T^{*}\mathbb{T}\times\mathbb{R}\big{|}$ $\displaystyle\quad
u>u_{\alpha}^{-}(x,s)\big{\\}}$

and

$\displaystyle\widehat{\Sigma}_{H}^{0}$ $\displaystyle:=$
$\displaystyle\big{\\{}(x,p,\bar{s},\alpha-f(s)u-H(x,p,s),u)\in T^{*}M\times
T^{*}\mathbb{T}\times\mathbb{R}\big{|}$ $\displaystyle\quad
u=u_{\alpha}^{-}(x,s)\big{\\}},$

then $\Omega:=\bigcap_{t\geq
0}\widehat{\varphi}_{H}^{t}(\widehat{\Sigma}_{H}^{-}\cup\widehat{\Sigma}_{H}^{0})$
is the maximal $\widehat{\varphi}_{H}^{t}-$invariant set, which satisfies:

    1. 1.

If the $p-$component of $\Omega$ is bounded, then the $u-$ and $I-$component
of $\Omega$ are also bounded.

    2. 2.

If $\Omega$ is compact, it has to be a global attractor in the sense that for
any point $(x,p,\bar{s},I,u)\in T^{*}M\times T^{*}\mathbb{T}\times\mathbb{R}$
and any open neighborhood $\mathcal{U}\supseteq\Omega$, there exists a
$T_{\Omega}(\mathcal{U})$ such that for all $t\geq T_{\Omega}(\mathcal{U})$,
$\widehat{\varphi}_{H}^{t}(x,p,\bar{s})\in\mathcal{U}$. Besides, the
followings hold:

       * –

$\Omega$ is a maximal attractor set, i.e. it isn’t strictly contained in any
other global attractor;

       * –

$\widehat{\mathcal{A}}$ is the maximal invariant set contained in
$\widehat{\Sigma}_{H}^{0}$, where

$\displaystyle\widehat{\mathcal{A}}$ $\displaystyle:=$
$\displaystyle\Big{\\{}\Big{(}\mathcal{L}(x,\partial_{x}u_{\alpha}^{-}(x,s),\bar{s}),\partial_{t}u_{\alpha}^{-}(x,s),u_{\alpha}^{-}(x,s)\Big{)}\in
TM\times$ $\displaystyle\quad
T^{*}\mathbb{T}\times\mathbb{R}\Big{|}(x,\bar{s})\in\mathcal{A}\Big{\\}}.$

######  Theorem 1.5 (main 2’)

For $f(t)$ satisfying (H00) and $H(x,p,t)$ satisfying (H1-H4), the Mañé
Critical Value $c(H)$ has an alternative expression

$\displaystyle-c(H)=\dfrac{\inf_{\tilde{\mu}\in\mathfrak{M}_{L}}\int_{TM\times\mathbb{T}}e^{F(t)}L(x,v,t)\mbox{d}\tilde{\mu}}{\int_{0}^{1}e^{F(t)}\mbox{d}t}.$
(13)

Moreover, the minimizer achieving the right side of (13) has to be a Mather
measure. Similarly we can define the Mather set $\widetilde{\mathcal{M}}$ as
the union of the support sets of all the Mather measures, which is Lipschitz-
graphic over the projected Mather set
$\mathcal{M}:=\pi\widetilde{\mathcal{M}}$.

###  1.2 Parametrized viscosity solutions and asymptotic dynamics

In this section we deal with two kinds of parametrized viscosity solutions
with practical meanings. The first case corresponds to a Hamiltonian

$\displaystyle\widehat{H}_{\delta}(x,p,t,I,u):=I+H(x,p,t)+f_{\delta}(t)u,$
(14)

with $(x,p,\bar{t},I,u)\in T^{*}M\times T^{*}\mathbb{T}\times\mathbb{R}$ and
$f_{\delta}\in C^{r}(\mathbb{T},\mathbb{R})$ continuous of
$\delta\in\mathbb{R}$. For suitable $\alpha\in\mathbb{R}$, we can seek the
weak KAM solution of

$\displaystyle\partial_{t}u_{\delta}(x,t)+H(x,\partial_{x}u_{\delta},t)+f_{\delta}(t)u_{\delta}=\alpha$
(15)

as we did in previous theorems. Consequently, it’s natural to explore the
convergence of viscosity solutions w.r.t. the parameter $\delta$:

######  Theorem 1.6 (main 3)

Suppose $f_{\delta}$ converges to $f_{0}$ w.r.t. the uniform norm as
$\delta\rightarrow 0_{+}$ such that $[f_{0}]=0$ and the right derivative of
$f_{\delta}$ w.r.t. $\delta$ exists at $0$, i.e.

$\displaystyle f_{1}(t):=\lim_{\delta\rightarrow
0_{+}}\frac{f_{\delta}(t)-f_{0}(t)}{\delta}>0.$ (16)

If $H(x,p,t)$ satisfying (H1-H4), then there exists a unique
$c(H)\in\mathbb{R}$ given by (13) and a $\delta_{0}>0$, such that the weak KAM
solution $u^{-}_{\delta}(x,t)$ of (15) associated with $f_{\delta}$ and
$\alpha_{\delta}\equiv c(H)$ for all $\delta\in(0,\delta_{0}]$ converges to a
uniquely identified viscosity solution of

$\displaystyle\partial_{t}u(x,t)+H(x,\partial_{x}u,t)+f_{0}(t)u=c(H),$ (17)

which equals

$\sup\Big{\\{}u\prec_{f_{0}}L+c(H)\Big{|}\int_{TM\times\mathbb{T}}e^{F_{0}(t)}f_{1}(t)\cdot
u(x,t)d\tilde{\mu}\leq 0,\
\forall\;\tilde{\mu}\in\mathfrak{M}_{m}(\delta=0)\Big{\\}}$

with $F_{0}(t)=\int_{0}^{t}f_{0}(\tau)\mbox{d}\tau$ and $\mathfrak{M}_{m}(0)$
being the set of Mather measures for the system with $\delta=0$.

###### Remark 1.7

The convergence of the viscosity solutions $1^{st}$ order PDEs was earlier
discussed in CCIZ; DFIZ; WYZ; Z, where the Comparison principle was used to
guarantee the uniqueness of viscosity solution for (15). However, in our case
we didn’t assume $f_{\delta}$ to be nonnegative, which invalidates this
principle and brings new new difficulties to prove the equi-boundedness and
equi-Lipschitzness of $\\{u_{\delta}^{-}\\}_{\delta>0}$. Fortunately, by
analyzing the properties of the Lax-Oleinik semigroups we can still overcome
these difficulties, see Sec. LABEL:s4 for more details.

The second parametrized problem we concern takes $M=\mathbb{T}$ and a
mechanical $H(x,p,t)$. We can involve a cohomology parameter $c\in
H^{1}(\mathbb{T},\mathbb{R})$ to

$\displaystyle\widehat{H}(x,p,t,I,u)=I+\underbrace{\frac{1}{2}(p+c)^{2}+V(x,t)}_{H(x,p,t)}+f(t)u$
(18)

of which $H(x,p,t)$ surely satisfies (H1-H4), then (4) becomes

$\displaystyle{\sf(Dissipative)\quad}\left\\{\begin{aligned} \dot{x}&=p+c\\\
\dot{p}&=-V_{x}-f(t)p\\\ \dot{t}&=1\\\ \dot{I}&=-V_{t}-f^{\prime}(t)u-f(t)I\\\
\dot{u}&=\frac{1}{2}(p^{2}-c^{2})-V(x,t)-f(t)u.\end{aligned}\right.$ (19)

In physical models, the former three equations of (19) is usually condensed
into a single equation

$\displaystyle\ddot{x}+V_{x}(x,t)+f(t)(\dot{x}-c)=0,\quad(x,t)\in
M\times\mathbb{T}.$ (20)

######  Theorem 1.8 (main 4)

For $f(t)$ satisfying (H0-), the following conclusions hold for equation (20):

  * •

For any $c\in H^{1}(\mathbb{T},\mathbb{R})$, there exists a unified rotation
number of $\widetilde{\mathcal{A}}(c)$, which is defined by

$\rho(c):=\lim_{T\rightarrow+\infty}\frac{1}{T}\int_{0}^{T}d\gamma,\quad\forall\
\text{globally calibrated curve }\gamma.$

  * •

$\rho(c)$ is continuous of $c\in H^{1}(\mathbb{T},\mathbb{R})$; Moreover, we
have

$\displaystyle|\rho(c)-c|\leq\varsigma([f])\cdot\|V(x,t)\|_{C^{1}}$ (21)

for some constant $\varsigma$ depending only on $[f]$. Consequently, for any
$p/q\in\mathbb{Q}$ irreducible, there always exists a $c_{p/q}$ such that
$\rho(c_{p/q})=p/q$.

  * •

There exists an compact maximal global attractor $\Omega\subset
T^{*}\mathbb{T}\times T^{*}\mathbb{T}\times\mathbb{R}$ of the flow
$\widehat{\varphi}_{H}^{t}$.

Organization of the article: The paper is organized as follows: In Sec. 2, we
exhibit a list of physical models with time-periodic damping. For these
models, we state some notable dynamic phenomena and show how these phenomena
can be linked to our main conclusions. In Sec. 3, we prove Theorem 1.2 and
Theorem 1.3. In Sec. LABEL:s3, we get an analogue Aubry Mather theory for
systems satisfying (H0-) condition, and prove Theorem 1.4. Besides, we also
prove Theorem 1.5 for systems satisfying (H00) condition. In Sec. LABEL:s4, we
discuss the parametrized viscosity solutions of (15), and prove the
convergence of them. In Sec. LABEL:s5, for 1-D mechanical systems with time
periodic damping, we prove Theorem 1.8, which is related to the dynamic
phenomena of the models in Sec 2. For the consistency of the proof, parts of
preliminary conclusions are postponed to the Appendix.

Acknowledgements: The first author is supported by Natural Scientific
Foundation of China (Grant No.11501437). The second author is supported by
National Natural Science Foundation of China (Grant No. 11631006, 11790272)
and Shanghai Science and Technology Commission (Grant No. 17XD1400500). The
third author is supported by the Natural Scientific Foundation of China (Grant
No. 11901560). All the authors are grateful to Prof. Wei Cheng for helpful
discussion about the details.

##  2 Zoo of practical models

In this section we display a bunch of physical models with time-periodic
damping, and introduce some practical problems (related with our main
conclusions) around them.

###  2.1 Conformally symplectic systems

For $f(t)\equiv\lambda>0$ being constant, we get a so called conformally
symplectic system (or discount system). The associated ODE becomes

$\displaystyle\left\\{\begin{aligned} \dot{x}&=\partial_{p}H(x,p,t),\\\
\dot{p}&=-\partial_{x}H(x,p,t)-\lambda p.\end{aligned}\right.$ (22)

This kind of systems has been considered in CCD; DFIZ; MS, although earlier
results on Aubry-Mather sets have been discussed by Le Calvez LC and Casdagli
Ca for $M=\mathbb{T}$. Besides, we need to specify that the Duffing equation
with viscous damping also conforms to this case, which concerns all kinds of
oscillations widely found in electromagnetics M and elastomechanics MH.

A significant property this kind of systems possess is that

$(\varphi_{H}^{1})^{*}dp\wedge dx=e^{\lambda}dp\wedge dx.$

When $H(x,p,t)$ is mechanical, the equation usually describes the low velocity
oscillation of a solid in a fluid medium (see Fig. 1), which can be formally
expressed as

$\displaystyle\ddot{x}+\lambda\dot{x}+\partial_{x}V(x,t)=0,\quad
x\in\mathbb{T},\;\lambda>0.$ (23)

Chaos and bifurcations topics of this setting has ever been rather popular in
1970s H.

Figure 1: A dissipative pendulum with $\lambda=1/5$ and $V(t,x)=1-\cos x$.

###  2.2 Tidal torque model

The tidal torque model was firstly introduced by P, describing the motion of a
rigid satellite $S$ under the gravitational influence of a point-mass planet
$P$. Due to the internal non-rigidity of the body, a tidal torque will causes
a time-periodic dissipative to the motion of $S$, which can be formalized by

$\displaystyle\ddot{x}+\varepsilon
V_{x}(x,e,t)+\kappa\eta(e,t)(\dot{x}-c(e))=0,\quad(x,t)\in\mathbb{T}^{2},$
(24)

with the parameter $e$ is the eccentricity of the elliptic motion $S$ around
$P$. Due to the astronomical observation, $\varepsilon$ is the equatorial
ellipticity of the satellite and

$\kappa\propto\frac{1}{a^{3}}\cdot\frac{m_{P}}{m_{S}},$

with $a$ being the semi-major and $m_{P}$ (resp. $m_{S}$) being the mass
respectively.

Figure 2: A tidal torque model for Moon-Earth and Mercury-Sun.

Although this model might seem very special, there are several examples in the
solar system for which such a model yields a good description of the motion,
at least in a first approximation, and anyhow represents a first step toward
the understanding of the problem. For instance, in the pairs Moon-Earth,
Enceladus-Saturn, Dione-Saturn, Rhea-Saturn even Mercury-Sun this model is
available. Besides, we need to specify that usually $\kappa\ll\varepsilon$ in
all these occasions.

A few interesting phenomena has been explained by numerical approaches, e.g.
the $1:1$ resonance for Moon-Earth system which make the people can only see
one side of the Moon from the Earth. However, the Mercury-Sun model shows a
different $3:2$ resonance because of the large eccentricity, see Fig. 2.

Due to Theorem 1.4 and Theorem 1.6, such a resonance seems to be explained by
the following aspect: any trajectory within the global attractor $\Omega$ of
(24) has a longtime stability of velocity, namely, the average velocity is
close to certain rotation number, or even asymptotic to it. In Sec. LABEL:s5
we will show that variational minimal trajectories indeed match this
description.

###### Remark 2.1

As a further simplification, a spin-orbit model with $\eta(e)$ being a
constant is also widely concerned, which is actually a conformally symplectic
system. In CCD they further discussed the existence of KAM torus for this
model and proved the local attraction of the KAM torus.

###  2.3 Pumping of the swing

The pumping of a swing is usually modeled as a rigid object forced to rotate
back and forth at the lower ends of supporting ropes. After a series of
approximations and reasonable simplifications, the pumping of the swing can be
characterized as a harmonic oscillator with driving and parametric terms Ca2.
Therefore, this model has a typical meaning in understanding the dynamics of
motors.

As shown in Fig. 3, the length of the ropes supporting the swinger is $l$, and
$s$ is the distance between the center of mass of the swinger to the lower
ends of the rope. The angle of the supporting rope to the vertical position is
denoted by $\phi$, and the angle between the symmetric axis of the swinger and
the rope is $\theta$, which varies as $\theta=\theta_{0}\cos\omega t$. So we
get the equation of the motion by

$\displaystyle(l^{2}-2ls\cos\theta+s^{2}+R^{2})\ddot{\phi}$ $\displaystyle=$
$\displaystyle-gl\sin\phi+gs\sin(\phi+\theta)-ls\sin\theta\dot{\theta}^{2}$
$\displaystyle+(ls\cos\theta-s^{2}-R^{2})\ddot{\theta}-2ls\sin\theta\dot{\theta}\dot{\phi},\quad\phi\in\mathbb{T}$

where $g$ is the gravity index and $mR^{2}$ is the moment of inertia of the
center ($m$ is the mass of swinger). We can see that by reasonable adjustment
of $l,s,\omega$ parameters, this system can be dissipative, accelerative or
critical.

Notice that numerical research of this equation for $|\phi|\ll 1$ has been
done by numerical experts in a bunch of papers, see PGDB for a survey of that.
Those results successfully simulate the swinging at small to modest
amplitudes. As the amplitude grows these results become less and less
accurate, and that’s why we resort to a theoretical analysis in this paper.

Figure 3: A simulation of the pumping of the swing

##  3 Weak KAM solution of (HJ+)

Due to the superlinearity of $L(x,v,t)$, for each $k\geq 0$, there exists
$C(k)\geq 0$, such that

$L(x,v,t)\geq k|v|-C(k),k>0,x\in M.$

Moreover, the compactness of $M$ implies that for each $k>0$, there exists
$C_{k}>0$ such that

$\max_{\begin{subarray}{c}(x,t)\in M\times\mathbb{T}\\\ |v|\leq
k\end{subarray}}L(x,v,t)\leq C_{k}.$

###  3.1 Weak KAM solution of (HJ+) in the condition (H0-)

Note that $[f]>0$. The following conclusion can be easily checked.

###### Lemma 3.1

Assume $t>s$, then

  1. 1.

$F(s)-F(t)\leq 2k_{0}-(t-s-1)[f];$

  2. 2.

$\int^{t}_{s}e^{F(\tau)-F(t)}\mbox{d}\tau\leq\frac{e^{2k_{0}+[f]}}{[f]}\big{(}1-e^{-(t-s)[f]}\big{)};$

  3. 3.

$\int^{t}_{-\infty}e^{F(\tau)-F(t)}\mbox{d}\tau\leq\frac{e^{2k_{0}+[f]}}{[f]},$

where $k_{0}=\max_{s\in[0,2]}\big{|}\int^{s}_{0}f(\tau)\mbox{d}\tau\big{|}$.

Now we define a function $u_{\alpha}^{-}:M\times\mathbb{R}\to\mathbb{R}$ by

$\displaystyle u_{\alpha}^{-}(x,t)$ $\displaystyle:=$
$\displaystyle\inf\int^{t}_{-\infty}e^{F(s)-F(t)}(L(\gamma(s),\dot{\gamma}(s),s)+\alpha)\mbox{d}s$
(26)

where the infimum is taken for all $\gamma\in
C^{ac}((-\infty,t],M)$222aboslutely continuous curves with $\gamma(t)=x$. We
can easily prove this function is bounded, since

$-|C(k=0)-\alpha|\cdot\frac{e^{2k_{0}+[f]}}{[f]}\leq
u_{\alpha}^{-}(x,t)\leq|C_{k=0}+\alpha|\frac{e^{2k_{0}+[f]}}{[f]},$

where $C(0)$ and $C_{0}$ have been defined in the beginning of Sec. 3.

###### Lemma 3.2

[(1) of Theorem 1.2] $u_{\alpha}^{-}(x,t)$ is 1-periodic with respect to $t$,
i.e.,

$u_{\alpha}^{-}(x,t+1)=u_{\alpha}^{-}(x,t).$

###### Proof

By the definition of $u_{\alpha}^{-}$,

$\displaystyle u_{\alpha}^{-}(x,t+1)$ $\displaystyle=$
$\displaystyle\inf_{\gamma(t+1)=x}\bigg{\\{}\int^{t+1}_{-\infty}e^{F(s)-F(t+1)}\big{(}L(\gamma(s),\dot{\gamma}(s),s)+\alpha\big{)}\mbox{d}s\bigg{\\}}$
$\displaystyle=$
$\displaystyle\inf_{\gamma(t+1)=x}\bigg{\\{}\int^{t}_{-\infty}e^{F(s+1)-F(t+1)}\big{(}L(\gamma(s+1),\dot{\gamma}(s+1),s+1)+\alpha\big{)}\mbox{d}s\bigg{\\}}$
$\displaystyle=$
$\displaystyle\inf_{\eta(t)=x}\bigg{\\{}\int^{t}_{-\infty}e^{F(s)-F(t)}\big{(}L(\eta(s),\dot{\eta}(s),s)+\alpha\big{)}\mbox{d}s\bigg{\\}}$
$\displaystyle=$ $\displaystyle u_{\alpha}^{-}(x,t)$

as desired.∎

###### Lemma 3.3

[(3) of Theorem 1.2] Let $\gamma:[s_{1},s_{2}]\to M$ be an absolutely
continuous curve. Then,

$\displaystyle
e^{F(s_{2})}u_{\alpha}^{-}(\gamma(s_{2}),s_{2})-e^{F(s_{1})}u_{\alpha}^{-}(\gamma(s_{1}),s_{1})$
$\displaystyle\leq$
$\displaystyle\int^{s_{2}}_{s_{1}}e^{F(\tau)}(L(\gamma(\tau),\dot{\gamma}(\tau),\tau)+\alpha)\mbox{d}\tau.$

###### Proof

Let $\\{\gamma_{n}\\}$ be a sequence of absolutely continuous curve from
$(-\infty,s_{1}]$ to $M$ with $\gamma_{n}(s_{1})=\gamma(s_{1})$, such that

$e^{F(s_{1})}u_{\alpha}^{-}(\gamma(s_{1}),s_{1})=\lim_{n\to\infty}\int_{-\infty}^{s_{1}}e^{F(\tau)}(L(\gamma_{n}(\tau),\dot{\gamma}_{n}(\tau),\tau)+\alpha)\mbox{d}\tau.$

Let $\hat{\gamma}_{n}=\gamma_{n}*\gamma$ for each $n\in\mathbb{N}$. Hence,

$\displaystyle e^{F(s_{2})}u_{\alpha}^{-}(\gamma(s_{2}),s_{2})$
$\displaystyle\leq$
$\displaystyle\int^{s_{2}}_{-\infty}e^{F(\tau)}(L(\hat{\gamma}_{n}(\tau),\dot{\hat{\gamma}}_{n}(\tau),\tau)+\alpha)\mbox{d}\tau$
$\displaystyle\leq$
$\displaystyle\int^{s_{2}}_{s_{1}}e^{F(\tau)}(L(\gamma(\tau),\dot{\gamma}(\tau),\tau)+\alpha)\mbox{d}\tau$
$\displaystyle+\int^{s_{1}}_{-\infty}e^{F(\tau)}(L(\gamma_{n}(\tau),\dot{\gamma}_{n}(\tau),\tau)+\alpha)\mbox{d}\tau.$

Taking the limit $n\to\infty$, we derive (3.3) is true.∎

###### Lemma 3.4

For each $(x,t)\in M\times\mathbb{R}$ and $s<t$, it holds

$\displaystyle e^{F(t)}u_{\alpha}^{-}(x,t)$ $\displaystyle=$
$\displaystyle\inf_{\begin{subarray}{c}\gamma\in C^{ac}([s,t],M)\\\
\gamma(t)=x\end{subarray}}\bigg{\\{}e^{F(s)}u_{\alpha}^{-}(\gamma(s),s)+\int^{t}_{s}e^{F(\tau)}(L(\gamma(\tau),\dot{\gamma}(\tau),\tau)+\alpha)\mbox{d}\tau\bigg{\\}}.$

Moreover, the infimum in (3.4) can be achieved by a $C^{r}$ smooth minimizer.

###### Proof

Due to Lemma 3.3,

$e^{F(t)}u_{\alpha}^{-}(x,t)\leq\inf_{\begin{subarray}{c}\gamma\in
C^{ac}([s,t],M)\\\
\gamma(t)=x\end{subarray}}\big{\\{}e^{F(s)}u_{\alpha}^{-}(\gamma(s),s)+\int^{t}_{s}e^{F(\tau)}(L(\gamma,\dot{\gamma},\tau)+\alpha)\mbox{d}\tau\big{\\}}.$

For each $\epsilon>0$, there exists an absolutely continuous curve
$\gamma:(-\infty,t]\to M$ with $\gamma(t)=x$, such that

$\displaystyle e^{F(t)}u_{\alpha}^{-}(x,t)+\epsilon$ $\displaystyle\geq$
$\displaystyle\int^{t}_{-\infty}e^{F(\tau)}(L(\gamma(\tau),\dot{\gamma}(\tau),\tau)+\alpha)\mbox{d}\tau$
$\displaystyle=$
$\displaystyle\int^{t}_{s}e^{F(\tau)}(L(\gamma(\tau),\dot{\gamma}(\tau),\tau)+\alpha)\mbox{d}\tau$
$\displaystyle+\int^{s}_{-\infty}e^{F(\tau)}(L(\gamma(\tau),\dot{\gamma}(\tau),\tau)+\alpha)\mbox{d}\tau$
$\displaystyle\geq$ $\displaystyle
e^{F(s)}u_{\alpha}^{-}(\gamma(s),s)+\int^{t}_{s}e^{F(\tau)}(L(\gamma(\tau),\dot{\gamma}(\tau),\tau)+\alpha)\mbox{d}\tau.$

Hence, (3.4) proves to be an equality. Therefore, we can find a sequence of
absolutely continuous curve $\\{\gamma_{n}\\}$ with $\gamma_{n}(t)=x$ such
that

$e^{F(t)}u_{\alpha}^{-}(x,t)=\lim_{n\to\infty}\bigg{\\{}e^{F(s)}u_{\alpha}^{-}(\gamma_{n}(s),s)+\int^{t}_{s}e^{F(\tau)}(L(\gamma_{n},\dot{\gamma}_{n},\tau)+\alpha)\mbox{d}\tau\bigg{\\}}.$

Hence, there exists a constant $c$ independent of $n$, such that

$\int^{t}_{s}e^{F(\tau)}(L(\gamma_{n}(\tau),\dot{\gamma}_{n}(\tau),\tau)+\alpha)\mbox{d}\tau\leq
c.$ (29)

Due to Dunford-Petti Theorem (Theorem 6.4 in DFIZ), there exists a subsequence
$\\{\gamma_{n_{k}}\\}$ converging to a curve $\gamma_{*}$ such that

$\int^{t}_{s}e^{F(\tau)}(L(\gamma_{*},\dot{\gamma}_{*},\tau)+\alpha)\mbox{d}\tau\leq\varliminf_{k\to\infty}\int^{t}_{s}e^{F(\tau)}(L(\gamma_{n_{k}},\dot{\gamma}_{n_{k}},\tau)+\alpha)\mbox{d}\tau.$
(30)

Hence, (3.4) can be achieved at $\gamma_{*}:[s,t]\to M$, which definitely
solves the Euler-Lagrange equation (E-L). Due to the Weierstrass Theorem in
Mat, $\gamma_{*}$ is $C^{r}$ smooth.∎

###### Lemma 3.5

[(4) of Theorem 1.2] For each $\alpha\in\mathbb{R}$ and $(x,t)\in
M\times\mathbb{R}$, there exists a curve $\gamma_{x,t}^{-}:(-\infty,t]\to M$
with $\gamma_{x,t}^{-}(t)=x$ such that for each $t_{1}<t_{2}\leq t$,

$\displaystyle
e^{F(t_{2})}u_{\alpha}^{-}(\gamma_{x,t}^{-}(t_{2}),t_{2})-e^{F(t_{1})}u_{\alpha}^{-}(\gamma_{x,t}^{-}(t_{1}),t_{1})$
$\displaystyle=$
$\displaystyle\int^{t_{2}}_{t_{1}}e^{F(\tau)}(L(\gamma_{x,t}^{-}(\tau),\dot{\gamma}_{x,t}^{-}(\tau),\tau)+\alpha)\mbox{d}\tau.$

###### Proof

By Lemma 3.4, for each $n\in\mathbb{N}$, there exists a sequence of $C^{r}$
curve $\gamma_{n}:[t-n,t]\to M$ with $\gamma_{n}(t)=x$ such that

$e^{F(t)}u_{\alpha}^{-}(x,t)=e^{F(t-n)}u_{\alpha}^{-}(\gamma_{n}(t-n),t-n)+\int^{t}_{t-n}e^{F(\tau)}(L(\gamma_{n}(\tau),\dot{\gamma}_{n}(\tau),\tau)+\alpha)\mbox{d}\tau.$

It is easy to see for each interval $[a,b]\subset[t-n,t]$

$\displaystyle
e^{F(b)}u_{\alpha}^{-}(\gamma_{n}(b),b)-e^{F(a)}u_{\alpha}^{-}(\gamma_{n}(a),a)$
$\displaystyle=$
$\displaystyle\int^{b}_{a}e^{F(\tau)}(L(\gamma_{n}(\tau),\dot{\gamma}_{n}(\tau),\tau)+\alpha)\mbox{d}\tau.$

Due to diagonal approach, we derive there exists a subsequence of
$\\{\gamma_{n}\\}$, denoted by $\\{\gamma_{n_{k}}\\}$ and a curve
$\gamma_{x,t}^{-}:(-\infty,t]\to M$ such that $\gamma_{n_{k}}$ converges
uniformly to $\gamma_{x,t}^{-}$ on each finite subinterval of $(-\infty,t]$.
Taking $k\to\infty$ in (3.1), due to Dunford Petti Theorem,

$\displaystyle
e^{F(b)}u_{\alpha}^{-}(\gamma_{x,t}^{-}(b),b)-e^{F(a)}u_{\alpha}^{-}(\gamma_{x,t}^{-}(a),a)$
$\displaystyle=\varliminf_{k\to\infty}\int^{b}_{a}e^{F(\tau)}(L(\gamma_{n_{k}},\dot{\gamma}_{n_{k}},\tau)+\alpha)\mbox{d}\tau\geq\int^{b}_{a}e^{F(\tau)}(L(\gamma_{x,t}^{-},\dot{\gamma}_{x,t}^{-},\tau)+\alpha)\mbox{d}\tau.$

Combining with (3.3), we get (3.5). Since $\gamma_{x,t}^{-}|_{[s,t]}$ is a
minimizer of (3.4) for each $s<t$, due to Lemma 3.4, $\gamma_{x,t}^{-}$ is
$C^{r}$ and solves (E-L).∎

###### Remark 3.6

Due to Lemma 3.7, if we take $t_{2}=t$ and make $t_{1}\rightarrow-\infty$ in
(3.5) we instantly get

$u_{\alpha}^{-}(x,t)=\int^{t}_{-\infty}e^{F(\tau)-F(t)}(L(\gamma_{x,t}^{-}(\tau),\dot{\gamma}_{x,t}^{-}(\tau),\tau)+\alpha)\mbox{d}\tau,$

i.e. the infimum in (26) is achieved at
$\gamma_{x,t}^{-}:(-\infty,t]\rightarrow M$.

###### Lemma 3.7

Suppose $\gamma_{x,\theta}^{-}:(-\infty,\theta]\to M$ is a backward calibrated
curve ending with $x$ of $u_{\alpha}^{-}(x,\theta)$, then

$|\dot{\gamma}_{x,\theta}^{-}(\tau)|\leq\kappa_{0},\quad\forall\ (x,\theta)\in
M\times\mathbb{T},\tau<\theta.$

for a constant $\kappa_{0}$ depending on $L$ and $\alpha$. That implies
$\gamma_{x,\theta}^{-}$ is actually Lipschitz on $(-\infty,\theta]$.

###### Proof

Let $s_{1},s_{2}\leq\theta$ and $s_{2}-s_{1}=1$. Due to Lemma 3.5,

$\displaystyle
e^{F(s_{2})}u^{-}_{\alpha}(\gamma^{-}_{x,\theta}(s_{2}),s_{2})-e^{F(s_{1})}u^{-}_{\alpha}(\gamma^{-}_{x,\theta}(s_{1}),s_{1})$
$\displaystyle=\int^{s_{2}}_{s_{1}}e^{F(\tau)}(L(\gamma^{-}_{x,\theta}(\tau),\dot{\gamma}^{-}_{x,\theta}(\tau),\tau)+\alpha)\mbox{d}\tau$
$\displaystyle\geq\int^{s_{2}}_{s_{1}}e^{F(\tau)}(|\dot{\gamma}^{-}_{x,\theta}(\tau)|-C(1)+\alpha)\mbox{d}\tau.$

On the other hand, let $\beta:[s_{1},s_{2}]\to M$ be a geodesic satisfying
$\beta(s_{1})=\gamma^{-}_{x,\theta}(s_{1}),\beta(s_{2})=\gamma^{-}_{x,\theta}(s_{2})$,
and $|\dot{\beta}(\tau)|\leq\mbox{diam(M)}=:k_{1}$. Then

$\displaystyle
e^{F(s_{2})}u^{-}_{\alpha}(\gamma^{-}_{x,\theta}(s_{2}),s_{2})-e^{F(s_{1})}u^{-}_{\alpha}(\gamma^{-}_{x,\theta}(s_{1}),s_{1})$
$\displaystyle\leq\int^{s_{2}}_{s_{1}}e^{F(\tau)}(L(\beta(\tau),\dot{\beta}(\tau),\tau)+\alpha)\mbox{d}\tau$
$\displaystyle\leq\int^{s_{2}}_{s_{1}}e^{F(\tau)}(C_{k_{1}}+\alpha)\mbox{d}\tau.$

Hence,

$\int^{s_{2}}_{s_{1}}e^{F(\tau)}|\dot{\gamma}^{-}_{x,\theta}(\tau)|\mbox{d}\tau\leq\int^{s_{2}}_{s_{1}}e^{F(\tau)}(C_{k_{1}}+C(1))\mbox{d}\tau.$

Due to the continuity of $\dot{\gamma}^{-}_{x,\theta}(\tau)$, there exists
$s_{0}\in(s_{1},s_{2})$ such that

$|\dot{\gamma}^{-}_{x,\theta}(s_{0})|\leq C_{k_{1}}+C(1).$ (33)

Note that $\gamma^{-}_{x,\theta}$ solves (E-L), so
$|\dot{\gamma}^{-}_{x,\theta}(\tau)|$ is uniformly bounded for $(x,\theta)\in
M\times\mathbb{T}$ and $\tau\in(-\infty,\theta]$.∎

###### Lemma 3.8

[(2) of Theorem 1.2] For each $\alpha\in\mathbb{R}$, $u_{\alpha}^{-}$ is
Lipschitz on $M\times\mathbb{T}$.