# Series expansion of the excess work using nonlinear response theory

Pierre Nazé<EMAIL_ADDRESS>Instituto de Física ‘Gleb Wataghin’,
Universidade Estadual de Campinas, 13083-859, Campinas, São Paulo, Brazil
Marcus V. S. Bonança<EMAIL_ADDRESS>Instituto de Física ‘Gleb
Wataghin’, Universidade Estadual de Campinas, 13083-859, Campinas, São Paulo,
Brazil

###### Abstract

The calculation of observable averages in non-equilibrium regimes is one of
the most important problems in statistical physics. Using the Hamiltonian
approach of nonlinear response theory, we obtain a series expansion of the
average excess work and illustrate it with specific examples of thermally
isolated systems. We report the emergence of non-vanishing contributions for
large switching times when the system is subjected to strong driving. The
problem is solved by using an adapted multiple-scale method to supress these
secular terms. Our paradigmatic examples show how the method is implemented
generating a truncated series that obeys the Second Law of Thermodynamics.

##  I Introduction

Linear response theory lacks access to regions far from equilibrium for
treating only weak perturbations. Nonlinear response theory appears then as
one of the few alternatives to go beyond the linear regime when exact
solutions are absent Kubo (1957); Bernard and Callen (1959); Tani (1964);
Tanaka  _et al._ (1967); Peterson (1967); Samokhin (1968); Kenkre (1971);
Samokhin (1972); Kenkre (1973); Murayama (1981). In this approach, the non-
equilibrium average of an observable is expressed as a series expansion that
is expected to be a good approximation of the exact result. In this work, we
consider the non-equilibrium average of the generalized force that appears in
the expression of the thermodynamic work

$W=\int_{0}^{\tau}\overline{\frac{\partial\mathcal{H}}{\partial\lambda}}\frac{\partial\lambda}{\partial
t}dt$ (1)

when a control parameter $\lambda$ is varied in time. By subtracting the
corresponding quasistatic work from the previous expression, one defines the
so-called excess work Allahverdyan and Nieuwenhuizen (2005, 2007); Sivak and
Crooks (2012); Acconcia and Bonança (2015) whose value is expected to vanish
as the time $\tau$, required to perform the change in $\lambda$, increases.
Such constraint is imposed by the Second Law of Thermodynamics in either
isothermal or adiabatic processes Jarzynski (2020).

Equation (1) tells us that by expanding the non-equilibrium average of
$\partial\mathcal{H}/\partial\lambda$, we will automatically have a series
expansion of the thermodynamic work whose behavior has to agree with the
above-mentioned constraint. However, it is well-known that perturbative
approaches to time-dependent problems often lead to asymptotic series and
present secular terms Arnold  _et al._ (2007); Lichtenberg and Lieberman
(2013); Logan (2013). To exemplify this ubiquitous behavior, consider the
Duffing equation

$\ddot{u}+u=-\epsilon u^{3},$ (2)

which can be understood as a rescaled harmonic oscillator subjected to an
external nonlinear force term. When $\epsilon\ll 1$, the solution could be
expanded in a series on the parameter $\epsilon$

$u(t)=u_{0}(t)+\epsilon u_{1}(t)+\epsilon^{2}u_{2}(t)+....$ (3)

Considering the initial conditions $u(0)=1$ and $\dot{u}(0)=0$, the solution
expanded up to first-order becomes

$u(t)=\cos{t}+\epsilon\left[\frac{1}{32}(\cos{3t}-\cos{t})-\frac{3}{8}t\sin{t}\right],$
(4)

which clearly diverges for large $t$, even if $\epsilon$ is small.

It is well known that what produces such divergence is the fact that the
perturbation affects not only the amplitudes but also the frequencies of the
oscillatory functions appearing in the expansion Arnold  _et al._ (2007);
Lichtenberg and Lieberman (2013); Logan (2013). In some situations, the
secular terms might appear in higher orders and the truncation of the series
to its lowest orders may lead to a good approximation. However, there are
techniques in which the use of multiple-time scale strategies to eliminate the
secular terms leads to a well-behaved series expansion. Arnold  _et al._
(2007); Lichtenberg and Lieberman (2013); Logan (2013). Along these lines, one
of the first methods ever used in this kind of problem is the Lindstedt-
Poincaré method Lindstedt (1882); Poincaré (1893) in the context of classical
perturbation theory in celestial mechanics Lichtenberg and Lieberman (2013);
Fasano and Marmi (2006). The idea behind it is simple and can be rephrased as
follows: besides the series expansion of the amplitude, as given by Eq. (3),
the instant of time $t$ must be also rescaled according to

$\tau=t(\omega_{0}+\epsilon\omega_{1}+\epsilon^{2}\omega_{2}+...),$ (5)

where $\omega_{0}$ is the frequency of oscillation of the unperturbed system
and the $\omega_{i}$, $i=1,2,...$, are not just frequency corrections but also
help us to eliminate the secular terms. In this work, we will use an adapted
Lindstedt-Poincaré method to systematically fix the series expansion of the
excess work upto an arbitrary order. In other words, the method presented here
supresses the emergence of secular terms and hence produces a meaningful
series expansion that agrees with the Second Law of Thermodynamics.

This work is organized in the following way. In Sec. II, we present the
Hamiltonian approach of nonlinear response theory. In Sec. III, we exemplify
our method using the thermally isolated harmonic oscillator subjected to a
strong linear driving in the stiffness parameter. We verify the emergence of
secular terms in the series expansion obtained from nonlinear response theory
and show how to suppress them employing the ideas of Lindstedt-Poincaré
method, both in classical and quantum cases. In Sec. IV we make our final
considerations.

##  II Nonlinear response theory

###  II.1 Preliminaries

Consider a classical system with a Hamiltonian
$\mathcal{H}(\mbox{\boldmath$z$}(\mbox{\boldmath$z_{0}$},t),\lambda(t))$,
where $z$ is a point in the phase space $\Gamma$ evolved from the initial
point $\mbox{\boldmath$z$}_{0}$ according to Hamiltonian dynamics of
$\mathcal{H}$ and $\lambda(t)$ is a control parameter. During a switching time
$\tau$, the control parameter is changed from $\lambda_{0}$ to
$\lambda_{0}+\delta\lambda$. Initially, this system is in contact with a heat
bath of inverse temperature $\beta\equiv{(k_{B}T)}^{-1}$, where $k_{B}$ is
Boltzmann’s constant and its Hamiltonian is
$\mathcal{H}_{0}=\mathcal{H}(\mbox{\boldmath$z_{0}$},\lambda_{0})$. Given some
observable $\vartheta(\mbox{\boldmath$z_{0}$},\lambda(t))$, our interest is to
calculate its non-equilibrium average

$\overline{\vartheta}(t)=\int_{\Gamma}\vartheta(\mbox{\boldmath$z$}(\mbox{\boldmath$z_{0}$},t),\lambda(t))\rho(\mbox{\boldmath$z$}(\mbox{\boldmath$z_{0}$},t))d{\mbox{\boldmath$z_{0}$}},$
(6)

where $\rho(\mbox{\boldmath$z$}(\mbox{\boldmath$z_{0}$},t))$ is the non-
equilibrium phase-space distribution of the system. Such quantity evolves
according to the Liouville equation

$\dot{\rho}=-\\{\rho,\mathcal{H}\\}\equiv{\mathcal{L}}\rho,\quad\rho(\mbox{\boldmath$z_{0}$},0)=\rho_{c}(\mbox{\boldmath$z_{0}$})$
(7)

where $\\{\cdot,\cdot\\}$ is the Poisson bracket, ${\mathcal{L}}$ the
Liouville operator, and $\rho_{c}:=\exp{(-\beta\mathcal{H}_{0})}/\mathcal{Z}$
is the canonical ensemble, where
$\mathcal{Z}=\int_{\Gamma}\exp{(-\beta\mathcal{H}_{0})}d\mbox{\boldmath$z_{0}$}$.
In the particular case where the Hamiltonian $\mathcal{H}=\mathcal{H}_{0}$,
the formal solution of the Liouville equation is

$\rho(\mbox{\boldmath$z$}(\mbox{\boldmath$z_{0}$},t))=e^{t\mathcal{L}_{0}}\rho_{c}(\mbox{\boldmath$z_{0}$}),$
(8)

where ${\mathcal{L}}_{0}$ is the Liouville operator associated to
$\mathcal{H}_{0}$. We remark also that the time evolution operator
$e^{-t\mathcal{L}_{0}}$ can act on observables
$\vartheta(\mbox{\boldmath$z_{0}$})$, indicating how they evolve in time

$\vartheta(\mbox{\boldmath$z$}(\mbox{\boldmath$z_{0}$},t))=e^{-t\mathcal{L}_{0}}\vartheta(\mbox{\boldmath$z_{0}$}),$
(9)

where $\mbox{\boldmath$z$}(\mbox{\boldmath$z_{0}$},t)$ is the solution of
Hamilton equations considering the Hamiltonian $\mathcal{H}_{0}$.

The control parameter can be expressed as

$\lambda(t)=\lambda_{0}+g(t)\delta\lambda,$ (10)

where the protocol $g(t)$ must satisfy the following boundary conditions

$g(0)=0,\quad g(\tau)=1.$ (11)

We also consider that $g(t)\equiv g(t/\tau)$, which means that the time
intervals are measured in units of the switching time $\tau$. Additionally,
the quantity $\delta\lambda/\lambda_{0}$ quantifies the driving strength of
the process.

In linear response theory, the response of the observable is restricted to
first-order in $\delta\lambda/\lambda_{0}$. This can be achieved by expanding
the instantaneous observable and the non-equilibrium ensemble in a power
series of $\delta\lambda/\lambda_{0}$, regrouping the corresponding terms and
keeping only the first-order term. In nonlinear response theory, the procedure
is the same, but the truncation goes further than that used in linear response
theory. Our goal is to calculate the average work performed by the switch of
the control parameter $\lambda$ beyond the linear-response regime. Its
expression is given by Eq. (1) and requires a series expansion of the
observable $\partial\mathcal{H}/\partial\lambda$. In the following, we develop
a systematic way of collecting the different orders given by non-linear
response theory. This method will be exemplified in Sec. III.

###  II.2 Expansion in higher-order terms

Our first step is to rewrite Eq. (6) as a power series of the driving strength

$\overline{\vartheta}(t)=\sum_{n=0}^{\infty}{\overline{\vartheta}}_{n}(t)\left(\frac{\delta\lambda}{\lambda_{0}}\right)^{n}.$
(12)

To find the functions $\overline{\vartheta}_{n}(t)$, consider then the
following expansions

$\vartheta(\mbox{\boldmath$z$}(\mbox{\boldmath$z_{0}$},t),\lambda(t))=\sum_{n=0}^{\infty}\vartheta_{n}({\bf
z_{0}},t)\left(\frac{\delta\lambda}{\lambda_{0}}\right)^{n},$ (13)

and

$\rho(\mbox{\boldmath$z$}(\mbox{\boldmath$z_{0}$},t))=\sum_{n=0}^{\infty}\rho_{n}({\bf
z_{0}},t)\left(\frac{\delta\lambda}{\lambda_{0}}\right)^{n}.$ (14)

Knowing the observable
$\vartheta(\mbox{\boldmath$z$}(\mbox{\boldmath$z_{0}$},t),\lambda(t))$, the
functions $\vartheta_{n}({\bf z_{0}},t)$ are the coefficients of the Taylor
expansion of $\vartheta$ in the control parameter $\lambda$ around
$\lambda_{0}$

$\vartheta_{n}({\bf
z_{0}},t)=\lambda_{0}^{n}\frac{g^{n}(t)}{n!}\frac{\partial^{n}\vartheta}{\partial\lambda^{n}}\Bigr{|}_{\lambda=\lambda_{0}},$
(15)

which have been multiplied by $\lambda_{0}^{n}$ to adjust the powers of the
driving strength. We remark that when the observable depends on the control
parameter, each term of the expansion is an instantaneous response due to the
variation of such quantity.

The terms $\rho_{n}({\bf z_{0}},t)$ are much more involved to obtain. They
require a manipulation of Eq. (7), which is explained in the following. First,
we consider the expansion of the Hamiltonian as a power series on the driving
strength

$\mathcal{H}=\sum_{n=0}^{\infty}\mathcal{H}_{n}(\mbox{\boldmath$z$}(\mbox{\boldmath$z_{0}$},t))g^{n}(t)\left(\frac{\delta\lambda}{\lambda_{0}}\right)^{n},$
(16)

whose coefficients are given by

$\mathcal{H}_{n}(\mbox{\boldmath$z$}(\mbox{\boldmath$z_{0}$},t))=\frac{\lambda_{0}^{n}}{n!}\frac{\partial^{n}\mathcal{H}}{\partial\lambda^{n}}\Bigr{|}_{\lambda=\lambda_{0}}.$
(17)

We can construct a Liouville operator $\mathcal{L}_{n}$ for each
$\mathcal{H}_{n}$

$\displaystyle\mathcal{L}(\bullet)$ $\displaystyle=-\\{\bullet,\mathcal{H}\\}$
(18)
$\displaystyle=\sum_{n=0}^{\infty}g^{n}(t)\left(\frac{\delta\lambda}{\lambda_{0}}\right)^{n}\left[-\\{\bullet,\mathcal{H}_{n}\\}\right]$
$\displaystyle\equiv\sum_{n=0}^{\infty}g^{n}(t)\left(\frac{\delta\lambda}{\lambda_{0}}\right)^{n}\mathcal{L}_{n}(\bullet).$

In particular, we distinghish those Liouville operators which depend on time
from the time-independent term,

$\mathcal{L}=\mathcal{L}_{0}+\sum_{n=1}^{\infty}\mathcal{L}_{n}g^{n}(t)\left(\frac{\delta\lambda}{\lambda_{0}}\right)^{n}\equiv\mathcal{L}_{0}+\mathcal{L}_{\text{ext}}(t).$
(19)

For the case where the initial ensemble is the canonical one, the integral
form of the Liouville equation reads Kubo  _et al._ (1985)

$\rho(\mbox{\boldmath$z$}(\mbox{\boldmath$z_{0}$},t),t)=e^{t\mathcal{L}_{0}}\rho_{c}({\bf
z_{0}})\\\
+\int_{0}^{t}e^{(t-s)\mathcal{L}_{0}}\mathcal{L}_{\text{ext}}(s)\rho(\mbox{\boldmath$z$}(\mbox{\boldmath$z_{0}$},s),s)ds,$
(20)

which, using the Liouville theorem, becomes

$\rho(\mbox{\boldmath$z$}(\mbox{\boldmath$z_{0}$},t),t)=\rho_{c}({\bf
z_{0}})\\\
+\int_{0}^{t}e^{(t-s)\mathcal{L}_{0}}\mathcal{L}_{\text{ext}}(s)\rho(\mbox{\boldmath$z$}(\mbox{\boldmath$z_{0}$},s),s)ds.$
(21)

Expanding now $\rho$ and $\mathcal{L}$ in a power series on the driving
strength, we have

$\sum_{n=0}^{\infty}\rho_{n}({\bf
z_{0}},t)\left(\frac{\delta\lambda}{\lambda_{0}}\right)^{n}=\rho_{c}({\bf
z_{0}})\\\
+\sum_{i=1}^{\infty}\sum_{j=0}^{\infty}\int_{0}^{t}e^{(t-s)\mathcal{L}_{0}}\mathcal{L}_{i}\rho_{j}({\bf
z_{0}},s)g^{i}(s)ds\left(\frac{\delta\lambda}{\lambda_{0}}\right)^{i+j}.$ (22)

Finally, regroupping the terms with the same powers
$\left(\delta\lambda/\lambda_{0}\right)^{n}$, we finally have

$\displaystyle\rho_{0}({\bf z_{0}},t)$ $\displaystyle=\rho_{c}({\bf z_{0}}),$
(23) $\displaystyle\rho_{n}({\bf z_{0}},t)$
$\displaystyle=\sum_{k=0}^{n-1}\int_{0}^{t}e^{(t-s)\mathcal{L}_{0}}\mathcal{L}_{n-k}\rho_{k}({\bf
z_{0}},s)g^{n-k}(s)ds,$

for $n\geq 1$. Therefore, to find $\rho_{n}$ one needs to know the previous
solutions $\rho_{n-k}$, for $1\leq k\leq n$. We remark that the $\rho_{n}$
cannot be considered individually as valid probability distributions since
their integral over the phase space is zero Kubo  _et al._ (1985).

Finally, using Eqs. (15) and (23), the term ${\overline{\vartheta}}_{n}(t)$
can be rearranged as

${\overline{\vartheta}}_{n}(t)=\sum_{k=0}^{n}\int_{\Gamma}\vartheta_{k}({\bf
z_{0}},t)\rho_{n-k}({\bf z_{0}},t)d{\bf z_{0}}.$ (24)

Therefore, the averaged observable becomes

${\overline{\vartheta}}(t)=\sum_{n=0}^{\infty}\left[\sum_{k=0}^{n}\int_{\Gamma}\vartheta_{k}({\bf
z_{0}},t)\rho_{n-k}({\bf z_{0}},t)d{\bf
z_{0}}\right]\left(\frac{\delta\lambda}{\lambda_{0}}\right)^{n}.$ (25)

In the particular case where the observable $\vartheta$ is independent of the
control parameter $\lambda(t)$, its nonequilibirum average becomes

${\overline{\vartheta}}(t)=\sum_{n=0}^{\infty}\left[\int_{\Gamma}\vartheta({\bf
z_{0}})\rho_{n}({\bf z_{0}},t)d{\bf
z_{0}}\right]\left(\frac{\delta\lambda}{\lambda_{0}}\right)^{n}.$ (26)

###  II.3 Calculating non-equilibrium averages of arbitrary order

To calculate non-equilibrium averages of arbitrary order, one may proceed as
follows:

  1. 1.

Choose an order $n$ to truncate Eq. (12);

  2. 2.

For $0\leq k\leq n$, calculate $\vartheta_{k}$ using Eq. (15);

  3. 3.

For $0\leq k\leq n$, calculate $\mathcal{H}_{k}$ using Eq. (16);

  4. 4.

Calculate the solution $\mbox{\boldmath$z$}(\mbox{\boldmath$z_{0}$},t)$ of
$\mathcal{H}_{0}$;

  5. 5.

For $0\leq k\leq n$, calculate $\rho_{k}$ using Eq. (22);

  6. 6.

For $0\leq k\leq n$, calculate ${\overline{\vartheta}}_{k}$ using Eq. (24);

  7. 7.

Calculate ${\overline{\vartheta}}$ using Eq. (25).

In what follows, we will illustrate this procedure with specific examples.

####  II.3.1 Recovering linear response theory

Let us use the above-mentioned sequence of steps to recover the standard
result of linear response theory, which corresponds to the truncation at
$n=1$. We restrict ourselves to the case in which the observable of interest
does not depend on the control parameter $\lambda(t)$

$\vartheta({\bf z_{0}},\lambda(t))=\vartheta_{0}({\bf z_{0}}).$ (27)

For the sake of simplicity, we also consider that the expansion of the
Hamiltonian $\mathcal{H}$ ends up exactly at first order,

$\mathcal{H}({\bf z_{0}},\lambda(t))=\mathcal{H}_{0}({\bf
z_{0}})+\mathcal{H}_{1}({\bf
z_{0}})g(t)\left(\frac{\delta\lambda}{\lambda_{0}}\right).$ (28)

According to the procedure described previously, we only need the first two
terms of Eq. (14), whose coefficients read

$\rho_{0}({\bf z_{0}},t)=\rho_{c}({\bf z_{0}}),$ (29) $\rho_{1}({\bf
z_{0}},t)=\int_{0}^{t}e^{(t-s)\mathcal{L}_{0}}\mathcal{L}_{1}\rho_{c}({\bf
z_{0}})g(s)ds.$ (30)

The term ${\overline{\vartheta}}_{0}$ is given by

${\overline{\vartheta}}_{0}(t)=\int_{\Gamma}\vartheta_{0}({\bf
z_{0}})\rho_{c}({\bf z_{0}})d{\bf
z_{0}}\equiv\left\langle\vartheta_{0}\right\rangle,$ (31)

while the term ${\overline{\vartheta}}_{1}$ reads

${\overline{\vartheta}}_{1}(t)=\int_{\Gamma}\vartheta_{0}({\bf
z_{0}})\rho_{1}({\bf z_{0}},t)d{\bf z_{0}},$ (32)

which can be rewritten as

$\displaystyle\int_{\Gamma}$ $\displaystyle\vartheta_{0}({\bf
z_{0}})\rho_{1}({\bf z_{0}},t)d{\bf z_{0}}=$ (33)
$\displaystyle=\int_{\Gamma}\vartheta_{0}({\bf
z_{0}})\left[\int_{0}^{t}e^{(t-s)\mathcal{L}_{0}}\mathcal{L}_{1}\rho_{c}({\bf
z_{0}})g(s)ds\right]d{\bf z_{0}}$
$\displaystyle=-\int_{0}^{t}\left[\int_{\Gamma}\mathcal{L}_{1}e^{-(t-s)\mathcal{L}_{0}}\vartheta_{0}({\bf
z_{0}})\rho_{c}({\bf z_{0}})d{\bf z_{0}}\right]g(s)ds$
$\displaystyle=\int_{0}^{t}\phi_{1}(t-s)g(s)ds,$

where we have used the antisymmetric property of the Liouville operator Kubo
_et al._ (1985). The term $\phi_{1}(t)$ is the so-called (first-order)
response function, defined by

$\phi_{1}(t)\equiv\left\langle\\{e^{-t\mathcal{L}_{0}}\vartheta_{0}({\bf
z_{0}}),\mathcal{H}_{1}({\bf z_{0}})\\}\right\rangle,$ (34)

where the symbol $\left\langle(...)\right\rangle$ denotes the equilibrium
average on $\rho_{c}$. Therefore the non-equilibrium average of the observable
upto its first-order is

${\overline{\vartheta}}(t)=\left\langle\vartheta_{0}\right\rangle+\frac{\delta\lambda}{\lambda_{0}}\int_{0}^{t}\phi_{1}(t-s)g(s)ds,$
(35)

which is the standard result of linear response theory.

####  II.3.2 Calculating the second-order

To illustrate the calculation of the second-order term, we restrict ourselves
again to the situation in which the observable does not depend on the control
parameter and the Hamiltonian is exactly given by Eq. (28) . In this case, the
non-equilibrium average of $\vartheta$ upto second order is given by Eq. (35)
plus the following term multiplied by
$\left(\delta\lambda/\lambda_{0}\right)^{2}$

${\overline{\vartheta}}_{2}(t)=\int_{\Gamma}\vartheta({\bf
z_{0}})\rho_{2}({\bf z_{0}},t)d{\bf z_{0}},$ (36)

where, according to Eq. (23), we have

$\displaystyle\rho_{2}({\bf z_{0}},t)=\int_{0}^{t}\int_{0}^{s_{1}}$
$\displaystyle
e^{(t-s_{1})\mathcal{L}_{0}}\mathcal{L}_{1}e^{(s_{1}-s_{2})\mathcal{L}_{0}}\mathcal{L}_{1}\rho_{0}({\bf
z_{0}})$ (37) $\displaystyle g(s_{2})g(s_{1})ds_{2}ds_{1}.$

Using the antisymmetric property of the Liouville operators as we did in the
linear response case, we arrive at

${\overline{\vartheta}}_{2}(t)=\int_{0}^{t}\int_{0}^{s_{1}}\phi_{2}(t-s_{1},s_{1}-s_{2})g(s_{1})g(s_{2})ds_{1}ds_{2},$
(38)

where

$\phi_{2}(t,s)=\langle\\{e^{-s\mathcal{L}_{0}}\\{e^{-t\mathcal{L}_{0}}\vartheta_{0}({\bf
z_{0}}),\mathcal{H}_{1}({\bf z_{0}})\\},\mathcal{H}_{1}({\bf
z_{0}})\\}\rangle,$ (39)

is the second-order response function.

We remark that, even in higher orders, the nonequilibrium average can always
be expressed in terms of response functions that depend only on equilibrium
correlation functions. In other words, this is not an exclusive feature of
linear response theory, but a consequence of the perturbative method used in
Eq.(21) with a zero-order term given by an equilibrium ensemble.

##  III Series expansion of the excess work

In what follows, we compare the exact solution of a simple but relevant
example and the corresponding perturbative expression from non-linear response
theory. Both results are obtained for the excess thermodynamic work performed
along a finite-time process whose duration or switching time we denote by
$\tau$ Acconcia and Bonança (2015). This comparison will show that the
perturbative expansion described in the Sec.II leads to an asymptotic series.
However, the divergences of higher-order terms can be removed by the
application of a modified Lindstedt-Poincaré method Lindstedt (1882); Poincaré
(1893).

###  III.1 Thermally isolated harmonic oscillator

Figure 1: (Color online) Excess work for different switching times for the
exact time-dependent dynamics given by the Hamiltonian (40) and the protocol
$\lambda(t)$ of Eq. (41). The curve was obtained from the exact result for the
excess work $W_{\rm ex}$. The driving strength $\delta\lambda/\lambda_{0}$ was
chosen equal to 0.5.

We will apply the results of Sec. II to a thermally isolated harmonic
oscillator whose Hamiltonian is given by

$\mathcal{H}(\lambda(t))=\frac{1}{2}\left[p^{2}+\lambda(t)q^{2}\right],$ (40)

where the time-depedent stiffness $\lambda(t)$ is the control parameter, which
is driven according to the linear protocol

$\lambda(t)=\lambda_{0}+\delta\lambda\frac{t}{\tau}.$ (41)

The oscillator is initially in equilibrium with a heat bath of inverse
temperature $\beta=(k_{B}T)^{-1}$, which is removed before the starting of the
non-equilibrium driving. Hence, the dynamics of the oscillator is Hamiltonian
during the whole process. The average work performed on the system along the
process was defined in Eq. (1), where the term
$\overline{\partial\mathcal{H}}/{\partial\lambda}$ is interpreted as a
generalized force. For thermally isolated systems, the energetic cost due to
finite-time driving can be measured by the excess work,

$W_{\text{ex}}\equiv W-W_{\text{qs}},$ (42)

defined as the difference between the thermodynamic work, given by Eq. (1),
and the quantity $W_{\rm qs}$. The term $W_{\text{qs}}$ is the quasistatic
work, whose value is the difference between the final and initial average
energies of the system when it is driven along a quasistatic process Jarzynski
(2020). This quantity is obtained from the adiabatic invariant Fasano and
Marmi (2006) which, for the one-dimensional harmonic oscillator, is nothing
but the action or, equivalently, the area in phase space enclosed by the curve
of constant energy. In Appendix A, we calculate the quasistatic work for our
system, which is given by

$W_{\text{qs}}=\frac{1}{\beta}\left[\left(1+\frac{\delta\lambda}{\lambda_{0}}\right)^{1/2}-1\right].$
(43)

Figure 1 depicts the excess work performed along the protocol (41) for
different switching times. Each point of the curve shown in Fig. 1 corresponds
to the average of the difference between final and initial values of the
Hamiltonian (40). These averages are taken over the initial canonical
ensemble. For the linear protocol (41), Hamilton equations are analytically
solvable (see Appendix B). We observe that the excess work vanishes for large
$\tau$, in agreement with the Second Law of Thermodynamics.

In what follows, we obtain a perturbative expression for $W_{\rm ex}$ based on
non-linear response theory and verify that the decay for large $\tau$ is not
observed for higher-order terms. We denote by $W_{\text{ex}}^{(n)}$ the $n$th
contribution in the series expansion of the excess work.

###  III.2 Perturbative expression via non-linear response theory

To obtain a perturbative expression for $W_{\text{ex}}$ using non-linear
response theory, we apply the sequence of steps described in Sec. II.3.
Following the definitions established in Sec. II.2 (see Eq. (16)), the
Hamiltonian (40) under the protocol (41) can be split into two terms

$\mathcal{H}_{0}=\frac{1}{2}\left[p^{2}+\lambda_{0}q^{2}\right],\quad\mathcal{H}_{1}=\lambda_{0}\frac{q^{2}}{2}.$
(44)

The expansion of the excess work goes through the expansion of the generalized
force, which is an average of the following observable

$\frac{\partial\mathcal{H}}{\partial\lambda}=\frac{q^{2}}{2}.$ (45)

The calculation of the $(n+1)$th term of the excess work is based on finding
the response functions $\phi_{n}$

$\phi_{n}(t_{1},...,t_{n})=\left\langle\\{e^{-t_{1}\mathcal{L}_{0}}\\{...\\{e^{-t_{n-1}\mathcal{L}_{0}}\times\right.\\\
\left.\\{e^{-t_{n}\mathcal{L}_{0}}\vartheta_{0}({\bf
z_{0}}),\mathcal{H}_{1}({\bf z_{0}})\\},\mathcal{H}_{1}({\bf
z_{0}})\\}...,\mathcal{H}_{1}({\bf z_{0}})\\}\right\rangle,$ (46)

which are going to be used to calculate the $n$th term of the generalized
force, as exemplified in Eq. (38).

To calculate those terms, one can implement the sequence of steps described in
Sec. II.3 as a computer algorithm to obtain response functions of arbitrary
order. In Appendix C we deduce via mathematical induction the expansion of the
quasistatic work, which is necessary as well to find the $(n+1)$th term of the
excess work.

It is then straightforward to verify the emergence of secular terms for large
$\tau$, something that is in stark contrast with the exact result (see Fig.
1). For instance, the fourth-order term is already non-vanishing for large
$\tau$,

$W_{\text{ex}}^{(4)}=...+\frac{1}{\beta}\frac{\cos{(2\lambda_{0}^{1/2}\tau)}}{128}{\left(\frac{\delta\lambda}{\lambda_{0}}\right)}^{4}+...,$
(47)

and the fifth-order term clearly diverges in the same limit

$W_{\text{ex}}^{(5)}=...-\frac{1}{\beta}\frac{\lambda_{0}^{1/2}\tau\sin{(2\lambda_{0}^{1/2}}\tau)}{768}{\left(\frac{\delta\lambda}{\lambda_{0}}\right)}^{5}+....$
(48)

We depict in Fig. 2 the perturbative expansion of the excess work calculated
upto the 7th order. We observe that the result starts to diverge after $n=4$
for large $\tau$. However, for relatively small values of $\tau$, the
agreement with the exact result improves as higher-order terms are included in
the sum (see Fig. 3).

The method we will present in the next section furnishes a new series in which
the secular terms are suppressed. This reformulation of the expansion provided
by nonlinear response theory will then be shown to successfully achieve a
meaningful physical behavior for large values of $\tau$.

Figure 2: (Color online) Large $\tau$ behavior of the perturbative expansion
of the excess work upto the $n$th order, with $4\leq n\leq 7$. The results are
compared with the exact result (see Fig. 1). The driving strength was taken as
$\delta\lambda_{0}/\lambda_{0}=0.5$.
Figure 3: (Color online) Small $\tau$ behavior of the perturbative expansion
of the excess work upto the $n$th order, with $4\leq n\leq 7$. The results are
compared with the exact result (see Fig. 1). The driving strength was taken as
$\delta\lambda_{0}/\lambda_{0}=0.5$.

###  III.3 Supressing secular terms

The appearance of secular terms in a series expansion can be solved by using
the Lindstedt-Poincaré technique. The main idea is to consider that the
switching time can be written as a series in the driving strength

$\tau=\tau^{\prime}\sum_{n=0}^{\infty}\frac{a_{n}}{n!}{\left(\frac{\delta\lambda}{\lambda_{0}}\right)}^{n},$
(49)

where $\tau^{\prime}$ is the rescaled switching time. The coefficients $a_{n}$
are obtained by demanding the removal of all divergent or non-vanishing terms
(for large $\tau$) at each order This becomes possible once Eq. (49) is
inserted in the perturbative expansion obtained from nonlinear response theory
and all functions of $\tau$ are expanded in powers of
$\delta\lambda/\lambda_{0}$. Then, new coefficients for each order have to be
found and the coefficients $a_{n}$ can be determined. Once the coefficients
$a_{n}$ are determined (upto a certain order), the whole power series is
rearranged. Below, we describe with some detail how the coefficients
$a_{0},a_{1},a_{2}$, and $a_{3}$ are obtained.

The series expansion of $W_{\rm ex}$ starts with a second-order term that is
well-behaved,

$W_{\text{ex}}^{(2)}=\frac{\sin^{2}\left(a_{0}\lambda_{0}^{1/2}\tau^{\prime}\right)}{8a_{0}^{2}\beta\lambda_{0}\tau^{\prime
2}}\left(\frac{\delta\lambda}{\lambda_{0}}\right)^{2}.$ (50)

As mentioned before, the first secular terms show up in $W^{(4)}_{\rm ex}$. To
remove them, we plug Eq. (49) truncated upto second order in the expressions
of $W^{(2)}_{\rm ex}$, $W^{(3)}_{\rm ex}$ and $W^{(4)}_{\rm ex}$. We then
expand all the terms upto the fourth-order and obtain the following
combination of secular terms

$W_{\text{ex}}^{(4)}=...-\frac{a_{1}^{2}\sin^{2}\left(a_{0}\lambda_{0}^{1/2}\tau^{\prime}\right)}{8a_{0}^{2}\beta}+\frac{a_{1}^{2}\cos^{2}\left(a_{0}\lambda_{0}^{1/2}\tau^{\prime}\right)}{8a_{0}^{2}\beta}\\\
+\frac{a_{1}\cos\left(2a_{0}\lambda_{0}^{1/2}\tau^{\prime}\right)}{16a_{0}\beta}+\frac{\cos\left(2a_{0}\lambda_{0}^{1/2}\tau^{\prime}\right)}{128\beta}+...,$
(51)

First, we observe that $a_{0}=1$, since $\tau$ is not rescaled at zeroth
order. Thus these terms can be suppressed if we choose $a_{1}=-1/4$. This is
our first example of how to obtain the coefficients of Eq. (49).

The expansion mentioned above upto the fourth-order also involves the
coefficient $a_{2}$. However, its value is determined only when we evaluate
the series expansion upto the sixth order and demand that the secular terms
appearing there vanish. Analogously, we have to go to the 9th order to
determine the value of $a_{3}$. This procedure yields to

$a_{0}=1,\quad a_{1}=-\frac{1}{4},\quad a_{2}=\frac{5}{24},\quad
a_{3}=-\frac{5}{16}.$ (52)

Substituting those values in the expression of the excess work expanded with
Eq. (49), we produce a well-behaved series of the excess work upto 5th order

$W_{\text{ex}}^{(5)}=\frac{\sin^{2}(\lambda_{0}^{1/2}\tau^{\prime})}{8\beta\lambda_{0}\tau^{\prime
2}}{\left(\frac{\delta\lambda}{\lambda_{0}}\right)}^{2}-\frac{\sin^{2}(\lambda_{0}^{1/2}\tau^{\prime})}{16\beta\lambda_{0}\tau^{\prime
2}}{\left(\frac{\delta\lambda}{\lambda_{0}}\right)}^{3}\\\
+\left(\frac{77}{768\beta\lambda_{0}\tau^{\prime
2}}-\frac{13}{128\beta\lambda_{0}^{2}\tau^{\prime
4}}+\frac{17\sin(2\lambda_{0}^{1/2}\tau^{\prime})}{256\beta\lambda_{0}^{3/2}\tau^{\prime
3}}\right.\\\
\left.-\frac{23\cos(2\lambda_{0}^{1/2}\tau^{\prime})}{768\beta\lambda_{0}\tau^{\prime
2}}+\frac{13\cos(2\lambda_{0}^{1/2}\tau^{\prime})}{128\beta\lambda_{0}^{2}\tau^{\prime
4}}\right){\left(\frac{\delta\lambda}{\lambda_{0}}\right)}^{4}\\\
+\left(-\frac{69}{512\beta\lambda_{0}\tau^{\prime
2}}+\frac{39}{256\beta\lambda_{0}^{2}\tau^{\prime
4}}-\frac{51\sin(2\lambda_{0}^{1/2}\tau^{\prime})}{512\beta\lambda_{0}^{3/2}\tau^{\prime
3}}\right.\\\
\left.+\frac{15\cos(2\lambda_{0}^{1/2}\tau^{\prime})}{512\beta\lambda_{0}\tau^{\prime
2}}-\frac{39\cos(2\lambda_{0}^{1/2}\tau^{\prime})}{256\beta\lambda_{0}^{2}\tau^{\prime
4}}\right){\left(\frac{\delta\lambda}{\lambda_{0}}\right)}^{5}.$ (53)

Figure 4 presents a comparison between the corrected perturbative expression
upto 5th order and the exact result. It was used a driving strength
$\delta\lambda/\lambda_{0}=0.5$, which is about the maximum value before a
deviation appears between both results. The value of $\tau^{\prime}$ is
calculated from Eq. (49) using the coefficients $a_{i}$ obtained and the value
of $\tau$ given as a boundary condition.

We remark that rescaling the switching time, whose value is one of the
boundary conditions of the driving, is just a practical way of rescaling the
frequency $\lambda_{0}^{1/2}$ since in all orders we have functions of the
product $\lambda_{0}^{1/2}\tau$ Lichtenberg and Lieberman (2013); Logan
(2013).

Figure 4: (Color online) Comparison between the exact result (black dots) and
the corrected series expansion (blue solid line) of the excess work upto fifth
order. It was used a driving strength of $\delta\lambda/\lambda_{0}=0.5$.

###  III.4 Extension to quantum case

The method developed previously can easily be extended to the quantum case.
Consider a quantum system with a Hamiltonian $\hat{\mathcal{H}}$. The first
modification occurs at the Liouville equation, which becomes the Liouville-von
Neumann equation

$\dot{\hat{\rho}}=\hat{\mathcal{L}}\hat{\rho}:=\frac{i}{\hbar}[\hat{\rho},\hat{H}],$
(54)

where $\hat{\mathcal{L}}$ is the quantum Liouville operator, $\hat{\rho}$ is
the density matrix and $[\cdot,\cdot]$ is the commutator. In the Heisenberg
picture, the quantum Liouville operator evolves a operator $\hat{\vartheta}$
in time in the following way

$e^{-t\hat{\mathcal{L}}}\hat{\vartheta}:=e^{i\frac{\hat{\mathcal{H}_{0}}t}{\hbar}}\hat{\vartheta}e^{-i\frac{\hat{\mathcal{H}_{0}}t}{\hbar}}.$
(55)

Also the equilibrium average of an observable $\hat{\vartheta}$ is defined as

$\langle\hat{\vartheta}\rangle:={\rm Tr}(\hat{\vartheta}\hat{\rho_{c}}),$ (56)

where $\hat{\rho_{c}}$ is the quantum canonical ensemble.

The development of a power series in the quantum case is based on a procedure
identical to the one developed in Sec. II, but using the quantum Liouville
operator. For example, for a quantum system of Hamiltonian

$\hat{\mathcal{H}}=\hat{\mathcal{H}}_{0}+\hat{\mathcal{H}}_{1}g(t)\left(\frac{\delta\lambda}{\lambda_{0}}\right),$
(57)

the quantum analogues to the first and second response functions are

$\phi_{1}(t)=\frac{i}{\hbar}{\rm
Tr}([e^{-t\hat{\mathcal{L}}}\hat{\vartheta}_{0},\hat{\mathcal{H}_{1}}]),$ (58)
$\phi_{2}(t)=-\frac{1}{\hbar^{2}}{\rm
Tr}([e^{-t\hat{\mathcal{L}}}[e^{-t\hat{\mathcal{L}}}\hat{\vartheta}_{0},\hat{\mathcal{H}_{1}}],\hat{\mathcal{H}_{1}}]).$
(59)

We applied this extension to the quantum case to calculate the thermodynamic
work given by Eq. 1 of the quantum analogue of Eq. (40). The time-dependent
Hamiltonian of such system is

$\hat{\mathcal{H}}(t)=\frac{\hat{p}^{2}}{2}+\left(\lambda_{0}+\frac{t}{\tau}\delta\lambda\right)\frac{\hat{q}^{2}}{2},$
(60)

where $\hat{q}$ and $\hat{p}$ are the position and momentum operators
respectively. The initial equilibrium situation was taken at temperature
$T=0$, which means that all equilibrium averages were taken in the ground
state of the unperturbed Hamiltonian.

In the present case, quantum nonlinear response theory leads to an asymptotic
series of the excess work that is identical to the classical case with the
only modification being

$\frac{1}{\beta}\rightarrow\frac{\hbar\sqrt{\lambda_{0}}}{2}:=E_{0}.$ (61)

For instance, the leading order reads,

$W_{\rm
ex}^{(2)}=\frac{\hbar\sin^{2}(\lambda_{0}^{1/2}\tau^{\prime})}{16\lambda_{0}^{1/2}\tau^{\prime
2}}{\left(\frac{\delta\lambda}{\lambda_{0}}\right)}^{2}.$ (62)

Thus the asymptotic series are corrected by the same coeffcients given in Eq.
(52). In Fig. 5 we compare the analytical result of the excess work with the
corrected series given by the nonlinear response theory and the Lindstedt-
Poincaré method. The fact that the quantum and classical series are the same
(upto a different pre-factor) can be understood as a consequence of the
existing analytical solution in both cases. In the case of a linear protocol
for the stiffness parameter, it has been shown that Heisenberg’s equations are
identical to the classical case Husimi (1953); Deffner and Lutz (2008).

Figure 5: (Color online) Comparison between the exact result (black dots) and
the corrected series expansion (blue solid line) of the excess work upto fifth
order for a time-dependent quantum harmonic oscillator. It was used a driving
strength of $\delta\lambda/\lambda_{0}=0.5$.

##  IV Final remarks

This work reports the existence of secular terms in the perturbative expansion
of the excess work obtained from nonlinear response theory. To recover a
meaningful result for large switching times, we employed a multiple-scale
method that successfully suppresses these secular terms. This approach implies
a rescaling of the frequency of oscillation observed in the behavior of the
excess work. The implementation of such rescaling was made through an
expansion of the switching time in powers of the driving strength. In
classical canonical perturbation theroy, this is done only after the
application of canonical transformations to action-angle variables which, in
the context of nonequilibrium statistical physics, can be an unwanted step.
Thus, our method has the advantage of keeping the formalism of nonlinear
response theory unchanged. This rescaling of the frequency affects
applications such as shortcuts to adiabaticity based on response theory
Acconcia  _et al._ (2015) and possible future applications on optimization of
the nonequilibrium work. In summary, we have shown how to obtain a meaningful
series expansion of the excess work introducing a simple additional tool to
the standard nonlinear response theory.

###### Acknowledgements.

P. N. gratefully thank Otavio L. Canton for the presentation which suggested
the key to solve this problem. This work was financially supported for FAPESP
(Fundação de Amparo à Pesquisa do Estado de São Paulo) (Brazil) (Grant No.
2018/06365-4 and No. 2018/21285-7) and for CNPq (Conselho Nacional de
Desenvolvimento Científico e Pesquisa) (Brazil) (Grant No. 141018/2017-8).

##  Appendix A Calculation of $W_{\rm qs}$

The quasistatic work is given by the difference of the equilibrium averages of
the energies calculated along the quasistatic process

$W_{\rm qs}=\left\langle E_{\rm f}^{\rm ad}\right\rangle-\left\langle E_{\rm
i}\right\rangle,$ (63)

which are obtained by the conservation of the adiabatic invariant
$\Omega(E,\lambda)$. This quantity is given by the area enclosed by the energy
shell, which can be written as

$\Omega(E,\lambda)=\int_{0}^{E}d\mathcal{H}\int_{0}^{2\pi}\mathcal{J}(\theta,\mathcal{H})d\theta,$
(64)

where $\mathcal{J}$ is the Jacobian of the action-angle transformation
$(q,p)\rightarrow(\theta,\mathcal{H})$. In the particular case of the harmonic
oscillator, the adiabatic invariant is

$\Omega(E,\lambda)=\frac{2\pi E}{\sqrt{\lambda}}.$ (65)

If the adiabatic invariant is conserved along the quasistatic process, it
holds

$E_{\rm f}^{\rm ad}=E_{\rm i}{\left(\frac{\lambda_{\rm f}}{\lambda_{\rm
i}}\right)}^{1/2}.$ (66)

Using Eq. (66) on Eq. (63), the quasistatic work becomes

$W_{\text{qs}}=\frac{1}{\beta}\left[\left(\frac{\lambda_{\rm f}}{\lambda_{\rm
i}}\right)^{1/2}-1\right].$ (67)

##  Appendix B Solutions of Eq. (40)

Given the initial conditions $q(0)=q_{0}$ and $p(0)=p_{0}$, the time-dependent
Hamiltonian of Eq. (40) has the following exact solution

$q(t)=\frac{\pi}{\sqrt[3]{\delta\lambda}}\left[\text{Ai}\left(-\frac{t\delta\lambda+\lambda_{0}\tau}{\delta\lambda^{2/3}\sqrt[3]{\tau}}\right)\left(\sqrt[3]{\tau}p_{0}\text{Bi}\left(-\frac{\lambda_{0}\tau^{2/3}}{\delta\lambda^{2/3}}\right)\right.\right.\\\
\left.+\sqrt[3]{\delta\lambda}q_{0}\text{Bi}^{\prime}\left(-\frac{\lambda_{0}\tau^{2/3}}{\delta\lambda^{2/3}}\right)\right)-\text{Bi}\left(-\frac{t\delta\lambda+\lambda_{0}\tau}{\delta\lambda^{2/3}\sqrt[3]{\tau}}\right)\times\\\
\left.\left(\sqrt[3]{\tau}p_{0}\text{Ai}\left(-\frac{\lambda_{0}\tau^{2/3}}{\delta\lambda^{2/3}}\right)+\sqrt[3]{\delta\lambda}q_{0}\text{Ai}^{\prime}\left(-\frac{\lambda_{0}\tau^{2/3}}{\delta\lambda^{2/3}}\right)\right)\right],$
(68)
$p(t)=\frac{\pi}{\sqrt[3]{\tau}}\left[\text{Bi}^{\prime}\left(-\frac{t\delta\lambda+\lambda_{0}\tau}{\delta\lambda^{2/3}\sqrt[3]{\tau}}\right)\left(\sqrt[3]{\tau}p_{0}\text{Ai}\left(-\frac{\lambda_{0}\tau^{2/3}}{\delta\lambda^{2/3}}\right)\right.\right.\\\
\left.+\sqrt[3]{\delta\lambda}q_{0}\text{Ai}^{\prime}\left(-\frac{\lambda_{0}\tau^{2/3}}{\delta\lambda^{2/3}}\right)\right)-\pi\text{Ai}^{\prime}\left(-\frac{t\delta\lambda+\lambda_{0}\tau}{\delta\lambda^{2/3}\sqrt[3]{\tau}}\right)\times\\\
\left.\left(\sqrt[3]{\tau}p_{0}\text{Bi}\left(-\frac{\lambda_{0}\tau^{2/3}}{\delta\lambda^{2/3}}\right)+\sqrt[3]{\delta\lambda}q_{0}\text{Bi}^{\prime}\left(-\frac{\lambda_{0}\tau^{2/3}}{\delta\lambda^{2/3}}\right)\right)\right],$
(69)

where $\text{Ai}(x)$ and $\text{Bi}(x)$ are respectively the Airy functions of
first and second type, which are defined as

$\text{Ai}(x)\equiv\frac{1}{\pi}\int_{0}^{\infty}\cos{\left(\frac{t^{3}}{3}+xt\right)}dt,$
(70)
$\text{Bi}(x)\equiv\frac{1}{\pi}\int_{0}^{\infty}\left[\exp{\left(-\frac{t^{3}}{3}+xt\right)}+\sin{\left(\frac{t^{3}}{3}+xt\right)}\right]dt.$
(71)

##  Appendix C Expansion of $W_{\text{qs}}$

We are going to show that the quasistatic work

$W_{\text{qs}}=\frac{1}{\beta}\left[\left(1+\frac{\delta\lambda}{\lambda_{0}}\right)^{1/2}-1\right]$
(72)

can be expressed as

$W_{\text{qs}}=\frac{1}{\beta}\sum_{n=1}^{\infty}{(-1)}^{n+1}\frac{(2n-3)!!}{2^{n}n!}\left(\frac{\delta\lambda}{\lambda_{0}}\right)^{n}.$
(73)

Consider, by analogy with the function of Eq. (72), the following function

$f(x)=(1+x)^{1/2}-1.$ (74)

Being $f(0)=0$, we state that the $n$-th derivative, for $n\geq 1$, is given
by

$f^{(n)}(x)=(-1)^{n+1}\frac{(2n-3)!!}{2^{n}}{(1+x)}^{-\frac{2n-1}{2}}.$ (75)

It is easy to see that the result holds for $n=1$. Taking the derivative of
Eq. (75), we have

$f^{(n+1)}(x)=(-1)^{(n+1)+1}\frac{(2(n+1)-3)!!}{2^{(n+1)}}{(1+x)}^{-\frac{2(n+1)-1}{2}},$
(76)

which proves, by mathematical induction, our statement. As the Taylor
expansion is given by

$f(x)=f(0)+\sum_{n=1}^{\infty}\frac{f^{(n)}(0)}{n!}x^{n},$ (77)

Eq. (73) holds by applying Eq. (75) in Eq. (77).

## References

  * Kubo (1957) Ryogo Kubo, “Statistical-mechanical theory of irreversible processes. i. general theory and simple applications to magnetic and conduction problems,” Journal of the Physical Society of Japan 12, 570–586 (1957).
  * Bernard and Callen (1959) William Bernard and Herbert B Callen, “Irreversible thermodynamics of nonlinear processes and noise in driven systems,” Reviews of Modern Physics 31, 1017 (1959).
  * Tani (1964) Kensuke Tani, “A formula of non-linear responses,” Progress of Theoretical Physics 32, 167–169 (1964).
  * Tanaka  _et al._ (1967) Tomoyasu Tanaka, Kishin Moorjani, and Tohru Morita, “Green’s-function theory of nonlinear transport coefficients,” Phys. Rev. 155, 388–392 (1967).
  * Peterson (1967) Robert L Peterson, “Formal theory of nonlinear response,” Reviews of Modern Physics 39, 69 (1967).
  * Samokhin (1968) AA Samokhin, “Theory of nonlinear response of an isolated spin system,” Physica 39, 541–559 (1968).
  * Kenkre (1971) V. M. Kenkre, “Integrodifferential equation for response theory,” Physical Review A 4, 2327 (1971).
  * Samokhin (1972) AA Samokhin, “Theory of nonlinear response of an isolated spin system. ii,” Physica 58, 26–36 (1972).
  * Kenkre (1973) V. M. Kenkre, “Equations for the theory of response and transport in statistical mechanics,” Phys. Rev. A 7, 772–781 (1973).
  * Murayama (1981) Yoshimasa Murayama, “Theory of nonlinear nonequilibrium response: Application to the esaki effect,” Physica A: Statistical Mechanics and its Applications 109, 251–264 (1981).
  * Allahverdyan and Nieuwenhuizen (2005) A. E. Allahverdyan and Th. M. Nieuwenhuizen, “Minimal work principle: Proof and counterexamples,” Phys. Rev. E 71, 046107 (2005).
  * Allahverdyan and Nieuwenhuizen (2007) A. E. Allahverdyan and Th. M. Nieuwenhuizen, “Minimal-work principle and its limits for classical systems,” Phys. Rev. E 75, 051124 (2007).
  * Sivak and Crooks (2012) David A. Sivak and Gavin E. Crooks, “Thermodynamic metrics and optimal paths,” Phys. Rev. Lett. 108, 190602 (2012).
  * Acconcia and Bonança (2015) T. V. Acconcia and M. V. S. Bonança, “Degenerate optimal paths in thermally isolated systems,” Physical Review E 91, 042141 (2015).
  * Jarzynski (2020) Christopher Jarzynski, “Fluctuation relations and strong inequalities for thermally isolated systems,” Physica A: Statistical Mechanics and its Applications 552, 122077 (2020).
  * Arnold  _et al._ (2007) Vladimir I Arnold, Valery V Kozlov, and Anatoly I Neishtadt,  _Mathematical aspects of classical and celestial mechanics_, Vol. 3 (Springer Science & Business Media, 2007).
  * Lichtenberg and Lieberman (2013) Allan J Lichtenberg and Michael A Lieberman,  _Regular and chaotic dynamics_, Vol. 38 (Springer Science & Business Media, 2013).
  * Logan (2013) J David Logan,  _Applied mathematics_ (John Wiley & Sons, 2013).
  * Lindstedt (1882) Anders Lindstedt, “Abh. k,” Akad. Wiss. St. Petersburg 31 (1882).
  * Poincaré (1893) Henri Poincaré,  _Les méthodes nouvelles de la mécanique céleste: Méthodes de MM. Newcomb, Glydén, Lindstedt et Bohlin. 1893_, Vol. 2 (Gauthier-Villars it fils, 1893).
  * Fasano and Marmi (2006) Antonio Fasano and Stefano Marmi,  _Analytical mechanics: an introduction_ (OUP Oxford, 2006).
  * Kubo  _et al._ (1985) R. Kubo, M. Toda, and N. Hashitsume,  _Statistical physics II: nonequilibrium statistical mechanics_, Vol. 31 (Springer-Verlag, Berlin, 1985).
  * Husimi (1953) Kôdi Husimi, “Miscellanea in elementary quantum mechanics, ii,” Progress of Theoretical Physics 9, 381–402 (1953).
  * Deffner and Lutz (2008) Sebastian Deffner and Eric Lutz, “Nonequilibrium work distribution of a quantum harmonic oscillator,” Phys. Rev. E 77, 021128 (2008).
  * Acconcia  _et al._ (2015) Thiago V. Acconcia, Marcus V. S. Bonança, and Sebastian Deffner, “Shortcuts to adiabaticity from linear response theory,” Phys. Rev. E 92, 042148 (2015).