Quantum Work Fluctuations in connection with Jarzynski Equality
Juan D. Jaramillo Jiawen Deng Jiangbin Gong




Department of Physics, National University of Singapore, Singapore 117546

NUS Graduate School for Integrative Science and Engineering, Singapore 117597
 phygj@nus.edu.sg
Department of Physics, National University of Singapore, Singapore 117546
NUS Graduate School for Integrative Science and Engineering, Singapore 117597




A result of great theoretical and experimental interest, Jarzynski equality predicts a free energy change Δ F of a system at inverse temperature β from an ensemble average of non-equilibrium exponential work, i.e., ⟨ e^-β W⟩ =e^-βΔ F.
The number of experimental work values needed to reach a given accuracy of Δ F is determined by the variance
of e^-β W, denoted var(e^-β W).

 We discover in this work that var(e^-β W) in both harmonic and an-harmonic Hamiltonian systems
 can systematically diverge in non-adiabatic work protocols, even when the adiabatic protocols do not suffer from such divergence. This divergence may be regarded as a type of dynamically induced phase transition in work fluctuations. For a quantum harmonic oscillator with time-dependent trapping frequency as a working example, any non-adiabatic work protocol is found to yield a diverging var(e^-β W) at sufficiently low temperatures, markedly different from the classical behavior. 
 The divergence of
var(e^-β W) indicates the too-far-from-equilibrium nature of a non-adiabatic work protocol and makes it compulsory
 to apply designed control fields to suppress the quantum work fluctuations in order to test Jarzynski equality.



Quantum Work Fluctuations in connection with Jarzynski Equality
    Jiangbin Gong
    December 30, 2023
===============================================================




























Introduction. Work fluctuation theorems are a central topic in modern non-equilibrium statistical mechanics <cit.>. One outstanding result is Jarzynski equality (JE) <cit.> which  makes use of an average of non-equilibrium exponential work, i.e., e^-β W (where W is the work and β is the inverse temperature) to obtain important equilibrium information.  It takes the form ⟨ e^-β W⟩ =e^-βΔ F, where ⟨·⟩ represents an ensemble average of all possible work values with the system initially prepared in a Gibbs state, Δ F is  the free energy difference between initial and final systems at the same β.  Regarded as a recent breakthrough in non-equilibrium statistical mechanics, JE holds no matter how the work protocol is executed, fast or slow, adiabatic or highly non-adiabatic.  The quantum version of JE  takes precisely the same form <cit.>, provided that the value of quantum work is interpreted as the energy difference obtained from two energy projective measurements.
 It is also convenient to define the so-called dissipated work  W_ dis= W-Δ F.  Then  JE gives ⟨ e^-β W_ dis⟩=1.


Early applications of JE focused on biomolecular systems, where it is a well-established method to compute Δ F <cit.>.  Proof-of-principle experiments were carried out <cit.> following an early proposal by Hummer and Szabo <cit.>.  Other experiments include those using classical oscillators <cit.>. On the quantum side, the first experiment testing JE was first reported for a quantum spin <cit.>, relying on interferometric schemes to reconstruct work distributions <cit.>.  While other settings have been proposed, e.g., circuit QED <cit.>, a more recent experimental confirmation, following the proposal <cit.>, uses a trapped ion subject to an external force <cit.>.  These experiments are by no means straightforward.  In both classical and quantum cases, the ensemble average of a highly nonlinear (in fact, exponential) function of work, as requested by JE, could become practically demanding in yielding a well-converged result for Δ F, as rare events can potentially dominate the average <cit.>.  This in fact has motivated a number of studies in the literature to investigate the errors in predicting Δ F based on JE  <cit.>.
As suggested by the central limit theorem (CLT), the errors are related to the variance in the exponential work, i.e., var(e^-β W)≡⟨ e^-2β W⟩- e^-2βΔ F.   The larger var(e^-β W) is, the more realizations of W we must collect to reach a given accuracy in Δ F.   To obtain Δ F within an error of k_BT, the number of W realizations needed was estimated as var(e^-β W_ dis)=  var(e^-β W)/e^-2βΔ F <cit.>.

 

 

 The quantity  var(e^-β W_ dis) represents the ensemble average of exponential quantities and the Boltzmann exponential cutoff for high-energy states might not always be effective in suppressing the contributions from the high-energy tail.
 
 Indeed, in contrast to JE itself which is always well behaved, var(e^-β W_ dis) can diverge, suggesting  higher-energy components from the initial Gibbs distribution may contribute even more to the variance.  Surprisingly, prior to our work <cit.>, this possibility of divergence was rarely mentioned, with one exception <cit.> treating an adiabatic work protocol.  According to the principle of minimal work fluctuations (PMWF) <cit.>, once var(e^-β W_ dis) diverges for an adiabatic protocol, then it will suffer from analogous divergence under arbitrary work protocols sharing the same initial and final system Hamiltonians. The accuracy in estimating Δ F  in such situations becomes problematic.   One can no longer rely on the CLT to predict how the error scales with the number of experimental runs.   A generalization of the CLT can prove useful <cit.> for a specific situation, but in general the error estimate under a diverging var(e^-β W_ dis) is unknown.

  
 This work discovers the divergence of var(e^-β W_ dis) that is only present for non-adiabatic dynamics.
  Such a possibility is consistent with PMWF. This represents intriguing and more complex situations, where  the “phase" boundary between a finite var(e^-β W_ dis) and a diverging var(e^-β W_ dis) is yet to be located case by case.
 
 Our working example below is mainly a quantum harmonic oscillator (QHO) with time-dependent trapping frequency ω because recent experiments to test quantum JE <cit.> and to construct quantum heat engines <cit.> were based on QHO (we consider an-harmonic Hamiltonians in the end).  Work is done to the system by varying ω.  For a QHO of fundamental and experimental importance, all quantum transition probabilities can be analytically obtained.  
To our great surprise, even in such a prototypic system with all dynamical aspects known for so many years,  divergence of var(e^-β W_ dis) may be induced by non-adiabatic work protocols. The domain of divergence as a function of temperature and other dynamical parameters is found to possess complicated structures.  In general, so long as the dynamics is non-adiabatic, quantum effects at lower temperatures tend to enlarge the domain of divergence in var(e^-β W_ dis). Indeed, even when the classical var(e^-β W_ dis) is well behaved, its quantum counterpart is bound to diverge at sufficiently low temperatures.  This itself constitutes a new aspect of quantum effects in non-equilibrium statistical mechanics. Our results indicate (i) that direct applicability of JE in free energy estimates can become much limited as a system of interest approaches the deep quantum regime, and (ii) that quantum non-adiabatic effects or “inner friction" <cit.> alone can induce critical changes in work fluctuations. This work shall also motivate parallel studies on several generalized quantum fluctuation theorems <cit.>.

Work characteristic function and model system. Throughout we focus on work applied to an isolated system.
 Let H(t) be the time-dependent Hamiltonian subject to a work protocol.
The work distribution function is assumed to be P(W).  The Fourier transformation of P(W) yields the so-called work characteristic function <cit.>

    G(μ)=∫ dW e^iμ WP(W).

For a work protocol starting at t=0 and ending at t=τ,  G(μ) can be expressed as a quantum correlation function G(μ)= Tr[e^iuH_H(τ)e^-iuH(0)ρ(0)], where ρ(0) is the initial Gibbs state and H_H(τ) is the final Hamiltonian in the Heisenberg representation.  While the value of G(μ) at μ=iβ yields JE, the value of G(μ) at μ=2iβ determines  var(e^-β W_ dis).  That is,

    var(e^-β W_ dis)=G(2iβ)/e^-2βΔ F-1.

This observation indicates that a divergence in G(2iβ) is equivalent to a divergence in var(e^-β W_ dis) or in var(e^-β W).  The quantum correlation function G(μ) on the imaginary axis is hence our central object. 
An operative expression of (<ref>) is given by










    G(μ)=∑_m,n exp[iμ(E_n^τ-E_m^0)]P_m,n^0→τP_m^0,


where P^0→τ_m,n represents the transition probability from mth state of H(0) to nth state of H(τ). Equations (<ref>) and (<ref>) indicate that, for a system with a finite-dimension Hilbert space such as a spin system,
var(e^-β W_ dis) is obtained from a finite summation and will be always well behaved.


To lay a solid ground for our surprising findings, we aim to present our main results with computational and analytical calculations supporting each other.  With this in mind, the parametric QHO system seems to be the best choice as a model with an infinite-dimensional Hilbert space.  QHO belongs to the algebraic class SU(1,1) of integrable systems, for which many solutions in the form of exact propagators can be found in the literature <cit.>.  Systems in the same algebraic class share common dynamical features.  For example, the Calogero-Sutherland model despite being an interacting system is dynamically analogous to QHO <cit.>.  In addition, in the history of quantum mechanics and statistical physics, QHO plays a pivotal role, providing general insights into quantum zero-point energy, low-temperature behavior of the heat capacity of a crystal, open quantum systems, and the relationship between classical mechanics and quantum mechanics, etc.  In particular, in understanding quantum-classical correspondence, tuning the dimensionless ratio of the thermal energy over the energy level spacing is essential to make transitions between classical and quantum regimes.


Consider then a QHO with driving in its trapping frequency, with H_QHO(ω(t))=p̂^2/2m+mω^2(t)q̂^2/2.
Work is done as the trapping frequency changes from ω_0=ω(0) to ω_1=ω(τ).
 






Thanks to Husimi, all P^0→τ_m,n  are available as a function of the Husimi coefficient  Q^∗≥ 1 <cit.>, which is defined as
Q^∗(t)=1/2(ω_0ω_1X^2+ω_0/ω_1Ẋ^2+ω_1/ω_0Y^2+1/ω_0ω_1Y^2), where X(t) and Y(t) are the two solutions of the equation of motion of the corresponding classical  oscillator <cit.>.
For adiabatic dynamics, Q^∗=1 and P^0→τ_m,n=δ_m,n.
Any non-adiabaticity is captured by deviations of Q^∗ from its adiabatic value Q^∗=1.   For later use, we define the compression ratio κ≡ω_0/ω_1.  For a sudden quench (sq) protocol τ→ 0, then Q^∗→ Q^∗_ sq=(κ+1/κ)/2.

High-temperature and low-temperature regimes.
It is of interest to first examine the behavior of G(μ=2iβ) in the strictly adiabatic case, i.e., P^0→τ_m,n=δ_m,n.  Because the convergence criteria of the resulting
geometric series in Eq. (<ref>) becomes iμω_1-(iμ+β)ω_0<0 with μ=2iβ, it is obvious that if κ≥ 2, then G(2iβ) and hence var(e^-β W_ dis) diverge.  With PMWF, one further concludes that the divergence occurs in all non-adiabatic work protocols.

Next we investigate the behavior of G(2iβ) under non-adiabatic work protocols, provided that adiabatic work protocols do not
yield divergence in G(2iβ), i.e., in the regime of κ< 2.
To that end, we partially resort to a compact expression for G(μ) <cit.>.  That is,


    G(μ)=√(2) (1-e^-βħω_0)e^iμħ(ω_1-ω_0)/2/√(Q^∗(1-e^2iμħω_1)(1-e^-2(iμ+β)ħω_0)+(1+e^2iμħω_1)(1+e^-2(iμ+β)ħω_0)-4e^iμħω_1 e^-(iμ+β)ħω_0).


 Extra care is needed when one extends such a non-analytic
 result based on integration assuming real μ to the imaginary axis of μ. The Jarzynski equality is always recovered from (<ref>) at μ=iβ <cit.>, i.e.,
G(iβ)=sinh(βħω_0/2)/sinh(βħω_1/2).
Encouraged by this, we cautiously use this compact expression to help us digest the possible critical boundary of divergence in G(2iβ) in the regime of κ<2.  Equation (<ref>) then takes us to the following compact expression for G(2iβ),


    G(2iβ)→√(2)sinh(βħω_0/2)/√(cosh(βħ(2ω_1-ω_0))-1-(Q^∗-1) sinh(2βħω_1)sinh(βħω_0)).


The real denominator of G(2iβ) in Eq. (<ref>) approaching zero signifies a boundary separating finite values of G(2iβ) from its diverging values.  We stress however, all such “shortcut"  conclusions are carefully checked against the original expression in Eq. (<ref>), with P^0→τ_m,n in the sum series truncated at some large values of m and n (see Supplementary Material <cit.>).



 For the high-temperature (classical) regime βħω_0,1≪ 1 we have

    G(2iβ) →κ/√((2-κ)^2-4κ(Q^∗-1)).

The boundary of divergence in G(2iβ) is identified at  
    κ_c=2(Q^∗-√(Q^∗ 2-1)).

Equation  (<ref>) for adiabatic cases (Q^∗=1) still reproduces the known κ_c=2 boundary.
In generic non-adiabatic cases with Q^∗>1, the divergence boundary is pushed to smaller values, i.e., κ_c<2, in agreement with the classical result (see Supplementary Material <cit.>).    Later on we consider specific work protocols to further digest this result.  One important feature is that
this phase boundary at κ_c is no longer dependent upon temperature. Indeed, Eq. (<ref>) shows that G(2iβ) itself is temperature-independent in the high-temperature (classical) regime.   Expanding the denominator of Eq. (<ref>) around κ_c under a given Q^∗, we find var(e^-β W_ dis) diverges as ∼ (κ_c-κ)^-1/2 as κ approaches the critical value κ_c from below.






 

The behavior of G(2iβ) is entirely different in the low-temperature (deep quantum) regime βħω_0,1≫ 1, where
 
 the compact expression for G(2iβ) becomes

    G(2iβ) →exp[βħ(ω_0-ω_1)]/√(1-(Q^∗-1/2) [exp(2βħω_0)-1]).

For strictly adiabatic cases, i.e., Q^∗=1, Eq. (<ref>) is well behaved.  However, for any work protocol with even slight non-adiabaticity, Q^∗ 1, one observes that the denominator in Eq. (<ref>) hits zero or becomes imaginary at sufficiently low temperatures. This indicates that var(e^-β W_ dis) must diverge as temperature decreases, so long as the protocol is not strictly adiabatic. The boundary of divergence in var(e^-β W_ dis), within the low-temperature regime, is located at β_c= ln(Q^∗+1/Q^∗-1)/(2ħω_0), which does not explicitly depend on κ.  Further, from Eq. (<ref>) we find that as β approaches β_c, var(e^-β W_ dis) diverges as

∼ (β_c-β)^-1/2.
To understand the divergence,  one can actually focus on the contribution to var(e^-β W_ dis) made by the transitions from 2n-th state of H(0) to the ground state of H(τ). After applying Sterling's
formula to P^0→τ_2n,0 for large n, this contribution is found to scale as
1/√(n)(Q^⋆-1/Q^⋆+1)^ne^2nβħω_0.
This predicts the same β_c as above, because for β>β_c, the more rare the initial state is (sampled from the initial Boltzmann distribution), the more contribution it makes to var(e^-β W_ dis).  In other words,
quantum transitions associated with very negative work values, though exponentially suppressed by (Q^⋆-1/Q^⋆+1)^n with n, still not suppressed as sharply as
in classical cases so as to lose the competition from the exponential increasing factor e^2nβħω_0 contained in e^-2β W at sufficiently low temperatures. This important insight has been confirmed by the behavior of classical and quantum deformed Jarzynzki equalities we just proposed <cit.>.






All these predicted features are checked in numerics, with some examples shown in Fig. <ref>. There the results are obtained based on Eq. (<ref>) directly as well as from Eq. (<ref>). In the high-temperature regime (small β),  all the plots var(e^-β W_ dis) vs ħβ (in dimensionless units) become flat, in agreement with Eq. (<ref>). This also indicates that the classical var(e^-β W_ dis) does not diverge in these cases.  For lower temperatures, however,  all the plotted curves tend to blow up, except for the strict adiabatic case (bottom curve) whose var(e^-β W_ dis) approaches zero as temperature decreases, in agreement with Eq. (<ref>).
 As seen from Fig. <ref>, local minima in var(e^-β W_ dis) as a function of β might emerge, reflecting a competition between quantum fluctuations and thermal fluctuations. In addition,  the divergent behavior of var(e^-β W_ dis) close to κ_c or β_c (obtained from direct numerics and not shown) is also found to agree with the scaling laws of (κ_c-κ)^-1/2 and (β_c-β)^-1/2.






Specific work protocols at intermediate temperatures. To further investigate the behavior of var(e^ -β W_ dis), we turn to a specific work protocol where the trapping frequency is suddenly quenched from ω_0 to ω_1. In this case Q^∗=(κ+ 1/κ)/2. In the high-temperature regime, Eq. (<ref>) yields a critical κ_c=√(2), namely, var(e^ -β W_ dis) diverges if  ω_0≥√(2)ω_1.  Fig. <ref> depicts the numerically obtained domain of divergence in var(e^ -β W_ dis) in terms of ω_0 and ω_1, with panel (a) exactly showing this high-temperature behavior.  Quantum effects however dramatically enlarge the domain of divergence (white area) in var(e^ -β W_ dis):  From panel (b) to (d), temperature decreases and  the domain of divergence in var(e^ -β W_ dis) gradually invades the (classical) domain of finite var(e^ -β W_ dis) in panel (a).   For even lower temperatures than shown in Fig. <ref>, the domain of convergence (gray area) collapse into the line ω_0=ω_1, suggesting any actual quench in ω will yield divergence in var(e^ -β W_ dis).






 We have also examined a finite-time work protocol, where the parameter dω/dt/ω^2(t) is chosen to be time-independent, with ω(t)=ω_1ω_0τ/[ω_0τ+t(ω_1-ω_0)] <cit.>.





 
 Compared with the sudden quench case, the divergence domain is now fragmented into multiple domains with subtle phase boundaries.
 
In addition, as the temperature decreases, the domain of divergence again grows.
 


Similar disconnected regions of divergence are observed along the time of driving, where convergent domains are found around the adiabatic times <cit.> for which Q^∗ is close to 1. Detailed results are presented in Supplementary Material <cit.>.
 
 
Finally, it is worth mentioning the possibility of work protocols with Q^∗>Q^∗_sq, see for example <cit.>. In these cases var(e^-β W_ dis) may even diverge under compression of the harmonic potential (κ<1) at both high and low temperatures.












Before ending this section, we also mention our numerical studies of work fluctuations in quantum anharmonic oscillators and other systems such as particle in a box. There we gain similar insights. In particular, we consider a time-dependent quartic potential
with the total Hamiltonian H_a(t)=H_QHO(ω(t))+J  cos(tπ/τ) q̂^4, where τ is the time of driving, as well as a time-independent quartic potential, with the total Hamiltonian H_b(t)=H_QHO(ω(t))+J q̂^4.
There quantum effects are also found to induce the divergence of var(e^-β W_ dis) even if it does not diverge in the classical (high-temperature) domain.
Fig. <ref> depicts a few such computational examples.  In all these cases, the effect of a quartic potential is to bring the onset of divergence at higher temperatures as compared with the harmonic cases. 




Conclusions. JE was hoped to take advantage of non-equilibrium work protocols to estimate the change of free energy in nanoscale and mesoscopic systems. The divergence in var(e^-β W_ dis)  presents a hurdle to a direct application of JE.  Note that if the divergence is solely induced by non-adiabatic quantum effects, then one promising solution is to use shortcuts to adiabaticity (STA) <cit.>, to go around this divergence and yet still realizing speedy work protocols (but with some price <cit.>).  Indeed, one can even use var(e^-β W_ dis) as a minimization target to design control fields <cit.>.
Because quantum effects are seen to enlarge the domain of divergence drastically in the low temperature regime, using JE for free energy estimates in the deep quantum regime will not be fruitful in the absence of a designed control field.  Nevertheless,  the divergence in var(e^-β W_ dis) offers a new angle to study quantum work fluctuations <cit.> and to characterize the too-far-from-equilibrium nature of a work protocol.  One fascinating question is to study how the divergence of G(μ=2iβ) may be reflected in the  behavior of G(μ) in the real-μ domain.
Finally, one might wonder if the recent experiment testing quantum JE using a trapped ion <cit.> already suffered from the divergence issue exposed here.  It turns out that this belongs to a fortunate case without any divergence in var(e^-β W_ dis) (see Supplementary Material <cit.>). Thus, analyzing possible divergence in exponential work fluctuations will be crucial in guiding the design of future quantum experiments testing JE.


Acknowledgements:  This work
is funded by Singapore MOE Academic Research
Fund Tier-2 project (Project No. MOE2014-T2-2- 119 with
WBS No. R-144-000-350-112).  J.G. also acknowledges encouraging discussions with many colleagues.















 
 Appendix  
 

As explained in the main text, all divergences in var(e^-β W) are traced back to divergences in

    ⟨ e^-2β W⟩=∑_m,n exp[-2β(E_n^τ-E_m^0)]P_m,n^0→τP_m^0,

where P_m,n^0→τ are the state-to-state transition probabilities and P_m^0=e^-β E^0_m/Z_0 is the Gibbs distribution of the initial state.  All our numeric calculations involve a truncation in the series (<ref>).





§ TRANSITION PROBABILITIES





The transition probabilities between instant eigenstates at t=0 and t=τ for the quantum harmonic oscillator (QHO) with arbitrary driving of the trapping frequency are given by <cit.>


    P^0→τ_2ν,2μ   =   (2/Q^∗+1)^1/2(Q^∗-1/Q^∗+1)^ν+μ(2μ)!(2ν)!/2^2μ+2ν {∑_λ=0^min(ν,μ)[-8/(Q^∗-1)]^λ/(2λ)!(μ-λ)!(ν-λ)!}^2,
    
    P^0→τ_2ν+1,2μ+1   =   (2/Q^∗+1)^3/2(Q^∗-1/Q^∗+1)^ν+μ(2μ+1)!(2ν+1)!/2^2μ+2ν {∑_λ=0^min(ν,μ)[-8/(Q^∗-1)]^λ/(2λ+1)!(μ-λ)!(ν-λ)!}^2;



 where, ν,μ∈ℕ, and Q^∗≥ 1 is the so called Husimi coefficient.  The selection rules prevent transitions between energy levels with different parity.  
 The definition of Q^∗ is <cit.>

    Q^∗(t)=1/2(ω_0ω_1X^2+ω_0/ω_1Ẋ^2+ω_1/ω_0Y^2+1/ω_0ω_1Y^2),

 where X(t) and Y(t) are the two solutions of the equation of motion of the corresponding classical  oscillator,

Ẍ(t)+ω^2(t)X(t)=0,


with initial conditions: X(0)=Ẏ(0)=0 and Ẋ(0)=Y(0)=1.  Any non-adiabaticity is captured by deviations of Q^∗ from its adiabatic value, Q^∗=1.   The Q^∗ used in the main text is assumed to be Q^∗(τ).







§ WORK FLUCTUATIONS FOR THE CLASSICAL HARMONIC OSCILLATOR




For classical systems the work distribution is given by

    P^c(W)=βω_0/2π∫ dp_0dq_0 e^-βℋ_0(p_0,q_0) δ[W-(ℋ_τ(p_0,q_0)-ℋ_0(p_0,q_0))],

where ℋ_0(p_0,q_0) is the initial classical Hamiltonian and ℋ_τ(p_0,q_0) represents the value of the final Hamiltonian ℋ_τ for the trajectory emanating from (p_0,q_0).
 For the classical harmonic oscillator (CHO) this integral results in two different expressions, associated to the cases: Q^∗≤ Q^∗_sq and Q^∗>Q^∗_sq, where Q^∗ is the Husimi coefficient defined in (<ref>) and Q^∗_sq stands for the Husimi coefficient in sudden quench cases. 


 










 We first consider the regime Q^∗≤ Q^∗_sq, which will be the case if, for example, we have a monotonically increasing or decreasing ω(t) from ω_0 to ω_1. This gives rise to   positive definite work or negative definite work.   For compression of the harmonic trap (ω_0<ω_1) the work distribution reads <cit.>

    P^c_<(W)=β √(κ/2(Q^∗_sq-Q^∗)) exp[β κ-Q^∗/2(Q^∗_sq-Q^∗) W] I_0(β √(Q^∗ 2-1)/2(Q_sq^∗-Q^∗) |W|)Θ(W),

where κ≡ω_0/ω_1 is the compression ratio, Θ(W) is the step function and I_0(z) is the modified Bessel function of the first kind.  This work distribution coincides with that of the QHO in the semiclassical regime <cit.>.
The asymptotic behavior P^c(W→∞) determines the convergence of

    ⟨ e^-2β W⟩_c=∫ dW e^-2β WP^c(W).

 One can then use the approximation: I_0(z)→ e^z/√(2π z), z→∞, to get the tail

    e^-2β WP^c_<(W)→√(βκ/π√(Q^∗ 2-1)) exp[β (κ-Q^∗/2(Q^∗_sq-Q^∗)-2+√(Q^∗ 2-1)/2(Q_sq^∗-Q^∗)) W].

Evaluation of ⟨ e^-2β W⟩ requires us to integrate e^-2β WP^c(W) along W from 0 to ∞.  A divergence of ⟨ e^-2β W⟩_c will take place when the exponent in (<ref>) is an increasing function of W.  Fortunately, the divergence does not show up for compression of the trap. However, for expansion of the trap the work W is always negative, the work function reading

    P^c_<(W)=β √(κ/2(Q^∗_sq-Q^∗)) exp[β κ-Q^∗/2(Q^∗_sq-Q^∗) W] I_0(β √(Q^∗ 2-1)/2(Q_sq^∗-Q^∗) |W|)Θ(-W).

Similarly, the tail is

    e^-2β WP^c_<(W)→√(βκ/π√(Q^∗ 2-1)) exp[β (κ-Q^∗/2(Q^∗_sq-Q^∗)-2-√(Q^∗ 2-1)/2(Q_sq^∗-Q^∗)) W].

Since the term ⟨ e^-2β W⟩ is an integration along W from -∞ to 0, the possibility of a negative coefficient, i.e.,

    (κ-Q^∗/2(Q^∗_sq-Q^∗)-2-√(Q^∗ 2-1)/2(Q_sq^∗-Q^∗))<0

allows divergence. In particular, for adiabatic expansion of the trap (Q^∗=1, κ>1) the criteria of divergence is ω_0>2ω_1.


Now consider the regime where the nonadiabaticity parameter Q^∗ is even larger than that associated with a sudden quench, i.e., Q^∗>Q^∗_sq. This is possible if we introduce multiple quenches to the trapping frequency ω.  In this case, the integral (<ref>) for compression or expansion becomes

    P^c(W)=β √(κ/2π^2(Q^∗-Q^∗_sq)) exp[-β κ-Q^∗/2(Q^∗-Q^∗_sq) W] K_0(β √(Q^∗ 2-1)/2(Q^∗-Q_sq^∗) |W|),

where K_0(x) is the Macdonald function or modified Bessel function of the third kind. 
 This result is consistent with that reported in <cit.> using a different approach starting from the characteristic function of the quantum harmonic oscillator (QHO) and taking a semiclassical limit. Similarly, with the approximation K_0(z)→√(π)e^-z/√(2 z), z→∞, one obtains the tail

    e^-2β WP^c(W)→√(βκ/π√(Q^∗ 2-1)) exp{-β [(κ-Q^∗/2(Q^∗-Q^∗_sq)+2) W+√(Q^∗ 2-1)/2(Q^∗-Q^∗_sq) |W|]}.

 Because W now can take both positive and negative values,  the tail must go to zero at both W→±∞ in order to have a finite integral over W.  The condition for divergence after integration over W is


    κ-Q^∗/2(Q^∗-Q^∗_sq)+2+√(Q^∗2-1)/2(Q^∗-Q_sq^∗)<0,    for   W>0;
    
    -κ-Q^∗/2(Q^∗-Q^∗_sq)-2+√(Q^∗2-1)/2(Q^∗-Q_sq^∗)<0,    for   W< 0.



Briefly,

    |κ-Q^∗/2(Q^∗-Q^∗_sq)+2|>√(Q^∗2-1)/2(Q^∗-Q_sq^∗).

Combining these results and after a few straightforward steps,  we finally reach a compact condition for divergence

    κ>2(Q^∗-√(Q^∗2-1)).



Classical characteristic function.  First we consider the case Q^∗<Q^∗_sq.  From the integral formula

    ∫_0^+∞ dx e^Ax I_0(Bx)=-A/|A| 1/√(A^2-B^2),   |(B)|+(A)≤ 0,

one obtains from (<ref>) the compact expression

    ⟨ e^-2β W⟩_c,<
    
    =√(2κ(Q^∗_sq-Q^∗))/√([κ-Q^∗-4(Q^∗_sq-Q^∗)]^2-|Q^∗ 2-1|).











With some algebra is easy to show that this result coincides with the high-temperature limit of the quantum characteristic function at μ=2iβ, depicted in Eq. (7) of the main text.    Next we consider the case Q^∗>Q^∗_sq.  From the integral

    1/π∫_-∞^+∞ dx e^A x K_0(B|x|)=1/√(-A^2+B^2),   |(A)|-(B)≤ 0,

and Eq. (<ref>) it is easy to see that ⟨ e^-2β W⟩_c,> takes the same form as ⟨ e^-2β W⟩_c,<, but in the domain Q^∗>Q^∗_sq.



 
 

§ WORK PROTOCOLS

 
 


 The existence of additional domains of divergence in var(e^-β W_ dis) away of adiabaticity is allowed by the PMWF.  We study such behavior for two characteristic work protocols: (i) Sudden quench, Q^∗→ Q^∗_sq, and (ii) a protocol with  dω/dt/ω^2(t) independent of the time of evolution, covering the regime 1<Q^∗≤ Q^∗_sq for different values of the time of driving τ.  On the other hand, protocols involving multiple quenches <cit.> report Q^∗ greater than sudden quench, Q^∗_sq.  Although not treated here, it is worth mentioning that such protocols with Q^∗>Q^∗_sq can lead to divergences in var(e^-β W_ dis) in the high-temperature regime not only for expansion protocols but also for compression.


 We briefly recall that
 
    var(e^-β W_ dis)=G(2iβ)/e^-2βΔ F-1,

where G(μ) is the characteristic function of work.  For the quantum harmonic oscillator

    G(μ)=(1-e^-βħω_0)e^iμħ(ω_1-ω_0)/2 [∑_m,n (e^iμħω_1)^m(e^-(iμ+β)ħω_0)^n P^0→τ_m,n],

where P^0→τ_m,n are the corresponding transition probabilities.  Our numerical approach is based on truncation of the quantum numbers n and m in the latter expression.  Using the generating function method one arrives at the compact expression (to be derived in the next section)

    G(μ)=√(2) (1-e^-βħω_0)e^iμħ(ω_1-ω_0)/2/√(Q^∗(1-e^2iμħω_1)(1-e^-2(iμ+β)ħω_0)+(1+e^2iμħω_1)(1+e^-2(iμ+β)ħω_0)-4e^iμħω_1 e^-(iμ+β)ħω_0).

In particular,

    G(2iβ)→√(2)sinh(βħω_0/2)/√(cosh(βħ(2ω_1-ω_0))-1-(Q^∗-1) sinh(2βħω_1)sinh(βħω_0)).

 For the high-temperature (classical) regime βħω_0,1≪ 1, we have

    G(2iβ) →κ/√((2-κ)^2-4κ(Q^∗-1)),

while for the low-temperature (deep quantum) regime βħω_0,1≫ 1,

    G(2iβ) →exp[βħ(ω_0-ω_1)]/√(1-(Q^∗-1/2) (exp(2βħω_0)-1)).

We now discuss the specific work protocols.

 

 §.§ Sudden quench


  For the case of sudden quench, Q^∗→ Q^∗_sq=(κ+1/κ)/2, the domain of divergence in var(e^-β W_ dis) for the high-temperature regime (classical) is ω_0>√(2)ω_1, as predicted by Eq. (<ref>).  As we depart from the high-temperature limit, an extra domain of divergence emerges from the regime of ω_1∼∞, reaching smaller values with lower temperatures, see panels (b) to (e) in Fig. <ref>.  This new domain divergence originates from quantum effects as they are absent in the limit of ħ→ 0.  We can further digest the divergence by looking for the zero and imaginary values of the denominator in the characteristic function in Eq. (<ref>), i.e.,
     
    cosh(βħ(2ω_1-ω_0))-1-(Q^∗_sq-1) sinh(2βħω_1) sinh(βħω_0)≤ 0.

 To gain useful insights we consider the limit ω_0→ 0.  Using  Q^∗_sq∼ 1/2κ,  Eq. (<ref>) reduces to
  
    2 sinh(βħω_1)-βħω_1 cosh(βħω_1)≤ 0.

Noting the obvious inequality

    2 sinh(βħω_1)-βħω_1 cosh(βħω_1)<(2-βħω_1) cosh(βħω_1),

the existence of a divergence requires

    (2-βħω_1)< 0.

This divergence is not present in the high-temperature (classical) regime, being associated to ħβ→ 0. Panels (c) to (e)  in Fig. <ref> show such boundary extending to ω_0>0, but always within the compression sector, ω_0<ω_1.  As an example, the inequality (<ref>) predicts an onset of divergence around ω_1=0.2 for the temperature ħβ=10, in the limit ω_0→ 0.   This is in agreement with Fig. <ref>(d).














 §.§ A protocol with constant dω/dt/ω^2(t).


The second protocol allows us to study the impact of a varying protocol duration τ, between sudden quench (τ→0) and the slow driving limit (τ→∞).  The explicit protocol is specified by <cit.>

    ω(t)=ω_1ω_0τ/ω_0τ+t(ω_1-ω_0).

Thanks to a time-independent  dω/dt/ω^2(t), one can readily obtain
the corresponding Husimi coefficient

    Q^∗={[  1+cosh(√(1-γ^2))ln(κ)/1-γ^2, if γ^2≤ 1,;                  1+1/2 ln(κ),, if γ^2=1,; 1+1-cos(√(γ^2-1))ln(κ)/γ^2-1, if γ^2≥ 1; ].

where γ= 2ω_0τ/(1-κ), and |dω/dt/ω^2(t)|=2/γ.  Depicted in Fig. <ref> is the behavior of var(e^-β W_ dis) for a time of driving τ=15 and temperatures ħβ=10, 15.












 Compared with the sudden quench case, the divergent domains (see Fig. <ref>(b) and Fig. <ref>(d)) are now fragmented into multiple domains with subtle phase boundaries.  This fragmentation can be partially traced back to the non-monotonic behavior of Q^∗ as a function of ω_0 and ω_1 (see dashed lines in Fig. <ref>(a) and Fig. <ref>(c)).  Indeed, for |ω/dt/ω^2(t)|≤ 2, then Q^∗=1+(1-cos(γ^2-1)ln(κ))/(γ^2-1), where γ=2ω_0τ/(1-κ).


In addition, as the temperature decreases, the domain of divergence tends to grow (the white area in Fig. <ref>(c) is greater than in Fig. <ref>(b)).  In the low temperature limit, the domain for having finite var(e^-β W_ dis) (gray area) shrinks to almost zero area,
leaving only zero-measure boundaries separating different domains of divergence.
Such zero-measure cases with finite var(e^-β W_ dis) arise from  finite-time adiabatic dynamics whose Q^∗=1, which is possible because for |ω/ω^2(t)|≤ 2, Q^∗=1 if  [1-cos(√(γ^2-1))ln(κ)]/(γ^2-1)=0.
As mentioned in the main text, similar disconnected regions of divergence are observed along the time of driving, where convergent domains are found around the adiabatic times <cit.>: τ_n=√(1+4π^2n^2/(lnκ)^2)|1-κ|/(2ω_0), n=1,2,3..., consistent with the solutions for Q^∗=1 in (<ref>).


As shown in Fig. <ref>(a), starting from the slow driving limit in a region of ω_0 and ω_1 where var(e^-β W_ dis) converges, there is a minimum time of driving τ_c above which no divergence is found. Below this time there is an irregular alternation of convergent and divergent behavior.  The local minimas in var(e^-β W_ dis)  are associated with finite-time adiabatic dynamics.   Note also that for small values of ω_0 similar conclusions as in the sudden quench case can be drawn for the present protocol (<ref>) because Q^∗→ Q_sq^∗ as ω_0→ 0.







 
















§ CHARACTERISTIC FUNCTION FOR THE QUANTUM HARMONIC OSCILLATOR UNDER CHANGE OF THE TRAP FREQUENCY







To compute the characteristic function we use the generating function method <cit.>.  The series in Eq. (3) of the main text has the form P(u,v)=∑_m,nu^mv^nP_m,n^τ, where u=e^iμħω_1 and v=e^-(iμ+β) ħω_0.  We review the procedure to arrive at a compact expression for P(u,v).  From the definition of quantum transition probability,

    P^0→τ_m,n   =   |∫ dx_0∫ dx (ϕ^τ_m(x))^∗U(x,x_0;τ)ϕ^0_n(x_0)|^2,

where each integral is along the real line.  The series becomes

    P(u,v)   =   ∑_m,n=0^∞u^mv^n(∫ dx_0∫ dx (ϕ^τ_m(x))^∗ U(x,x_0;τ)ϕ^0_n(x_0))(∫ dy_0∫ dy (ϕ^τ_m(y))^∗ U(y,y_0;τ)ϕ^0_n(y_0))^∗
       =   ∫∫∫∫ dx_0dxdy_0dy (∑_m=0^∞u^mϕ^τ_m(y)(ϕ^τ_m(x))^∗) U(x,x_0;τ)U(y_0,y;-τ) (∑_n=0^∞v^nϕ^0_n(x_0)(ϕ^0_n(y_0))^∗),

where U^∗(y,y_0;τ)=U(y_0,y;-τ).  From Mehler's Hermite polynomial formula, we have


    ∑_m=0^∞u^mϕ^τ_m(y)(ϕ^τ_m(x))^∗   =   √(mω_1/ħπ(1-u^2)) exp(-mω_1/2ħ (1+u^2)(x^2+y^2)-4uxy/(1-u^2)),
    ∑_n=0^∞v^nϕ^0_n(x_0)(ϕ^0_n(y_0))^∗   =   √(mω_0/ħπ(1-v^2)) exp(-mω_0/2ħ (1+v^2)(x_0^2+y_0^2)-4vx_0y_0/(1-v^2)).



The latter equations require (u)^2 and (v)^2 to be lesser than 1.  Using the well known propagators <cit.>


    U(x,x_0;τ)   =   √(m/2π iħ X) exp(im/2ħ X(Ẋx^2-2xx_0+Yx_0^2)),
    
    U(y_0,y;-τ)   =   √(im/2πħ X) exp(-im/2ħ X(Ẋy^2-2yy_0+Yy_0^2));



one arrives at the Gaussian integral

    P(u,v)
       =   √(mω_1/ħπ(1-u^2)) √(mω_0/ħπ(1-v^2)) √(m/2π iħ X) √(im/2πħ X)∫∫∫∫ dx_0dxdy_0dy exp(-  v^ T M  v),

where

    M=m/2ħ([ 1+u^2/1-u^2ω_1 - iẊ/X                   i/X          -2uω_1/1-u^2                     0;                   i/X 1+v^2/1-v^2ω_0 - iY/X                     0          -2vω_0/1-v^2;          -2uω_1/1-u^2                     0   1+u^2/1-u^2ω_1+iẊ/X                  -i/X;                     0          -2vω_0/1-v^2                  -i/X 1+v^2/1-v^2ω_0 + iY/X ])    and     v=([   x; x_0;   y; y_0 ]).

Using (<ref>) and the identity ẊY-XẎ=1, one gets

    det( M)=(m/2ħ)^4 2X^2ω_0ω_1 [Q^∗(1-u^2)(1-v^2)+(1+u^2)(1+v^2)-4uv].

From (<ref>) it follows that <cit.>

    P(u,v)=√(2)/√(Q^∗(1-u^2)(1-v^2)+(1+u^2)(1+v^2)-4uv).

It is worth noting that numerics supports a broader domain of validity to (<ref>) than (u)^2 and (v)^2 to be lesser than 1, as required by the generating function method.  This is easy to see in the case of adiabatic dynamics, where P(u,v)=∑_n(uv)^n and the convergence condition reduces to uv<1, for u,v∈ℝ.





§ WORK FLUCTUATIONS FOR A QUANTUM HARMONIC OSCILLATOR UNDER AN EXTERNAL FORCE





We illustrate how a special work protocol, essentially the one used in the recent proof-of-principle experiment <cit.>,  can lead to total absence of divergence in var(e^-β W_ dis).  In the experiment <cit.>, work is done by applying an external force to an ion while keeping the trapping frequency fixed.   The general driving in this case is given by the time-dependent Hamiltonian

    H(t)=ħω a^†a+f^∗(t)a+f(t)a^†,

where a and a^† are the creation and annihilation operators of a QHO with trapping frequency ω; the complex mechanical parameters are f and f^∗, with initial values f(0)=f^∗(0)=0.  We note that the qubit in <cit.>, modifying (<ref>) by a coupling to the external force (a+a^†) via the Pauli matrix σ̂_x, is not initialized at thermal equilibrium but in one of its eigenstates; thus only the phonon is initialized at thermal equilibrium and the qubit is solely used to apply the corresponding work.  The work characteristic function associated with the Hamiltonian in Eq. (<ref>) is given by <cit.>


    G(μ)=exp(iμ |f(τ)|^2/ħω) exp[(e^iμħω-1)|z|^2] ∑_n=0^∞P_n^0L_n(4|z|^2sin^2(ħωμ/2)),

where L_n(x) are Laguerre polynomials and again we assume an initial state at thermal equilibrium, P^0_n=e^-β E^0_n/Z_0.  The non-adiabatic dynamics is captured by the rapidity parameter

    z(τ)=1/ħω∫_0^τdt ḟ(t)exp(iω t).

 Using again ⟨ e^-2βW⟩=G(2iβ), one finds that G(2iβ) is always finite.

 Indeed,

from Perron's formula 

    L_n(x)=2^-1π^-1/2e^x/2(-x)^-1/4n^-1/4e^2√(-nx)(1+𝒪(n^-1/2)),        x∈ℂ\ℝ_+,


it is clear that the uniform convergence of the series in (<ref>) is dominated by the probability P_n^0 whose exponential is linear in n.  

  The compact expression of (<ref>) is given by <cit.>

    G(μ)=exp(iμ |f(τ)|^2/ħω+(e^iμħω-1)|z|^2-4|z|^2 sin^2(ħωμ/2)/e^βħω-1).


  Notice that evaluation at μ=iβ renders the Jarzynski equality such that ΔF=|f(τ)|^2/(ħω).  Also, consistently with the PMWF the minimum of G(2iβ) with respect to the rapidity is obtained at |z|^2=0.  For this type of work protocol there is no local minima in var(e^-βW_ dis), in contrast with the behavior depicted in Fig. 1 in the main text.  Indeed, using (<ref>), we find

    var(e^-βW_ dis)=e^2|z|^2sinh(βħω)-1,

which does not suffer from divergence at any given β.





§ WORK FLUCTUATIONS FOR A DRIVEN INFINITE SQUARE-WELL POTENTIAL









We present in Fig. <ref> results for an infinite square-well potential, i.e., H_box(t)=-1/2Mp̂^2+V_t(q̂), where

    V_t(q)={[                0, -L(t)/2<q<L(t)/2,;                ∞,        otherwise; ].

 with a constant velocity of compression, dL/dt.  Again, analytical results are not available yet numerics shows the qualitative behavior characteristic of divergence in var[exp(-β W_ dis)].  Although not shown in the figures we confirmed the following observations: (i) benchmark for the numerics with the adiabatic divergence, which can be easily predicted at L_1>√(2) L_0. (ii) For larger velocities we observe a shift of the divergent behavior towards higher temperatures.  We used the reported transition probability from Ref. <cit.>:

    P_m,n^0→τ=|∑_l=1^∞ {2/L_0∫_0^L_0 exp(-iMv/2ħ L_0x^2) sin(lπ x/L_0) sin(mπ x/L_0)dx..
    ..×  exp[-il^2 π^2ħ(L_1-L_0)/2MvL_0L_1] 2/L_1∫_0^L_1 exp(iMv/2ħ L_1y^2) sin(lπ y/L_1) sin(nπ y/L_1)dy} |^2,

where we set M=1 and ħ=1.  The n-th energy level of the system being E^t_n=n^2π^2ħ^2/(2ML^2(t)).   Our numerical calculation is based on truncation of the series in Eq. (<ref>).







As shown in Fig. <ref>, for work protocols with finite var[exp(-β W_ dis)] at high temperatures, divergence will emerge as temperature decreases.  This again confirms our general insights. 



























 

 





unsrt

10

Hanggireview1
Michele Campisi, Peter Hänggi, and Peter Talkner.
Colloquium : Quantum fluctuation relations: Foundations and
  applications.
Rev. Mod. Phys., 83:771–791, 2011; Erratum Rev. Mod. Phys., 83:1653, 2011.

Hanggireview2
Peter Talkner and Peter Hänggi.
Aspects of quantum work.
Phys. Rev. E, 93:022131, 2016.

Jarzynski1
C. Jarzynski.
Nonequilibrium equality for free energy differences.
Phys. Rev. Lett., 78:2690–2693, 1997.

Kurchan1
J. Kurchan.
A quantum fluctuation theorem.
arXiv:cond-mat/0007360, 2000.

Tasaki1
H. Tasaki.
Jarzynski relations for quantum systems and some applications.
arXiv:cond-mat/0009244, 2000.

Dellago2
Christoph Dellago and Gerhard Hummer.
Computing equilibrium free energies using non-equilibrium molecular
  dynamics.
Entropy, 16:41–61, 2014.

Bustamante1
Jan Liphardt, Sophie Dumont, Steven B. Smith, Ignacio Tinoco, and Carlos
  Bustamante.
Equilibrium information from nonequilibrium measurements in an
  experimental test of Jarzynski's equality.
Science, 296(5574):1832–1835, 2002.

Kiang1
Nolan C. Harris, Yang Song, and Ching-Hwa Kiang.
Experimental free energy surface reconstruction from single-molecule
  force spectroscopy using Jarzynski's equality.
Phys. Rev. Lett., 99:068101, 2007.

Ritort1
Anna Alemany, Alessandro Mossa, Ivan Junier, and Felix Ritort.
Experimental free-energy measurements of kinetic molecular states
  using fluctuation theorems.
Nat. Phys., 8:688–694, 2012.

Szabo1
Gerhard Hummer and Attila Szabo.
Free energy reconstruction from nonequilibrium single-molecule
  pulling experiments.
Proc. Natl. Acad. Sci. USA, 98(7):3658–3661, 2001.

Rabbiosi1
F. Douarche, S. Ciliberto, A. Petrosyan, and I. Rabbiosi.
An experimental test of the Jarzynski equality in a mechanical
  experiment.
EPL (Europhysics Letters), 70(5):593, 2005.

Bechinger1
V. Blickle, T. Speck, L. Helden, U. Seifert, and C. Bechinger.
Thermodynamics of a colloidal particle in a time-dependent
  nonharmonic potential.
Phys. Rev. Lett., 96:070603, 2006.

Serra1
Tiago B. Batalhão, Alexandre M. Souza, Laura Mazzola, Ruben Auccaise,
  Roberto S. Sarthour, Ivan S. Oliveira, John Goold, Gabriele De Chiara, Mauro
  Paternostro, and Roberto M. Serra.
Experimental reconstruction of work distribution and study of
  fluctuation relations in a closed quantum system.
Phys. Rev. Lett., 113:140601, 2014.

Vedral1
R. Dorner, S. R. Clark, L. Heaney, R. Fazio, J. Goold, and V. Vedral.
Extracting quantum work statistics and fluctuation theorems by
  single-qubit interferometry.
Phys. Rev. Lett., 110:230601, 2013.

Paternostro2
L. Mazzola, G. De Chiara, and M. Paternostro.
Measuring the characteristic function of the work distribution.
Phys. Rev. Lett., 110:230602, 2013.

Campisi1
Michele Campisi, Ralf Blattmann, Sigmund Kohler, David Zueco, and Peter
  Hänggi.
Employing circuit qed to measure non-equilibrium work fluctuations.
New Journal of Physics, 15:105028, 2013.

Lutz3
Gerhard Huber, Ferdinand Schmidt-Kaler, Sebastian Deffner, and Eric Lutz.
Employing trapped cold ions to verify the quantum Jarzynski equality.
Phys. Rev. Lett., 101:070403, 2008.

Kihwan1
Shuoming An, Jing-Ning Zhang, Mark Um, Dingshun Lv, Yao Lu, Junhua Zhang,
  Zhang-Qi Yin, H. T. Quan, and Kihwan Kim.
Experimental test of the quantum Jarzynski equality with a
  trapped-ion system.
Nat. Phys., 11(2):193–199, 2015.

Jarzynski2
O. Mazonka and C. Jarzynski.
Exactly solvable model illustrating far-from-equilibrium predictions.
arXiv:cond-mat/9912121, 1999.

Zuckerman1
Daniel M. Zuckerman and Thomas B. Woolf.
Theory of a systematic computational error in free energy
  differences.
Phys. Rev. Lett., 89:180602, 2002.

Bustamante2
J. Gore, F. Ritort, and C. Bustamante.
Bias and error in estimates of equilibrium free-energy differences
  from nonequilibrium measurements.
Proc. Natl. Acad. Sci. USA, 100(22):12564, 2003.

Zuckerman2
Daniel M. Zuckerman and Thomas B. Woolf.
Systematic finite-sampling inaccuracy in free energy differences and
  other nonlinear quantities.
Journal of Statistical Physics, 114(5):1303–1323, 2004.

Dellago1
Wolfgang Lechner and Christoph Dellago.
On the efficiency of path sampling methods for the calculation of
  free energies from non-equilibrium simulations.
Journal of Statistical Mechanics: Theory and Experiment,
  2007(04):P04001, 2007.

Ritort2
Matteo Palassini and Felix Ritort.
Improving free-energy estimates from unidirectional work
  measurements: Theory and experiment.
Phys. Rev. Lett., 107:060601, 2011.

Zuckerman3
F. Marty Ytreberg, Robert H. Swendsen, and Daniel M. Zuckerman.
Comparison of free energy methods for molecular systems.
The Journal of Chemical Physics, 125(18):184114, 2006.

Yi1
Seongjin Kim, Yong Woon Kim, Peter Talkner, and Juyeon Yi.
Comparison of free-energy estimators and their dependence on
  dissipated work.
Phys. Rev. E, 86:041130, 2012.

Kirkpatrick1T. R. Kirkpatrick, J. R. Dorfman, and J. V. Sengers.
Work, work fluctuations, and the work distribution in a thermal nonequilibrium steady state.
Phys. Rev. E, 94: 052128, 2016.

Jarzynski3
Christopher Jarzynski.
Rare events and the convergence of exponentially averaged work
  values.
Phys. Rev. E, 73:046105, 2006.

Deng1
Jiawen Deng, Alvis Mazon Tan, Peter Hänggi, and Jiangbin Gong.
Merits and qualms of work fluctuations in classical fluctuation
  theorems.
Phys. Rev. E, 95:012106, 2017.

Gaoyang1
Gaoyang Xiao and Jiangbin Gong.
Principle of minimal work fluctuations.
Phys. Rev. E, 92:022130, 2015.

Singer1
Johannes Roßnagel, Samuel T. Dawkins, Karl N. Tolazzi, Obinna Abah, Eric
  Lutz, Ferdinand Schmidt-Kaler, and Kilian Singer.
A single-atom heat engine.
Science, 352 (6283):325–329, 2016.

italy1 F. Plastina, A. Alecce, T. J. G. Apollaro, G. Falcone, G. Francica, F. Galve, N. Lo Gullo, and R. Zambrini.
Irreversible work and inner friction in quantum thermodynamic processes.
Phys. Rev. Lett., 113;260601, 2014.

genFT1A. E. Allahverdyan.
Nonequilibrium quantum
fluctuations of work.
Phys. Rev. E, 90;032137,2014.

genFT2G. Watanabe, B. Prasanna Venkatesh, P. Talkner, M. Campisi, and P. Hänggi.
Quantum fluctuation theorems and generalized measurements during the force protocol.
Phys. Rev. E , 89:032114,2014.

genFT3 D. Kafri and S. Deffner.
Holevo's bound from a general quantum fluctuation theorem.
Phys. Rev. A, 86:044302,2012.

genFT4 A. E. Rastegin.
Non-equilibrium equalities with unital quantum channels.
J. Stat. Mech., P06016, 2013.

genFT5 T. Albash, D. A. Lidar, M. Marvian, and P. Zanardi.
Fluctuation theorems for quantum processes.
Phys. Rev. E 88:032146, 2013.

Talkner1
Peter Talkner, Eric Lutz, and Peter Hänggi.
Fluctuation theorems: Work is not an observable.
Phys. Rev. E, 75:050102, 2007.

Nikonov
V. V. Dodonov, V. I. Man'ko and D. E. Nikonov.
Exact propagators for time-dependent Coulomb, delta and other potentials.
Physical Letters A 162 (1992) 359–364.

Lohe
M. A. Lohe.
Exact time dependence of solutions to the time-dependent Schrdinger equation.
J. Phys. A: Math. Theor. 42 (2009) 035307.

Jaramillo
J. Jaramillo, M. Beau and A. del Campo.
Quantum supremacy of many-particle thermal machines.
New J. Phys. 18:075019, 2016.

Husimi1
K. Husimi.
Miscellanea in elementary quantum mechanics, ii.
Progress of Theoretical Physics, 9:381, 1953.

Lutz1
Sebastian Deffner and Eric Lutz.
Nonequilibrium work distribution of a quantum harmonic oscillator.
Phys. Rev. E, 77:021128, 2008.

supp For quantum transition probabilities, classical work fluctuations, derivation of G(μ), and additional results on a harmonic oscillator under an exteral force, as well as a particle in a box, see Supplementary Material.
dJE J. Deng, J. D. Jaramillo, P.  Hänggi, and J. Gong, Deformed Jarzynski equality, arxiv 1707:07393.

delCampo1
Mathieu Beau, Juan Jaramillo, and Adolfo del Campo.
Scaling-up quantum heat engines efficiently via shortcuts to
  adiabaticity.
Entropy, 18(5), 2016.

Uzdin1                                                                                                                      
Raam Uzdin, Emanuele G. Dalla Torre, Ronnie Kosloff, and Nimrod Moiseyev.
Effects of an exceptional point on the dynamics of a single particle
  in a time-dependent harmonic trap.
Phys. Rev. A, 88:022505, 2003.



Rezek1                                                                                                                     
Y. Rezek, P. Salamon, K. H. Hoffmann, and R. Kosloff.
The quantum refrigerator: The quest for absolute zero.
EPL (Europhysics Letters), 85:30008, 2009.

Lutz2                                                                                                                
Fernando Galve and Eric Lutz.
Nonequilibrium thermodynamic analysis of squeezing.
Phys. Rev. A, 79:055804, 2009.

adam1                                                                                                             
J. Janszky and P. Adam.
Strong squeezing by repeated frequency jumps.
Phys. Rev. A, 46:6091–6092, 1992.

Deng2
Jiawen Deng, Qing-hai Wang, Zhihao Liu, Peter Hänggi, and Jiangbin Gong.
Boosting work characteristics and overall heat-engine performance via
  shortcuts to adiabaticity: Quantum and classical systems.
Phys. Rev. E, 88:062122, 2013.

Paternostro1
A. del Campo, J. Goold, and M. Paternostro.
More bang for your buck: Super-adiabatic quantum engines.
Scientific Reports, 4:6208, 2014.

Muga1
E. Torrontegui, S. Ibáñez, S. Martínez-Garaot, M. Modugno, A. del Campo,
  D. Guéry-Odelin, A. Ruschhaupt, X. Chen, and J. G. Muga.
Shortcuts to adiabaticity.
Adv. At. Mol. Opt. Phys., 62:117–169, 2013.

Deffner1
Steve Campbell and Sebastian Deffner.
Trade-off between speed and cost in shortcuts to adiabaticity.
Phys. Rev. Lett., 118:100601, 2017.

Ueda1
Ken Funo, Jing-Ning Zhang, Cyril Chatou, Kihwan Kim, Masahito Ueda, and Adolfo
  del Campo.
Universal work fluctuations during shortcuts to adiabaticity by
  counterdiabatic driving.
Phys. Rev. Lett., 118:100602, 2017.

gaoyang2
Gaoyang Xiao and Jiangbin Gong.
Suppression of work fluctuations by optimal control: An approach
  based on Jarzynski's equality.
Phys. Rev. E, 90:052132, 2014.

Quan1
Long Zhu, Zongping Gong, Biao Wu, and H. T. Quan.
Quantum-classical correspondence principle for work distributions in
  a chaotic system.
Phys. Rev. E, 93:062108, 2016.

Rahav1
Christopher Jarzynski, H. T. Quan, and Saar Rahav.
Quantum-classical correspondence principle for work distributions.
Phys. Rev. X, 5:031038, 2015.



vanZon1                                                                                                          
Ramses van Zon, Lisandro Hernández de la Peña, Gilles H. Peslherbe, and
  Jeremy Schofield.
Quantum free-energy differences from nonequilibrium path integrals.
  ii. convergence properties for the harmonic oscillator.
Phys. Rev. E, 78:041104, 2008.




Talkner2                                                                                                           
Peter Talkner, P. Sekhar Burada, and Peter Hänggi.
Statistics of work performed on a forced quantum oscillator.
Phys. Rev. E, 78:011115, 2008; Erratum Phys. Rev. E, 79:039902, 2009.

Quan2
H. T. Quan and C. Jarzynski.
Validity of nonequilibrium work relations for the rapidly expanding quantum piston.
Phys. Rev. E 85:031102, 2012.











