Quantum Work Fluctuations in connection with Jarzynski Equality

TL;DR

This paper investigates the divergence of work fluctuations in quantum and classical systems during non-adiabatic processes, revealing a phase transition-like behavior that impacts the testing of Jarzynski equality.

Contribution

It uncovers the systematic divergence of work fluctuation variance in non-adiabatic protocols for quantum and classical systems, highlighting the need for control fields to verify Jarzynski equality.

Findings

Divergence of variance in non-adiabatic protocols at low temperatures.

Quantum harmonic oscillator shows distinct divergence behavior from classical systems.

Divergence indicates far-from-equilibrium conditions requiring control fields.

Abstract

A result of great theoretical and experimental interest, Jarzynski equality predicts a free energy change $\Delta F$ of a system at inverse temperature $\beta$ from an ensemble average of non-equilibrium exponential work, i.e., $\langle e^{-\beta W}\rangle =e^{-\beta\Delta F}$. The number of experimental work values needed to reach a given accuracy of $\Delta F$ is determined by the variance of $e^{-\beta W}$, denoted ${\rm var}(e^{-\beta W})$. We discover in this work that ${\rm var}(e^{-\beta 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 ${\rm var}(e^{-\beta W})$ at sufficiently low temperatures, markedly different from the classical behavior. The divergence of ${\rm var}(e^{-\beta 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.

Equations

$$G(\mu )=\int dW{e}^{i\mu W}P(W).$$

$$\mathrm{var}(e^{-\beta W_{\mathrm{dis}}})=\frac{G(2i\beta )}{e^{-2\beta \Delta F}}-1.$$

$$G(\mu )=\sum _{m,n}\exp \left[i\mu (E_n^{\tau }-E_m^0)\right]P_{m,n}^{0\to \tau }{P}_m^0,$$

$$G(\mu )=\frac{\sqrt{2}(1-e^{-\beta \hslash {\omega }_0})e^{i\mu \hslash ({\omega }_1-{\omega }_0)/2}}{\sqrt{Q^{\ast }(1-e^{2i\mu \hslash {\omega }_1})(1-e^{-2(i\mu +\beta )\hslash {\omega }_0})+(1+e^{2i\mu \hslash {\omega }_1})(1+e^{-2(i\mu +\beta )\hslash {\omega }_0})-4e^{i\mu \hslash {\omega }_1}{e}^{-(i\mu +\beta )\hslash {\omega }_0}}}.$$

$$G(2i\beta )\to \frac{\sqrt{2}\sinh (\beta \hslash {\omega }_0/2)}{\sqrt{\cosh (\beta \hslash (2{\omega }_1-{\omega }_0))-1-(Q^{\ast }-1)\sinh (2\beta \hslash {\omega }_1)\sinh (\beta \hslash {\omega }_0)}}.$$

$$G(2i\beta )\to \frac{\kappa }{\sqrt{(2-\kappa {)}^2-4\kappa (Q^{\ast }-1)}}.$$

$${\kappa }_c=2\left(Q^{\ast }-\sqrt{Q^{\ast 2}-1}\right).$$

$$G(2i\beta )\to \frac{\exp [\beta \hslash ({\omega }_0-{\omega }_1)]}{\sqrt{1-\left(\frac{Q^{\ast }-1}{2}\right)\left[\exp (2\beta \hslash {\omega }_0)-1\right]}}.$$

$$⟨{\mathrm{e}}^{-2\beta W}⟩=\sum _{m,n}\exp \left[-2\beta (E_n^{\tau }-E_m^0)\right]P_{m,n}^{0\to \tau }{P}_m^0,$$

$$P_{2\nu ,2\mu }^{0\to \tau }$$

$$P_{2\nu +1,2\mu +1}^{0\to \tau }$$

$$Q^{\ast }(t)=\frac{1}{2}\left({\omega }_0{\omega }_1{X}^2+\frac{{\omega }_0}{{\omega }_1}{\dot{X}}^2+\frac{{\omega }_1}{{\omega }_0}{Y}^2+\frac{1}{{\omega }_0{\omega }_1}{Y}^2\right),$$

$$P^c(W)=\frac{\beta {\omega }_0}{2\pi }\int d{p}_{0}d{q}_0{e}^{-\beta {\mathcal{H}}_0(p_0,q_0)}\delta [W-({\mathcal{H}}_{\tau }(p_0,q_0)-{\mathcal{H}}_0(p_0,q_0))],$$

$$P_{<}^c(W)=\beta \sqrt{\frac{\kappa }{2(Q_{sq}^{\ast }-Q^{\ast })}}\exp \left[\beta \frac{\kappa -Q^{\ast }}{2(Q_{sq}^{\ast }-Q^{\ast })}W\right]I_0\left(\beta \frac{\sqrt{Q^{\ast 2}-1}}{2(Q_{sq}^{\ast }-Q^{\ast })}|W|\right)\Theta (W),$$

$$⟨e^{-2\beta W}{⟩}_c=\int dW{e}^{-2\beta W}{P}^c(W).$$

$$e^{-2\beta W}{P}_{<}^c(W)\to \sqrt{\frac{\beta \kappa }{\pi \sqrt{Q^{\ast 2}-1}}}\exp \left[\beta \left(\frac{\kappa -Q^{\ast }}{2(Q_{sq}^{\ast }-Q^{\ast })}-2+\frac{\sqrt{Q^{\ast 2}-1}}{2(Q_{sq}^{\ast }-Q^{\ast })}\right)W\right].$$

$$P_{<}^c(W)=\beta \sqrt{\frac{\kappa }{2(Q_{sq}^{\ast }-Q^{\ast })}}\exp \left[\beta \frac{\kappa -Q^{\ast }}{2(Q_{sq}^{\ast }-Q^{\ast })}W\right]I_0\left(\beta \frac{\sqrt{Q^{\ast 2}-1}}{2(Q_{sq}^{\ast }-Q^{\ast })}|W|\right)\Theta (-W).$$

$$e^{-2\beta W}{P}_{<}^c(W)\to \sqrt{\frac{\beta \kappa }{\pi \sqrt{Q^{\ast 2}-1}}}\exp \left[\beta \left(\frac{\kappa -Q^{\ast }}{2(Q_{sq}^{\ast }-Q^{\ast })}-2-\frac{\sqrt{Q^{\ast 2}-1}}{2(Q_{sq}^{\ast }-Q^{\ast })}\right)W\right].$$

$$\left(\frac{\kappa -Q^{\ast }}{2(Q_{sq}^{\ast }-Q^{\ast })}-2-\frac{\sqrt{Q^{\ast 2}-1}}{2(Q_{sq}^{\ast }-Q^{\ast })}\right)<0$$

$$P^c(W)=\beta \sqrt{\frac{\kappa }{2{\pi }^2(Q^{\ast }-Q_{sq}^{\ast })}}\exp \left[-\beta \frac{\kappa -Q^{\ast }}{2(Q^{\ast }-Q_{sq}^{\ast })}W\right]K_0\left(\beta \frac{\sqrt{Q^{\ast 2}-1}}{2(Q^{\ast }-Q_{sq}^{\ast })}|W|\right),$$

$$e^{-2\beta W}{P}^c(W)\to \sqrt{\frac{\beta \kappa }{\pi \sqrt{Q^{\ast 2}-1}}}\exp \left\{-\beta \left[\left(\frac{\kappa -Q^{\ast }}{2(Q^{\ast }-Q_{sq}^{\ast })}+2\right)W+\frac{\sqrt{Q^{\ast 2}-1}}{2(Q^{\ast }-Q_{sq}^{\ast })}|W|\right]\right\}.$$

$$\frac{\kappa -Q^{\ast }}{2(Q^{\ast }-Q_{sq}^{\ast })}+2+\frac{\sqrt{Q^{\ast 2}-1}}{2(Q^{\ast }-Q_{sq}^{\ast })}<0,\text{for}W>0;$$

$$-\frac{\kappa -Q^{\ast }}{2(Q^{\ast }-Q_{sq}^{\ast })}-2+\frac{\sqrt{Q^{\ast 2}-1}}{2(Q^{\ast }-Q_{sq}^{\ast })}<0,\text{for}W<0.$$

$$\left|\frac{\kappa -Q^{\ast }}{2(Q^{\ast }-Q_{sq}^{\ast })}+2\right|>\frac{\sqrt{Q^{\ast 2}-1}}{2(Q^{\ast }-Q_{sq}^{\ast })}.$$

$$\kappa >2(Q^{\ast }-\sqrt{Q^{\ast 2}-1}).$$

$${\int }_0^{+\infty }dx{e}^{Ax}{I}_0(Bx)=-\frac{A}{|A|}\frac{1}{\sqrt{A^2-B^2}},|\Re (B)|+\Re (A)\le 0,$$

$$⟨e^{-2\beta W}{⟩}_{c,<}=\frac{\sqrt{2\kappa (Q_{sq}^{\ast }-Q^{\ast })}}{\sqrt{[\kappa -Q^{\ast }-4(Q_{sq}^{\ast }-Q^{\ast }){]}^2-|Q^{\ast 2}-1|}}.$$

$$\frac{1}{\pi }{\int }_{-\infty }^{+\infty }dx{e}^{Ax}{K}_0(B|x|)=\frac{1}{\sqrt{-A^2+B^2}},|\Re (A)|-\Re (B)\le 0,$$

$$\mathrm{var}(e^{-\beta W_{\mathrm{dis}}})=\frac{G(2i\beta )}{e^{-2\beta \Delta F}}-1,$$

$$G(\mu )=(1-e^{-\beta \hslash {\omega }_0})e^{i\mu \hslash ({\omega }_1-{\omega }_0)/2}\left[\sum _{m,n}(e^{i\mu \hslash {\omega }_1}{)}^m(e^{-(i\mu +\beta )\hslash {\omega }_0}{)}^n{P}_{m,n}^{0\to \tau }\right],$$

$$G(\mu )=\frac{\sqrt{2}(1-e^{-\beta \hslash {\omega }_0})e^{i\mu \hslash ({\omega }_1-{\omega }_0)/2}}{\sqrt{Q^{\ast }(1-e^{2i\mu \hslash {\omega }_1})(1-e^{-2(i\mu +\beta )\hslash {\omega }_0})+(1+e^{2i\mu \hslash {\omega }_1})(1+e^{-2(i\mu +\beta )\hslash {\omega }_0})-4e^{i\mu \hslash {\omega }_1}{e}^{-(i\mu +\beta )\hslash {\omega }_0}}}.$$

Full text

Introduction. Work fluctuation theorems are a central topic in modern non-equilibrium statistical mechanics Hanggireview1 ; Hanggireview2 . One outstanding result is Jarzynski equality (JE) Jarzynski1 which makes use of an average of non-equilibrium exponential work, i.e., $e^{-\beta W}$ (where $W$ is the work and $\beta$ is the inverse temperature) to obtain important equilibrium information. It takes the form $\langle e^{-\beta W}\rangle=e^{-\beta\Delta F}$, where $\langle\cdot\rangle$ represents an ensemble average of all possible work values with the system initially prepared in a Gibbs state, $\Delta F$ is the free energy difference between initial and final systems at the same $\beta$. 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 Kurchan1 ; Tasaki1 , 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_{\rm dis}=W-\Delta F$. Then JE gives $\langle e^{-\beta W_{\rm dis}}\rangle=1$.

Early applications of JE focused on biomolecular systems, where it is a well-established method to compute $\Delta F$ Dellago2 . Proof-of-principle experiments were carried out Bustamante1 ; Kiang1 ; Ritort1 following an early proposal by Hummer and Szabo Szabo1 . Other experiments include those using classical oscillators Rabbiosi1 ; Bechinger1 . On the quantum side, the first experiment testing JE was first reported for a quantum spin Serra1 , relying on interferometric schemes to reconstruct work distributions Vedral1 ; Paternostro2 . While other settings have been proposed, e.g., circuit QED Campisi1 , a more recent experimental confirmation, following the proposal Lutz3 , uses a trapped ion subject to an external force Kihwan1 . 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 $\Delta F$, as rare events can potentially dominate the average Jarzynski1 ; Jarzynski2 . This in fact has motivated a number of studies in the literature to investigate the errors in predicting $\Delta F$ based on JE Zuckerman1 ; Bustamante2 ; Zuckerman2 ; Dellago1 ; Ritort2 ; Zuckerman3 ; Yi1 ; Kirkpatrick1 . As suggested by the central limit theorem (CLT), the errors are related to the variance in the exponential work, i.e., ${\rm var}(e^{-\beta W})\equiv\langle e^{-2\beta W}\rangle-e^{-2\beta\Delta F}$. The larger ${\rm var}(e^{-\beta W})$ is, the more realizations of $W$ we must collect to reach a given accuracy in $\Delta F$. To obtain $\Delta F$ within an error of $k_{B}T$, the number of $W$ realizations needed was estimated as ${\rm var}(e^{-\beta W_{\rm dis}})={\rm var}(e^{-\beta W})/e^{-2\beta\Delta F}$ Jarzynski3 ; Dellago1 .

The quantity ${\rm var}(e^{-\beta W_{\rm 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, ${\rm var}(e^{-\beta W_{\rm dis}})$ can diverge, suggesting higher-energy components from the initial Gibbs distribution may contribute even more to the variance. Surprisingly, prior to our work Deng1 , this possibility of divergence was rarely mentioned, with one exception Dellago1 treating an adiabatic work protocol. According to the principle of minimal work fluctuations (PMWF) Gaoyang1 ; Deng1 , once ${\rm var}(e^{-\beta W_{\rm 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 $\Delta 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 Deng1 for a specific situation, but in general the error estimate under a diverging ${\rm var}(e^{-\beta W_{\rm dis}})$ is unknown.

This work discovers the divergence of ${\rm var}(e^{-\beta W_{\rm 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 ${\rm var}(e^{-\beta W_{\rm dis}})$ and a diverging ${\rm var}(e^{-\beta W_{\rm 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 $\omega$ because recent experiments to test quantum JE Kihwan1 and to construct quantum heat engines Singer1 were based on QHO (we consider an-harmonic Hamiltonians in the end). Work is done to the system by varying $\omega$. 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 ${\rm var}(e^{-\beta W_{\rm 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 ${\rm var}(e^{-\beta W_{\rm dis}})$. Indeed, even when the classical ${\rm var}(e^{-\beta W_{\rm 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” italy1 alone can induce critical changes in work fluctuations. This work shall also motivate parallel studies on several generalized quantum fluctuation theorems genFT1 ; genFT2 ; genFT3 ; genFT4 ; genFT5 .

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 Talkner1

[TABLE]

For a work protocol starting at $t=0$ and ending at $t=\tau$, $G(\mu)$ can be expressed as a quantum correlation function $G(\mu)={\rm Tr}\left[e^{iuH_{H}(\tau)}e^{-iuH(0)}\rho(0)\right]$, where $\rho(0)$ is the initial Gibbs state and $H_{H}(\tau)$ is the final Hamiltonian in the Heisenberg representation. While the value of $G(\mu)$ at $\mu=i\beta$ yields JE, the value of $G(\mu)$ at $\mu=2i\beta$ determines ${\rm var}(e^{-\beta W_{\rm dis}})$. That is,

[TABLE]

This observation indicates that a divergence in $G(2i\beta)$ is equivalent to a divergence in ${\rm var}(e^{-\beta W_{\rm dis}})$ or in ${\rm var}(e^{-\beta W})$. The quantum correlation function $G(\mu)$ on the imaginary axis is hence our central object. An operative expression of (1) is given by

[TABLE]

where $P^{0\rightarrow\tau}_{m,n}$ represents the transition probability from $m$th state of $H(0)$ to $n$th state of $H(\tau)$. Equations (2) and (3) indicate that, for a system with a finite-dimension Hilbert space such as a spin system, ${\rm var}(e^{-\beta W_{\rm 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 Nikonov ; Lohe . 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 Jaramillo . 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_{\text{QHO}}(\omega(t))=\hat{p}^{2}/2m+m\omega^{2}(t)\hat{q}^{2}/2$. Work is done as the trapping frequency changes from $\omega_{0}=\omega(0)$ to $\omega_{1}=\omega(\tau)$. Thanks to Husimi, all $P^{0\rightarrow\tau}_{m,n}$ are available as a function of the Husimi coefficient $Q^{\ast}\geq 1$ Husimi1 , which is defined as $Q^{\ast}(t)=\frac{1}{2}\left(\omega_{0}\omega_{1}X^{2}+\frac{\omega_{0}}{\omega_{1}}\dot{X}^{2}+\frac{\omega_{1}}{\omega_{0}}Y^{2}+\frac{1}{\omega_{0}\omega_{1}}Y^{2}\right)$, where $X(t)$ and $Y(t)$ are the two solutions of the equation of motion of the corresponding classical oscillator Husimi1 . For adiabatic dynamics, $Q^{\ast}=1$ and $P^{0\rightarrow\tau}_{m,n}=\delta_{m,n}$. Any non-adiabaticity is captured by deviations of $Q^{\ast}$ from its adiabatic value $Q^{\ast}=1$. For later use, we define the compression ratio $\kappa\equiv\omega_{0}/\omega_{1}$. For a sudden quench (sq) protocol $\tau\rightarrow 0$, then $Q^{\ast}\rightarrow Q^{\ast}_{\rm sq}=(\kappa+1/\kappa)/2$.

High-temperature and low-temperature regimes. It is of interest to first examine the behavior of $G(\mu=2i\beta)$ in the strictly adiabatic case, i.e., $P^{0\rightarrow\tau}_{m,n}=\delta_{m,n}$. Because the convergence criteria of the resulting geometric series in Eq. (3) becomes $i\mu\omega_{1}-(i\mu+\beta)\omega_{0}<0$ with $\mu=2i\beta$, it is obvious that if $\kappa\geq 2$, then $G(2i\beta)$ and hence ${\rm var}(e^{-\beta W_{\rm 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\beta)$ under non-adiabatic work protocols, provided that adiabatic work protocols do not yield divergence in $G(2i\beta)$, i.e., in the regime of $\kappa<2$. To that end, we partially resort to a compact expression for $G(\mu)$ Husimi1 . That is,

[TABLE]

Extra care is needed when one extends such a non-analytic result based on integration assuming real $\mu$ to the imaginary axis of $\mu$. The Jarzynski equality is always recovered from (4) at $\mu=i\beta$ Lutz1 , i.e., $G(i\beta)=\frac{\sinh(\beta\hbar\omega_{0}/2)}{\sinh(\beta\hbar\omega_{1}/2)}$. Encouraged by this, we cautiously use this compact expression to help us digest the possible critical boundary of divergence in $G(2i\beta)$ in the regime of $\kappa<2$. Equation (4) then takes us to the following compact expression for $G(2i\beta)$,

[TABLE]

The real denominator of $G(2i\beta)$ in Eq. (5) approaching zero signifies a boundary separating finite values of $G(2i\beta)$ from its diverging values. We stress however, all such “shortcut” conclusions are carefully checked against the original expression in Eq. (3), with $P^{0\rightarrow\tau}_{m,n}$ in the sum series truncated at some large values of $m$ and $n$ (see Supplementary Material supp ).

For the high-temperature (classical) regime $\beta\hbar\omega_{0,1}\ll 1$ we have

[TABLE]

The boundary of divergence in $G(2i\beta)$ is identified at

[TABLE]

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

The behavior of $G(2i\beta)$ is entirely different in the low-temperature (deep quantum) regime $\beta\hbar\omega_{0,1}\gg 1$, where the compact expression for $G(2i\beta)$ becomes

[TABLE]

For strictly adiabatic cases, i.e., $Q^{\ast}=1$, Eq. (8) is well behaved. However, for any work protocol with even slight non-adiabaticity, $Q^{\ast}\neq 1$, one observes that the denominator in Eq. (8) hits zero or becomes imaginary at sufficiently low temperatures. This indicates that ${\rm var}(e^{-\beta W_{\rm dis}})$ must diverge as temperature decreases, so long as the protocol is not strictly adiabatic. The boundary of divergence in ${\rm var}(e^{-\beta W_{\rm dis}})$, within the low-temperature regime, is located at $\beta_{c}=\ln\left(\frac{Q^{\ast}+1}{Q^{\ast}-1}\right)/(2\hbar\omega_{0})$, which does not explicitly depend on $\kappa$. Further, from Eq. (8) we find that as $\beta$ approaches $\beta_{c}$, ${\rm var}(e^{-\beta W_{\rm dis}})$ diverges as $\sim(\beta_{c}-\beta)^{-1/2}$. To understand the divergence, one can actually focus on the contribution to ${\rm var}(e^{-\beta W_{\rm dis}})$ made by the transitions from $2n$-th state of $H(0)$ to the ground state of $H(\tau)$. After applying Sterling’s formula to $P^{0\rightarrow\tau}_{2n,0}$ for large $n$, this contribution is found to scale as $\frac{1}{\sqrt{n}}\left(\frac{Q^{\star}-1}{Q^{\star}+1}\right)^{n}e^{2n\beta\hbar\omega_{0}}$. This predicts the same $\beta_{c}$ as above, because for $\beta>\beta_{c}$, the more rare the initial state is (sampled from the initial Boltzmann distribution), the more contribution it makes to ${\rm var}(e^{-\beta W_{\rm dis}})$. In other words, quantum transitions associated with very negative work values, though exponentially suppressed by $\left(\frac{Q^{\star}-1}{Q^{\star}+1}\right)^{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\beta\hbar\omega_{0}}$ contained in $e^{-2\beta W}$ at sufficiently low temperatures. This important insight has been confirmed by the behavior of classical and quantum deformed Jarzynzki equalities we just proposed dJE .

All these predicted features are checked in numerics, with some examples shown in Fig. 1. There the results are obtained based on Eq. (3) directly as well as from Eq. (5). In the high-temperature regime (small $\beta$), all the plots ${\rm var}(e^{-\beta W_{\rm dis}})$ vs $\hbar\beta$ (in dimensionless units) become flat, in agreement with Eq. (6). This also indicates that the classical ${\rm var}(e^{-\beta W_{\rm 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 ${\rm var}(e^{-\beta W_{\rm dis}})$ approaches zero as temperature decreases, in agreement with Eq. (8). As seen from Fig. 1, local minima in ${\rm var}(e^{-\beta W_{\rm dis}})$ as a function of $\beta$ might emerge, reflecting a competition between quantum fluctuations and thermal fluctuations. In addition, the divergent behavior of ${\rm var}(e^{-\beta W_{\rm dis}})$ close to $\kappa_{c}$ or $\beta_{c}$ (obtained from direct numerics and not shown) is also found to agree with the scaling laws of $(\kappa_{c}-\kappa)^{-1/2}$ and $(\beta_{c}-\beta)^{-1/2}$.

Specific work protocols at intermediate temperatures. To further investigate the behavior of ${\rm var}(e^{-\beta W_{\rm dis}})$, we turn to a specific work protocol where the trapping frequency is suddenly quenched from $\omega_{0}$ to $\omega_{1}$. In this case $Q^{\ast}=(\kappa+1/\kappa)/2$. In the high-temperature regime, Eq. (7) yields a critical $\kappa_{c}=\sqrt{2}$, namely, ${\rm var}(e^{-\beta W_{\rm dis}})$ diverges if $\omega_{0}\geq\sqrt{2}\omega_{1}$. Fig. 2 depicts the numerically obtained domain of divergence in ${\rm var}(e^{-\beta W_{\rm dis}})$ in terms of $\omega_{0}$ and $\omega_{1}$, with panel (a) exactly showing this high-temperature behavior. Quantum effects however dramatically enlarge the domain of divergence (white area) in ${\rm var}(e^{-\beta W_{\rm dis}})$: From panel (b) to (d), temperature decreases and the domain of divergence in ${\rm var}(e^{-\beta W_{\rm dis}})$ gradually invades the (classical) domain of finite ${\rm var}(e^{-\beta W_{\rm dis}})$ in panel (a). For even lower temperatures than shown in Fig. 2, the domain of convergence (gray area) collapse into the line $\omega_{0}=\omega_{1}$, suggesting any actual quench in $\omega$ will yield divergence in ${\rm var}(e^{-\beta W_{\rm dis}})$.

We have also examined a finite-time work protocol, where the parameter $\frac{d{\omega}}{dt}/\omega^{2}(t)$ is chosen to be time-independent, with $\omega(t)=\omega_{1}\omega_{0}\tau/[\omega_{0}\tau+t(\omega_{1}-\omega_{0})]$ delCampo1 . 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 delCampo1 ; Uzdin1 ; Rezek1 for which $Q^{\ast}$ is close to 1. Detailed results are presented in Supplementary Material supp . Finally, it is worth mentioning the possibility of work protocols with $Q^{\ast}>Q^{\ast}_{sq}$, see for example adam1 ; Lutz2 . In these cases ${\rm var}(e^{-\beta W_{\rm dis}})$ may even diverge under compression of the harmonic potential ($\kappa<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_{\text{QHO}}(\omega(t))+J\ {\rm cos}(t\pi/\tau)\hat{q}^{4}$, where $\tau$ is the time of driving, as well as a time-independent quartic potential, with the total Hamiltonian $H_{b}(t)=H_{\text{QHO}}(\omega(t))+J\ \hat{q}^{4}$. There quantum effects are also found to induce the divergence of ${\rm var}(e^{-\beta W_{\rm dis}})$ even if it does not diverge in the classical (high-temperature) domain. Fig. 3 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 ${\rm var}(e^{-\beta W_{\rm 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) Deng2 ; delCampo1 ; Paternostro1 ; Paternostro1 ; Muga1 , to go around this divergence and yet still realizing speedy work protocols (but with some price Deffner1 ; Ueda1 ). Indeed, one can even use ${\rm var}(e^{-\beta W_{\rm dis}})$ as a minimization target to design control fields Gaoyang1 ; gaoyang2 . 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 ${\rm var}(e^{-\beta W_{\rm dis}})$ offers a new angle to study quantum work fluctuations Hanggireview2 ; Quan1 ; Rahav1 and to characterize the too-far-from-equilibrium nature of a work protocol. One fascinating question is to study how the divergence of $G(\mu=2i\beta)$ may be reflected in the behavior of $G(\mu)$ in the real-$\mu$ domain. Finally, one might wonder if the recent experiment testing quantum JE using a trapped ion Kihwan1 already suffered from the divergence issue exposed here. It turns out that this belongs to a fortunate case without any divergence in ${\rm var}(e^{-\beta W_{\rm dis}})$ (see Supplementary Material supp ). Thus, analyzing possible divergence in exponential work fluctuations will be crucial in guiding the design of future quantum experiments testing JE.

**Appendix **

As explained in the main text, all divergences in ${\rm var}(e^{-\beta W})$ are traced back to divergences in

[TABLE]

where $P_{m,n}^{0\rightarrow\tau}$ are the state-to-state transition probabilities and $P_{m}^{0}=e^{-\beta E^{0}_{m}}/Z_{0}$ is the Gibbs distribution of the initial state. All our numeric calculations involve a truncation in the series (9).

Appendix A Transition Probabilities

The transition probabilities between instant eigenstates at $t=0$ and $t=\tau$ for the quantum harmonic oscillator (QHO) with arbitrary driving of the trapping frequency are given by Husimi1

[TABLE]

where, $\nu,\mu\in\mathbb{N}$, and $Q^{\ast}\geq 1$ is the so called Husimi coefficient. The selection rules prevent transitions between energy levels with different parity. The definition of $Q^{\ast}$ is Husimi1

[TABLE]

where $X(t)$ and $Y(t)$ are the two solutions of the equation of motion of the corresponding classical oscillator, $\ddot{X}(t)+\omega^{2}(t)X(t)=0,$ with initial conditions: $X(0)=\dot{Y}(0)=0$ and $\dot{X}(0)=Y(0)=1$. Any non-adiabaticity is captured by deviations of $Q^{\ast}$ from its adiabatic value, $Q^{\ast}=1$. The $Q^{\ast}$ used in the main text is assumed to be $Q^{\ast}(\tau)$.

Appendix B Work fluctuations for the classical harmonic oscillator

For classical systems the work distribution is given by

[TABLE]

where $\mathcal{H}_{0}(p_{0},q_{0})$ is the initial classical Hamiltonian and $\mathcal{H}_{\tau}(p_{0},q_{0})$ represents the value of the final Hamiltonian $\mathcal{H}_{\tau}$ 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^{\ast}\leq Q^{\ast}_{sq}$ and $Q^{\ast}>Q^{\ast}_{sq}$, where $Q^{\ast}$ is the Husimi coefficient defined in (11) and $Q^{\ast}_{sq}$ stands for the Husimi coefficient in sudden quench cases.

We first consider the regime $Q^{\ast}\leq Q^{\ast}_{sq}$, which will be the case if, for example, we have a monotonically increasing or decreasing $\omega(t)$ from $\omega_{0}$ to $\omega_{1}$. This gives rise to positive definite work or negative definite work. For compression of the harmonic trap $(\omega_{0}<\omega_{1})$ the work distribution reads Deng2

[TABLE]

where $\kappa\equiv\omega_{0}/\omega_{1}$ is the compression ratio, $\Theta(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 vanZon1 . The asymptotic behavior $P^{c}(W\rightarrow\infty)$ determines the convergence of

[TABLE]

One can then use the approximation: $I_{0}(z)\rightarrow e^{z}/\sqrt{2\pi z},\ z\rightarrow\infty$, to get the tail

[TABLE]

Evaluation of $\langle e^{-2\beta W}\rangle$ requires us to integrate $e^{-2\beta W}P^{c}(W)$ along $W$ from [math] to $\infty$. A divergence of $\langle e^{-2\beta W}\rangle_{c}$ will take place when the exponent in (15) 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

[TABLE]

Similarly, the tail is

[TABLE]

Since the term $\langle e^{-2\beta W}\rangle$ is an integration along $W$ from $-\infty$ to [math], the possibility of a negative coefficient, i.e.,

[TABLE]

allows divergence. In particular, for adiabatic expansion of the trap ($Q^{\ast}=1,\ \kappa>1$) the criteria of divergence is $\omega_{0}>2\omega_{1}$.

Now consider the regime where the nonadiabaticity parameter $Q^{\ast}$ is even larger than that associated with a sudden quench, i.e., $Q^{\ast}>Q^{\ast}_{sq}$. This is possible if we introduce multiple quenches to the trapping frequency $\omega$. In this case, the integral (12) for compression or expansion becomes

[TABLE]

where $K_{0}(x)$ is the Macdonald function or modified Bessel function of the third kind. This result is consistent with that reported in Lutz1 ; vanZon1 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)\rightarrow\sqrt{\pi}e^{-z}/\sqrt{2z},\ z\rightarrow\infty$, one obtains the tail

[TABLE]

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

[TABLE]

Briefly,

[TABLE]

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

[TABLE]

Classical characteristic function. First we consider the case $Q^{\ast}<Q^{\ast}_{sq}$. From the integral formula

[TABLE]

one obtains from (14) the compact expression

[TABLE]

With some algebra is easy to show that this result coincides with the high-temperature limit of the quantum characteristic function at $\mu=2i\beta$, depicted in Eq. (7) of the main text. Next we consider the case $Q^{\ast}>Q^{\ast}_{sq}$. From the integral

[TABLE]

and Eq. (19) it is easy to see that $\langle e^{-2\beta W}\rangle_{c,>}$ takes the same form as $\langle e^{-2\beta W}\rangle_{c,<}$, but in the domain $Q^{\ast}>Q^{\ast}_{sq}$.

Appendix C Work protocols

The existence of additional domains of divergence in ${\rm var}(e^{-\beta W_{\rm dis}})$ away of adiabaticity is allowed by the PMWF. We study such behavior for two characteristic work protocols: (i) Sudden quench, $Q^{\ast}\rightarrow Q^{\ast}_{sq}$, and (ii) a protocol with $\frac{d{\omega}}{dt}/\omega^{2}(t)$ independent of the time of evolution, covering the regime $1<Q^{\ast}\leq Q^{\ast}_{sq}$ for different values of the time of driving $\tau$. On the other hand, protocols involving multiple quenches Lutz2 ; adam1 report $Q^{\ast}$ greater than sudden quench, $Q^{\ast}_{sq}$. Although not treated here, it is worth mentioning that such protocols with $Q^{\ast}>Q^{\ast}_{sq}$ can lead to divergences in ${\rm var}(e^{-\beta W_{\rm dis}})$ in the high-temperature regime not only for expansion protocols but also for compression.

We briefly recall that

[TABLE]

where $G(\mu)$ is the characteristic function of work. For the quantum harmonic oscillator

[TABLE]

where $P^{0\rightarrow\tau}_{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)

[TABLE]

In particular,

[TABLE]

For the high-temperature (classical) regime $\beta\hbar\omega_{0,1}\ll 1$, we have

[TABLE]

while for the low-temperature (deep quantum) regime $\beta\hbar\omega_{0,1}\gg 1$,

[TABLE]

We now discuss the specific work protocols.

C.1 Sudden quench

For the case of sudden quench, $Q^{\ast}\rightarrow Q^{\ast}_{sq}=(\kappa+1/\kappa)/2$, the domain of divergence in ${\rm var}(e^{-\beta W_{\rm dis}})$ for the high-temperature regime (classical) is $\omega_{0}>\sqrt{2}\omega_{1}$, as predicted by Eq. (31). As we depart from the high-temperature limit, an extra domain of divergence emerges from the regime of $\omega_{1}\sim\infty$, reaching smaller values with lower temperatures, see panels (b) to (e) in Fig. 4. This new domain divergence originates from quantum effects as they are absent in the limit of $\hbar\rightarrow 0$. We can further digest the divergence by looking for the zero and imaginary values of the denominator in the characteristic function in Eq. (30), i.e.,

[TABLE]

To gain useful insights we consider the limit $\omega_{0}\rightarrow 0$. Using $Q^{\ast}_{sq}\sim 1/2\kappa$, Eq. (33) reduces to

[TABLE]

Noting the obvious inequality

[TABLE]

the existence of a divergence requires

[TABLE]

This divergence is not present in the high-temperature (classical) regime, being associated to $\hbar\beta\rightarrow 0$. Panels (c) to (e) in Fig. 4 show such boundary extending to $\omega_{0}>0$, but always within the compression sector, $\omega_{0}<\omega_{1}$. As an example, the inequality (36) predicts an onset of divergence around $\omega_{1}=0.2$ for the temperature $\hbar\beta=10$, in the limit $\omega_{0}\rightarrow 0$. This is in agreement with Fig. 4(d).

C.2 A protocol with constant $\frac{d{\omega}}{dt}/\omega^{2}(t)$.

The second protocol allows us to study the impact of a varying protocol duration $\tau$, between sudden quench ($\tau\rightarrow 0$) and the slow driving limit ($\tau\rightarrow\infty$). The explicit protocol is specified by delCampo1

[TABLE]

Thanks to a time-independent $\frac{d{\omega}}{dt}/\omega^{2}(t)$, one can readily obtain the corresponding Husimi coefficient

[TABLE]

where $\gamma=2\omega_{0}\tau/(1-\kappa)$, and $|\frac{d{\omega}}{dt}/\omega^{2}(t)|=2/\gamma$. Depicted in Fig. 5 is the behavior of ${\rm var}(e^{-\beta W_{\rm dis}})$ for a time of driving $\tau=15$ and temperatures $\hbar\beta=10,\ 15$.

Compared with the sudden quench case, the divergent domains (see Fig. 5(b) and Fig. 5(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^{\ast}$ as a function of $\omega_{0}$ and $\omega_{1}$ (see dashed lines in Fig. 5(a) and Fig. 5(c)). Indeed, for $\left|\frac{\omega}{dt}/\omega^{2}(t)\right|\leq 2$, then $Q^{\ast}=1+\left(1-\cos\left(\gamma^{2}-1\right)\ln(\kappa)\right)/\left(\gamma^{2}-1\right)$, where $\gamma=2\omega_{0}\tau/\left(1-\kappa\right)$.

In addition, as the temperature decreases, the domain of divergence tends to grow (the white area in Fig. 5(c) is greater than in Fig. 5(b)). In the low temperature limit, the domain for having finite ${\rm var}(e^{-\beta W_{\rm dis}})$ (gray area) shrinks to almost zero area, leaving only zero-measure boundaries separating different domains of divergence. Such zero-measure cases with finite ${\rm var}(e^{-\beta W_{\rm dis}})$ arise from finite-time adiabatic dynamics whose $Q^{\ast}=1$, which is possible because for $|{\omega}/\omega^{2}(t)|\leq 2$, $Q^{\ast}=1$ if $[1-\cos(\sqrt{\gamma^{2}-1})\ln(\kappa)]/(\gamma^{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 delCampo1 ; Uzdin1 ; Rezek1 : $\tau_{n}=\sqrt{1+\frac{4\pi^{2}n^{2}}{(\ln\kappa)^{2}}}|1-\kappa|/(2\omega_{0}),\ n=1,2,3...$, consistent with the solutions for $Q^{\ast}=1$ in (38).

As shown in Fig. 6(a), starting from the slow driving limit in a region of $\omega_{0}$ and $\omega_{1}$ where ${\rm var}(e^{-\beta W_{\rm dis}})$ converges, there is a minimum time of driving $\tau_{c}$ above which no divergence is found. Below this time there is an irregular alternation of convergent and divergent behavior. The local minimas in ${\rm var}(e^{-\beta W_{\rm dis}})$ are associated with finite-time adiabatic dynamics. Note also that for small values of $\omega_{0}$ similar conclusions as in the sudden quench case can be drawn for the present protocol (37) because $Q^{\ast}\rightarrow Q_{sq}^{\ast}$ as $\omega_{0}\rightarrow 0$.

Appendix D Characteristic function for the quantum harmonic oscillator under change of the trap frequency

To compute the characteristic function we use the generating function method Husimi1 . The series in Eq. (3) of the main text has the form $P(u,v)=\sum_{m,n}u^{m}v^{n}P_{m,n}^{\tau}$, where $u=e^{i\mu\hbar\omega_{1}}$ and $v=e^{-(i\mu+\beta)\hbar\omega_{0}}$. We review the procedure to arrive at a compact expression for $P(u,v)$. From the definition of quantum transition probability,

[TABLE]

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

[TABLE]

where $U^{\ast}(y,y_{0};\tau)=U(y_{0},y;-\tau)$. From Mehler’s Hermite polynomial formula, we have

[TABLE]

The latter equations require $\Re(u)^{2}$ and $\Re(v)^{2}$ to be lesser than $1$. Using the well known propagators Husimi1

[TABLE]

one arrives at the Gaussian integral

[TABLE]

where

[TABLE]

Using (11) and the identity $\dot{X}Y-X\dot{Y}=1$, one gets

[TABLE]

From (43) it follows that Husimi1

[TABLE]

It is worth noting that numerics supports a broader domain of validity to (46) than $\Re(u)^{2}$ and $\Re(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)=\sum_{n}(uv)^{n}$ and the convergence condition reduces to $uv<1$, for $u,v\in\mathbb{R}$.

Appendix E 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 Kihwan1 , can lead to total absence of divergence in ${\rm var}(e^{-\beta W_{\rm dis}})$. In the experiment Kihwan1 , 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

[TABLE]

where $a$ and $a^{{\dagger}}$ are the creation and annihilation operators of a QHO with trapping frequency $\omega$; the complex mechanical parameters are $f$ and $f^{\ast}$, with initial values $f(0)=f^{\ast}(0)=0$. We note that the qubit in Kihwan1 , modifying (47) by a coupling to the external force $(a+a^{{\dagger}})$ via the Pauli matrix $\hat{\sigma}_{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. (47) is given by Talkner2

[TABLE]

where $L_{n}(x)$ are Laguerre polynomials and again we assume an initial state at thermal equilibrium, $P^{0}_{n}=e^{-\beta E^{0}_{n}}/Z_{0}$. The non-adiabatic dynamics is captured by the rapidity parameter

[TABLE]

Using again $\langle e^{-2\beta{W}}\rangle={G}(2i\beta)$, one finds that ${G}(2i\beta)$ is always finite. Indeed, from Perron’s formula

[TABLE]

it is clear that the uniform convergence of the series in (48) is dominated by the probability $P_{n}^{0}$ whose exponential is linear in $n$. The compact expression of (48) is given by Talkner2

[TABLE]

Notice that evaluation at $\mu=i\beta$ renders the Jarzynski equality such that $\Delta{F}=|f(\tau)|^{2}/(\hbar\omega)$. Also, consistently with the PMWF the minimum of ${G}(2i\beta)$ with respect to the rapidity is obtained at $|z|^{2}=0$. For this type of work protocol there is no local minima in ${\rm var}(e^{-\beta{W}_{\rm dis}})$, in contrast with the behavior depicted in Fig. 1 in the main text. Indeed, using (51), we find

[TABLE]

which does not suffer from divergence at any given $\beta$.

Appendix F Work fluctuations for a driven infinite square-well potential

We present in Fig. 7 results for an infinite square-well potential, i.e., $H_{box}(t)=-\frac{1}{2M}\hat{p}^{2}+V_{t}(\hat{q})$, where

[TABLE]

with a constant velocity of compression, $dL/dt$. Again, analytical results are not available yet numerics shows the qualitative behavior characteristic of divergence in ${\rm var}[\exp(-\beta W_{\rm 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}>\sqrt{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. Quan2 :

[TABLE]

where we set $M=1$ and $\hbar=1$. The $n$-th energy level of the system being $E^{t}_{n}=n^{2}\pi^{2}\hbar^{2}/(2ML^{2}(t))$. Our numerical calculation is based on truncation of the series in Eq. (9). As shown in Fig. 7, for work protocols with finite ${\rm var}[\exp(-\beta W_{\rm dis})]$ at high temperatures, divergence will emerge as temperature decreases. This again confirms our general insights.