Quantum work statistics close to equilibrium

We study the statistics of work, dissipation, and entropy production of a quantum quasi-isothermal process, where the system remains close to the thermal equilibrium along the transformation. We derive a general analytic expression for the work distribution and the cumulant generating function. All work cumulants split into a classical (non-coherent) and quantum (coherent) term, implying that close to equilibrium there are two independent channels of dissipation at all levels of the statistics. For non-coherent or commuting protocols, only the first two cumulants survive, leading to a Gaussian distribution with its first two moments related through the classical fluctuation-dissipation relation. On the other hand, quantum coherence leads to positive skewness and excess kurtosis in the distribution, and we demonstrate that these non-Gaussian effects are a manifestation of asymmetry in relation to the resource theory of thermodynamics. Furthermore, we also show that the non-coherent and coherent contributions satisfy independently the Evans-Searles fluctuation theorem, which sets strong bounds on the statistics, with negative values of the dissipation being exponentially suppressed. Our findings are illustrated in a driven two-level system and an Ising chain, where quantum signatures of the work distribution in the macroscopic limit are discussed.


I. INTRODUCTION
The statistics of work, heat, and dissipation play a central role in the study of the non-equilibrium thermodynamics of small systems, both classical [1,2] and quantum [3][4][5]. They are known to satisfy fluctuation theorems [1,[6][7][8][9], which impose some general restrictions on the form of such distributions. Yet, a complete characterisation of such thermodynamic distributions is highly non-trivial, depending heavily on the details of the specific thermodynamic protocol and the nature (either classical or quantum) of the system under consideration. This complexity contrasts with equilibrium thermodynamics, where it is known that an isothermal reversible transformation outputs a deterministic amount of work, given by the difference of free energy evaluated at the endpoints of the protocol. This observation motivates the study of quasi -isothermal processes, where the transformation is slow enough for the system to always remain close to thermal equilibrium. In this regime, one might hope to find some universal and simple behaviour for the thermodynamic fluctuations. This, for quantum systems, is the aim of the article.
The quasi -isothermal, or slow driving, regime has been thoroughly studied for classical systems [10][11][12][13][14][15]. Notably, the work distribution is known to become Gaussian [11], at least for typical work values [13]. Using Jarzynski equality, this result implies the Fluctuation-Dissipation-Relation (FDR) [6]: Here, σ 2 w ≡ w 2 − w 2 is the variance of the work distribution p(w), w dis the average dissipated work along the protocol, and β = 1/k B T with T the temperature of the environment. This relation means that the work probability distribution p(w) for slowly driven classical systems is completely characterised by average quantities.
When moving to quantum systems this picture appears to be incomplete whenever the driving produces coherences between different energy levels [16]. In fact, it has been shown in [16] that a correction to Eq. (1) is needed in order to account for the fluctuations arising from additional quantum uncertainty in the system. This result also implies that the probability distribution will deviate from a Gaussian whenever coherences are produced, meaning that one needs more information than average quantities alone to characterise p(w).
Starting from this observation, reviewed in section III, a complete study of the irreversibility of quantum systems close to equilibrium is presented. We consider a process in which the Hamiltonian of the system is modified by a sequence of N 1 discrete quenches; after each quench, the system is allowed enough time to thermalise. This kind of protocols, which has been extensively used to describe quasi-isothermal processes [10,12,[17][18][19], is particularly appealing both from the interpretational and the analytical point of view: in fact, providing a way to clearly define work and heat production (being the change of internal energy during the quench or the thermalisation procedure, respectively), we can avoid any reference to the actual equilibration mechanism. Indeed, discrete protocols were introduced in classical thermody-namics as a way to isolate the features of a slowly driven continuous protocol from the details of the relaxation dynamics [10]. In section VII, we show that this intuition carries out to the quantum regime, and we explain how to generalise our method to more general dynamics. Therefore, both for its pedagogical value and the fact that almost no generality is lost, we focus our attention on discrete protocols.
In this context, we are able to identify a number of universal properties of the quantum work statistics close to equilibrium. First we show that the cumulant generating function, and hence all cumulants of p(w), decouple into a classical (non-coherent) and a quantum (coherent) contribution, generalising previous considerations for average quantities [20][21][22][23]. For non-coherent or commuting protocols, only the first two cumulants survive, leading to a Gaussian distribution with its first two moments related through the classical fluctuation-dissipation relation (1). In the presence of quantum coherence, Gaussianity is lost as a consequence of positive skewness and excess kurtosis, witnessing a tendency of the system for extreme deviations above the average dissipation. All such higher work cumulants (beyond the 2nd one) can be obtained by differentiating the quantum skew information (Eq. (31)), a measure of quantum uncertainties [24][25][26]. We also prove that the probability distribution of the dissipation during a protocol equals the one for its time reversed. This fact, together with Crooks relation, leads to the Evans-Searles fluctuation theorem [27], which implies that negative values of the dissipation are exponentially damped (Eq. (26)). These results are explained in section IV, where we also numerically study the validity of the slow driving approximation.
The main technical tool we use is the quasistatic expansion of the cumulant generating function (CGF), see Eqs. (4) and (24). The connection between the CGF and the Fourier transform of p(w) (Eq. (45)) yields a systematic procedure to reconstruct the probability distribution, which we show in section V. In particular, in section V B we consider a protocol in which the eigenbasis of the system Hamiltonian is changed, while keeping the spectrum fixed: we obtain that the probability distribution concentrates on a discrete set even in the quasistatic limit, which is a purely quantum effect reflecting the additional freedom provided by the possibility of creating coherence between energy levels. In section V D we review how the central limit theorem relates to our results, and we study the thermodynamic limit of an Ising chain, showing that despite the Gaussianity of the distribution, the breakdown of the FDR still witness the underlying quantum character of the process.
Finally, in section VI A we identify the origin of the entropy production with the degradation of athermality and asymmetry resources. In particular, we show how in the quasistatic regime the two channels of entropy production decouple at all levels of statistics (Eq. (63,64)). This result also has implications for the resource theory of thermodynamics [28]. We prove that the family of second laws of [29], accounting for the non equilibrium resources in the energy spectrum, collapse into a single law in this regime, which only contributes with a Gaussian term to the probability distribution. On the other hand, we show that the additional constraints on the coherence discussed in [21,30] explicitly account for the non Gaussian effects in the entropy production. These results constitute an additional step in the direction of binding together the quantum work statistics with the resource theoretical description of thermodynamics [31].
Due to the number of results and the technicality of the derivations, many of the discussions and proofs are deferred to the appendices, which are integral part of the work. We always assume bounded spectra and non degeneracy in the energy levels. Dropping both assumptions would not qualitatively change the results, at the price of complicating the exposition. Lastly, a guide to the notations used is given in Appendix J.

II. FRAMEWORK
Consider a thermodynamic protocol where a system is driven while being in contact with a thermal bath at inverse temperature β. In classical thermodynamics the amount of irreversibility produced during the process can be quantified in terms of the Clausius inequality Σ := ∆S −β∆Q ≥ 0, where ∆Q is the heat transferred to the environment, ∆S the change of the system's entropy. This motivates the introduction of the entropy production Σ. Equivalently, one can also define the dissipated work w diss := w − ∆F to be the difference between the work needed to complete the process with respect to the average minimum value given by the free energy change ∆F . The duality between the two formulations is a consequence of the first law of thermodynamics: which connects the two quantities via the relation Σ = βw diss (see e.g. [32] for a recent discussion for quantum systems). This identification should be kept in mind in the rest of the text, as we will often pass from a description in terms of work dissipation rate to one in terms of entropy production rate. We model such thermodynamic processes in the quantum regime by transforming the Hamiltonian of the system by a series of instantaneous quenches between two fixed endpoints H A and H B , and leaving the state relax to thermal equilibrium after each quench. Specifically, we consider protocols performed of N steps, which consists of [17]: 1. Quench on the Hamiltonian: a very fast process in which the Hamiltonian of the system changes as H i → H i+1 , while the state remains unaffected; 2. Equilibration procedure: in which the Hamiltonian is kept fixed, while the system is allowed enough time to perfectly thermalise (ρ i → ρ i+1 ≡ π(H i+1 )).
Since the initial and final point of the process are fixed, increasing the number of steps corresponds to an increasingly slower process. In particular, in the limit N → ∞, the system is always in thermal equilibrium, so that the work is a function of state. We are interested in finite N corrections to this limit, so in the following we will perturbatively expand the quantities of interest over 1/N and keep only first order corrections.
In quantum thermodynamics the work depends on the measurement scheme chosen [33][34][35][36]. We will use here the standard two projective measurement scheme (TPM), which consists in measuring the energy at the beginning and at the end of each quench, and identifying the difference between the two with work [37]. This is justified by the fact that the system is isolated during the quench, so that any change of the internal energy arises from the work performed on the system. From this definition the where we denoted the eigenvalues and eigenvectors of the i-th Hamiltonian by E i and |E i , respectively.
In order to obtain the full probability distribution we can use the fact that the work at each step is an independent random variable, as the thermalisation processes erase any memory of the previous states. Then, the full work distribution p(w) can be obtained by convoluting the work distributions at each step, which is in general an untreatable task. For this reason, it is more convenient to consider the cumulant generating function (CGF), which is additive under independent random processes. The probability distribution p(w) can then be obtained by an appropriate inverse Fourier transform of (4), whereas the cumulants of the work can be directly computed by differentiation of the CGF, i.e.: In the case at hand the CGF is given by (Appendix A): log Tr e −βλHi+1 e βλHi π i .
It is insightful to rewrite Eq. (6) as [38]: where we isolated the contribution coming from the free energy of equilibrium F (H) = −β −1 log Z(H), and we made use of the definition of λ-Rényi divergence S λ ( ||σ) = 1 λ−1 log Tr λ σ 1−λ . In this way we can split the CGF in a deterministic part which is independent on the particular process, and which only shifts the average work by a constant, and a contribution which explicitly depends on the protocol, accounting for the dissipation arising during the process. This splitting also justifies the restriction to the study of the dissipative CGF, defined as: Moreover, this expression highlights the relation between the cumulants of dissipated work and the second laws of thermodynamics [29], as pointed out in [31].
Before continuing, it is worth to point out a number of remarks: (i) it is important to note that the TPM scheme is not invasive in the processes considered; indeed, the first measurement is performed on a diagonal state π i , whereas the second one does dephase the state but this would anyway take place due to the thermalisation process; (ii) the previous observation also implies that the same p(w) would be obtained by other schemes to estimate work such as weak or continuous measurements [39][40][41], so that our results do not depend on the particular measurement scheme used to measure work [34]; (iii) as it was pointed out in the introduction, while we have characterised slow processes by a discrete model, our results can be extended to more general continuous dynamics (e.g., described by time-dependent master equations) as it is sketched in Sec. VII; in this case the derivations and results become more cumbersome so we prefer to keep the more pedagogical and physically insightful discrete model along most of the manuscript.

III. QUANTUM FLUCTUATION DISSIPATION RELATIONS
We can pass to review how the FDR in Eq. (1) changes when passing to quantum systems. These considerations are the natural extension to discrete processes of the work in [16], and constitute the starting point for the study of the probability distribution in the next sections.
We focus on the first cumulants of the distribution: w − ∆F and σ 2 diss = κ (2) w . From (6) and (5), we obtain (Appendix B): where S( ||σ) is the usual relative entropy, and V ( ||σ) is the relative entropy variance, defined as V ( ||σ) := Tr (log − log σ) 2 − Tr [ (log − log σ)] 2 . In the limit in which N 1, both S(.||.) and V (.||.) go to zero as O 1/N 2 in the leading order. Therefore, in this regime one can expand the average dissipation and the work fluctuations as (Appendix B): where we moved to the continuous description H t with t ∈ (0, 1), i.e., we approximated the discrete trajectory (H j with N 1) by a continuous one, see Appendix J for details. The two superoperators in (11) and (12) are defined by: and  11) and (12), this condition means that whenever no coherence in the energy basis is created during the protocol ([H t ,Ḣ t ] = 0 at all times), we have the fluctuationdissipation relation: in complete analogy with the classical case of Eq. (1). As anticipated in the introduction, though, in the general case we obtain a modified FDR of the form where a non-negative quantum correction is present. This takes the form where we have introduced the Wigner-Yanase-Dyson skew information: which is a quantifier of quantum uncertainty in the observable L, as measured in the state [43]. In particular, the skew information is positive I y ( , L) ≥ 0 (it vanishes iff [L, ρ] = 0), decreases under classical mixing, and reduces to the usual variance if ρ is pure. In this sense the breakdown of the FDR (15) is a witness of the presence of purely quantum effects in the protocol, in particular, the creation of coherence between different energy levels. This point is made more precise in section VI A, where we connect the presence of the correction (17) to the creation and degradation of asymmetry during each step of the process.
As explained in [16], the quantum FDR (15) can be interpreted as follows: in a slow process, the work fluctuations σ 2 diss can be expressed as the sum of a thermal contribution (given by 2 w diss /β) and a quantum one coming from the presence of quantum coherence (given by 2Q/β). As an illustration, in Fig. 1 we present the dissipated work and the work fluctuations for a commuting and a non commuting protocol. For lower temperature the two quantities qualitatively differ in the presence of coherence, while in the high temperature limit one always regains the classical picture as the thermal fluctuations dominate.
Interestingly, Eq. (16) also reflects the emergence of non Gaussian behaviour in the work distribution. Indeed, if we take the logarithm of Jarzynski equality we can isolate the first two cumulants and obtain [16]: From Eq. (16), we can identify the sum in the last equation as −βQ, that is: From the properties of the skew information it can be deduced that Q = 0 as soon as coherence is generated at any point along a protocol. This fact, together with Eq. (20), implies that higher-order cumulants must be non-zero. We thus conclude that coherence implies a non-Gaussian work distribution. In fact, we will later prove a stronger statement, namely that the work distribution is Gaussian if and only if Q = 0: In this way we see that in the slow driving regime non-Gaussianity of the work distribution provides a direct witness of quantum coherence. In order to illustrate this point we will first give a closed expression for the CGF, w to the dissipated work, showing how the sum asymptotically converges to the quantum contribution Q. In figure (c) we compare the exact value of Q (dots) with the quasistatic limit presented here (line); in the inset we use the same convention to study the convergence of higher cumulants. All the values have been multiplied by N, and β is set to one. The protocol considered is the same as in (b). (c) Exact numerics study higher cumulants, and then give a characterisation of the probability distribution in the quasistatic regime.

IV. THE CUMULANT GENERATING FUNCTION
In this section we discuss the slow driving approximation (N 1) of Eq. (8). Due to the technicality of the derivations, in order not to over complicate the exposition, we defer the main proofs to Appendix C and D, limiting ourselves here to the qualitative discussion of the results.
The main tool we are going to use is the expansion of the λ-Rényi divergence, given by (Appendix C): where we have defined the y-covariance as: The y-covariance represents a non-commutative generalisation of the classical covariance, reducing to the usual form AB − A B for commuting observables [44]. Using Eq. (22) to expand the sum present in the CGF in Eq. (8), and passing to the continuous limit through the definition of Riemann sums, we finally obtain: where we used the simplified notation cov y t ≡ cov y πt . The y-covariance can be rewritten as the connected two point correlation function betweenḢ t (0) and its time evolved in the imaginary timeḢ t (iβy) (Appendix D). This identification shows the similarity between our approach and linear response theory [45].
As a sanity check for the validity of the approximation, we show in Appendix D that Eq. (24) satisfies both the normalisation condition (K diss (0) = log 1 = 0) and the Jarzynski equality (K diss (1) = log e −β(w−∆F ) = 0). For the derivation of the latter condition, one has to explicitly use the identity: which can be linked in a precise manner with the KMS condition. In this way, one can understand the necessity of a thermal initial state for the Jarzynski equality to hold (Appendix D). The symmetry in the y-covariance Eq. (25) also implies that the CGF satisfies the relation K diss (λ) = K diss (1 − λ). This should be contrasted with the general case, in which K diss (1 − λ) = K diss rev (λ), where we implicitly defined the CGF for the time reversed process. Putting the two conditions together, we can deduce that the probability of having a dissipation w diss is the same both for the forward and the backwards protocol. Then, using Crooks fluctuation theorem, we see that (details in Appendix D): meaning that negative values of the fluctuations are exponentially suppressed. This relation takes the name of Evans-Searles fluctuation theorem, and it will be further analysed in section VI B. In the case in which the protocol does not create coherences ([H t ,Ḣ t ] = 0), the y-covariance reduces to the usual variance, cov y t (Ḣ t ,Ḣ t ) ≡ Var t [Ḣ t ]. Then, one can explicitly carry out the x and y integral in the CGF Eq. (24), which gives: Both Gaussianity and the classical FDR can be directly inferred from this form of the cumulant generating function. In general, since the appearance of finite cumulants of order three or higher is a quantum witness, we can expect that their expression can be directly connected with a measure of coherence, similarly to Eq. (17). This is indeed true: in fact, the CGF can be split in the form: As a consequence, from (5) we see that all cumulants decouple into a classical (i.e., commuting) and a quantum contribution. The connection between this expression and the creation of asymmetry across the protocol will be investigated in section VI A. We can now proceed to the study of the functional form of higher cumulants.

A. Characterisation of higher cumulants
Eq. (28) can be used to give a particularly simple expression for the cumulants. In particular, as it was anticipated in section III, the first two cumulants are given by: in which we see how Q naturally appears in Eq. (29). On the other hand, since the commuting CGF contributes only to the first two cumulants, Eq. (28) highlights how higher cumulants can be expressed in function of the Wigner-Yanase-Dyson skew information alone. In fact, by differentiation we obtain the formula: This shows that all higher moments of p(w) can be directly inferred by taking derivatives of the quantum skew information, a measure of purely quantum uncertainties. Hence from the knowledge of I λ (π t ,Ḣ t ) we can fully obtain all the cumulants of the quantum work distribution in slow discrete processes.
Using the definition of I λ (π t ,Ḣ t ) given in Eq. (18) we can also express the cumulants in a compact form in terms of nested commutators between the Hamiltonian andḢ t (Appendix E): where n ∈ N * , and we recursively defined the family of operators C n by the two properties: C 0 :=Ḣ t , and C n := [H t , C n−1 ]. For example, the first higher cumulants take the simple form: In Appendix E , we show that all the cumulants (odd and even) are positive in the presence of coherence. That is, for protocols with quantum coherence ([H t ,Ḣ t ] = 0 for some t), we have, The general proof of (36) is quite technical being based on a coordinate expression for the CGF. Yet note that from Eq. (33) one can immediately deduce the positivity of all even cumulants. In the case, n = 3, the positivity can be deduced by noting that the skew information is positive for λ ∈ (0, 1), but identically zero for λ = 0 (and hence the first derivative which gives κ (3) w must be positive). From (36) we can infer some information about the shape of the distribution: indeed, the skewness γ 1 and the excess kurtosis γ 2 are connected to the cumulants by the relations: The positivity of κ w and κ (4) w then means that the probability distribution has a fat tail on the right. That is, compared to a normal distribution, values of the dissipation which are bigger than the average by five or more standard deviations are more likely to occur due to quantum fluctuations.

B. Explicit form of the CGF and numerical verifications
Before going further with the analysis, in order to gain some intuitions on the form of the CGF for physical systems we present here the CGF for a two level system and a quantum Ising chain in a transverse field.
We parametrise the Hamiltonian of the two level system by spherical coordinates: H(r, θ, φ) = r cos φ sin θσ x + r sin φ sin θσ y + r cos θσ z ; notice that we can neglect any term proportional to the identity since this would only correspond to a shift in the ground state energy. In this setting, assuming external control over the parameters (r, θ, φ), we can write the CGF as (Appendix F): where we separated the classical contribution and the quantum one, which respectively read: As it can be noticed, the first eigenvalue is associated to a change in the energy spacing only, while the second one corresponds to a change in the eigenbasis orientation. The similarity between Eq. (27) and Eq. (41) is not a coincidence: since changing the energy levels does not create coherence, this part of the protocol will contribute only to the first two cumulants, and will ultimately behave classically.
The second model we consider is an Ising chain in a transverse field, whose Hamiltonian is given by: where we assume control only on the magnetic field h.
Since this model can be mapped into free fermions via a Jordan-Wigner transformation (Appendix G), it is possible to exactly diagonalise it. This allows us to give a closed expression of the y-covariance close to the thermodynamic limit L 1, which reads: where the explicit definition of the function C(k, y, h) is provided in Appendix G. Plugging this expression in Eq. (24) gives the CGF of the model. It should be noticed that thanks to the factor L in Eq. (44), both the y-covariance and the CGF are extensive. At this point we can numerically investigate the validity of the slow driving approximation we used to obtain Eq. (24) from Eq. (8). For this reason, in Fig. 1 we compare the first four cumulants obtained differentiating Eq. (24) and the ones coming from the exact evolution. As it can be seen, the approximation behaves well already for a small number of steps (of the order N ∼ 10). Fast convergences of higher cumulants is also guaranteed by the plot of Q.

V. RECONSTRUCTION OF THE PROBABILITY DISTRIBUTION
The expression for the CGF obtained in Eq. (24) yields a tool to qualitatively characterise the entropy produc-tion distribution, both as a consequence of the symmetries of K diss (λ) (e.g., as we will show with the modified Crooks relations Eq. (26)), and by providing an algorithm to compute the cumulants, which in turn are used to characterise the shape of the probability density.
Additionally to these qualitative results one can use the expression for the CGF to obtain via an inverse Fourier transform the probability distribution. Analytically continuing K diss (λ) to imaginary values (λ ≡ iω/β) gives the identity: Hence, in order to reconstruct the probability distribution it is sufficient to inverse Fourier transform the relation just obtained. We illustrate this procedure in the two opposite limits of a protocol in which no coherence between the energy levels is created and one in which only a change of basis is performed. The qualitatively different results we obtain will be connected with the different origin of the dissipation.

A. Gaussianity of the distribution: commuting protocols
We first consider the simple case in which the protocol does not produce coherence between energy levels. This means that only the commuting part of the CGF Eq. (27) has to be considered. Thanks to the quadratic structure of K diss comm (λ) it is straightforward to perform the inverse Fourier transform, which simply gives: where we defined the average dissipated work as: From this expression one can also directly read off the fluctuation dissipation relation Eq. (15).

B. Non Gaussianity of the distribution: purely coherent protocol
On the opposite limit as the one considered in the previous subsection, we can now pass to study the case in which the process does not affect the spectrum, but it only creates coherences between different energy levels. In order to illustrate this point we will consider here the example of a qubit for whichṙ ≡ 0. Since the eigenvalue e q (λ) only depends on r and not on the coordinates φ and θ, the CGF Eq. (40) takes the form: Figure 2. Illustration of how the slow driving limit affects a process for a classical protocol and a purely coherent one. While in the first case the work output becomes infinitesimal, giving rise to a continuous distribution in the quasistatic limit, for a purely coherent process only the transition probability from one state to the other is affected, while the work output at each step will be either ±∆ or 0. The two protocols are respectively: (a) In this way the time integration reduces to a path dependent constant. For concreteness we choose here to study a protocol in which φ is constant and θ changes as θ = 2πt, for t ∈ [0, 1]. In this case the integral reduces to 4π 2 , and the shape of the distribution is completely characterised by e q only. Writing down the explicit formula for the probability distribution we get: where the CGF reads: Before passing to actually work out the integral in Eq. (48), it is interesting to perform the change of variables ω → ω + π/r. This gives the condition: which tells us that a non zero probability is possible only for w diss ∼ 2nr. In fact, since e q is a periodic function in ω, we can express the exponential in Eq. (48) in terms of the Fourier series: where the factors c n are given by: Plugging this decomposition into Eq. (48), we see that the probability distribution is simply given by a sum of Dirac deltas of the form: this result is compatible with Eq. (50). The distribution is presented in Fig. 2, together with the purely commutative protocol, which leads to a Gaussian distribution.
The persistence of a discrete distribution in the large N limit is a purely quantum effect: in fact, in classical systems there is no way to produce work (or dissipation) without affecting the energy levels. At each step one will produce an energy output of the form {∆ i /N }, resulting in a continuous distribution in the regime N 1 (see Fig. 2 for a illustrative depiction). Contrary, for quantum systems there also exists the freedom of manipulating the system without changing the energy levels, by creating coherence between different eigenvectors. It should be noticed that in this case at each step there is a finite probability of producing an energy output given by ±∆, where ∆ is the spectral gap between the two levels. This quantity does not scale with N , so that even in the limit of infinite number of steps one can only obtain a distribution concentrated on a discrete set. This is illustrated in Fig. 2. This discrete behaviour of the work distribution in slow processes is purely quantum in nature, and the presence of δ functions in the work distribution of a slowly driven system can be understood as a quantum witness.

C. Non Gaussianity of the distribution: qualitative behaviour
For generic transformations of the Hamiltonian in which both the eigenvalues and the eigenvectors are modified, one will obtain work distributions which interpolate between the two extreme cases described above: indeed, p(w) will tend to a continuous distribution but non-Gaussian, with peaks in certain work values coming from quantum fluctuations. To qualitatively describe the deviation from Gaussianity, it is enlightening to look at the sensitivity of the y-covariance on y. We know in fact, from Eq. (27), that if the y-covariance does not depend at all on y, the resulting CGF will be quadratic and, consequently, the distribution will be Gaussian. On the opposite side, we have seen how a non-polynomial dependence of the CGF on λ can lead to exotic behaviour. We will illustrate here how this intuitions can be applied to qualitatively explore the high and low temperature limit for the qubit and the Ising chain.
First, notice that if we plot the y-covariance as a function of β and y (Fig. 3) we can see how the non Gaussianity of the distribution is affected by the temperature. In fact, both for the qubit and for the Ising chain in the high temperature limit the y dependence becomes flat. In this way we regain the expected result that at high temperature the system will behave classically. On the other hand, for higher values of β a non trivial dependence on y manifests, signalling the appearance of quantum effects.
In particular, it is interesting to notice how for the Ising chain the system is more responsive to changes in y for values of the transverse field h ≤ 1. This phenomenon can be understood as a signal of the zero temperature phase transition to a long ordered phase.
Let us now discuss the dependence of e c/q (iω/β) on the temperature, which characterises the work distribution for the driven qubit. First notice that in the high temperature limit we have: witnessing the emergence of a Gaussian behaviour. The fact that both the classical and the quantum contribution converges to the same limit is a consequence of the fact that for β ∼ 0 regardless of the orientation of the eigenbasis or the energy spacing between the two levels the Gibbs state will be given by 1/2. For this reason, [H,Ḣ] ∼ 0 for any change of parameters, making the classical and quantum contribution equal.
If we look at the opposite limit, the case in which β → ∞, we can first notice that: meaning that in the low temperature limit changing the energy spacing will not affect the work distribution. Indeed, most of the population of the system lives in the ground state, so that any manipulation of the excited state will leave the system mostly unaffected. On the other hand, if we look at the same regime for e q a more exotic behaviour emerges. Taking the limit along the imaginary axis we obtain: where the periodicity of the function signals that the distribution becomes concentrated on a discrete set of points, in analogy to what happens for purely coherent protocols. The analysis just presented shows how y-covariances can be used not only to infer statistical properties of a distribution with a level of detail higher than the averages alone, but it can also provide a tool to infer the physics of the underlying system.

D. Gaussianity of the distribution: central limit theorem
The results obtained in the previous sections could seem in contradiction with the central limit theorem: looking at the definition of Q, which, we remind the reader, witness a breakdown of the Gaussianity in the distribution, it is easy to prove that this quantity is extensive. This means that if one considers a system σ made up of L non interacting copies of the same subsystem (σ = ⊗L ), the quantum correction will behave as: and similarly w diss and σ 2 w , so that all terms in the FDR are extensive. This condition, together with Jarzynski equality, implies that the probability distribution of the work output for any finite L, however big, will deviate from a Gaussian distribution.
At the same time though, the central limit theorem says that the standardised sum Σ of L i.i.d. random variables X i , defined by Σ : L, converges in distribution to a Gaussian as L → ∞. For this reason one could be lead to think that Q should converge to zero, due to considerations in the spirit of Eq. (20).
This apparent contradiction comes from an erroneous interpretation of Jarzynski equality: in fact, it should be noticed that it applies only to the work distribution, and not to its rescaled version, which is the one treated by the central limit theorem. For this reason, for any L the work distribution output by σ will deviate from Gaussianity whenever coherences are present. On the other hand, since cumulants of order n are homogeneous of degree n, for the rescaled work output w(σ)/ √ L we have: which is a simple demonstration of the central limit theorem. Using this formula, we have that the rescaled distribution will converge to a normal of the form: with N (µ, σ 2 ) = e −(w−µ) 2 /2σ 2 / √ 2πσ 2 . In this way we see that, even if the Gaussianity of the standardised sum still holds, one can deduce the underlying production of coherence by the breakdown of the classical FDR. In other words, the work distribution of large macroscopic quantum systems (L 1) which are slowly driven will tend to a Gaussian distribution with a larger variance than that predicted by the classical FDR.
In this context, it is also interesting to study what happens to the work distribution of the Ising model in the thermodynamic limit. In fact, if we consider scales larger than the correlation length, the system will behave as the sum of L eff independent copies, meaning that the central limit theorem should hold. Indeed, the dissipative CGF for the rescaled work w/L (where the scaling is chosen so to make the average dissipation finite for L → ∞) takes the form (Appendix G): up to corrections of order O 1/L 2 . We can recognise from this formula the average dissipation and the fluctuations defined in Eq. (11) and Eq. (12). We also recognise the CFG of a Gaussian distribution, with mean independent of L and variance ∝ L −1/2 . If we now take the thermodynamic limit L → ∞, we see that the fluctuations term goes to zero, leaving us with a δ-distribution centred in w diss . This result is a consequence of the equivalence between the canonical and the microcanonical ensemble in the thermodynamic limit [46].

VI. CONSEQUENCES OF THE SLOW DRIVING REGIME
As it was pointed out in the introduction, the characterisation of the entropy production for processes arbitrary out of equilibrium is in general difficult, and it is foreseeable that any universal result won't be sufficient to constrain the statistical properties of the dissipation. On the other hand, we have seen that in the linear response regime the probability distribution of the entropy production can be characterised with relative ease, making reference only to the instantaneous thermal state and to the driving speed. The simplicity of the expressions obtained can be partially ascribed to the decoupling of different channels of the entropy production (section VI A), and to the time reversal symmetry which arises for slow driving protocols (section VI B). These effects and their consequences will be described in the following sections.
A. Channels of entropy production: non Gaussianity and asymmetry The second law of thermodynamics, together with Eq. (2), implies that the entropy production can be interpreted as the deterioration of the ability of a system to perform work. This motivates the study of thermodynamics as a resource theory [28], where in particular one interprets a system out of equilibrium as a resource that can be expended to generate work. This intuition was first investigated by Lieb and Yngvason in [47], where the uniqueness of the entropy functional was proved for equilibrium states.
In the context of the resource theory of thermodynamics, a single law is not sufficient to characterise the irreversibility of a thermodynamic process. In fact, it has been shown in [29] that a necessary condition for the transition → σ to happen through thermal operations is that the following family of second laws is satisfied: Since no work can be extracted from a Gibbs ensemble, the Rényi divergences measure how statistically different a state is from the zero of the theory. The corresponding resource is called athermality and it is progressively lost under thermodynamic evolution.
Moreover, if the state presents off-diagonal terms in the energy eigenbasis an additional family of constraints have to be satisfied by any thermal operations [21,30,48]: where D H is the dephasing map in the H eigenbasis. The corresponding resource is called asymmetry, as it is connected with the breaking of the time translation symmetry of the state. From this set of operations it appears that one cannot increase the coherence between different energy levels with thermal operations alone. Coherence does not come for free, as it could have been guessed from the fact that it can be converted into work in the presence of a coherent bath [49].
We can now pass to analyse the work extraction protocol in this framework. Focussing on a step of the process only, we see that the system starts in a thermal state (which is automatically time symmetric). The quench in the Hamiltonian provides work to the system, which brings the state out of equilibrium and, at the same time, breaks its symmetry: hence, part of the work is converted in athermality, part in asymmetry. Right after the quench a perfectly thermalising operation is applied, which dissipates both resources, bringing the system back to a symmetric equilibrium state.
In general the entropy production Σ can be split in a classical contribution Σ c and a purely coherent one Σ q only on average [21,30]. In this case, the second law is measured by the relative entropy, which corresponds to the limit λ → 1 of the λ-Rényi divergence in Eq. (61,62), where in Eq. (61) is substituted by D H ( ) [21]. Notably this situation changes in the quasistatic regime, and the two channels of entropy production decouple at all levels of statistics. In fact, the y-covariance, which arises in the expansion of the Rényi divergences Eq. (22), can be split in a dephased and a coherent contribution: where we introduced the bookkeeping notationḢ c t := H t − D t (Ḣ t ). This also leads to naturally define the dephased and coherent CGF in the quasistatic regime as: In Appendix H we show how each contribution can be linked with the expansion of terms akin to Eq. (61,62). We can then interpret Eq. (64) as the fact that in the slow driving regime Σ c and Σ q become independent random variables. This situation is exemplified in Fig. 4. In this context, since cov y t (D t (Ḣ t ), D t (Ḣ t )) = Var[D t (Ḣ t )], we can see that the diagonal family of second laws in Eq. (61) collapse into a single constraint, in the spirit of [47]. This result can be thought as a consequence of the adiabatic theorem: in the slow driving regime, in fact, the most probable transitions are of the In the quasistatic regime the entropy production Σ splits in two additive contribution, one which accounts only for the dissipation associated with the athermality created at each step (Σc, red line), and a part coming solely from a change in the energy basis (Σq, orange line).
i+1 , so that the work output becomes quasi-continuous (this discussion should be compared with the one in section V B). Since the only difference between a diagonal quantum state and a classical one is the discreteness of the spectrum, it can be intuitively understood that when this discreteness is smeared out one regains the classical result.
The CGF corresponding to the degradation of athermality takes the form: On the other hand, the CGF for the dissipation of coherence is given by: The splitting in Eq. (64) implies that all the cumulants decompose as κ asymm , which, in particular, means: This result confirms the intuition that the additional channel of entropy production provided by the dissipation of asymmetry raises the average dissipation. Indeed, thanks to the fact that all the cumulants derived from Eq. (66) are positive (section IV), it also holds that κ deph , for any n. More interestingly, the two terms in Eq. (64) independently satisfy the Jarzynski equality, as it can be verified by evaluating Eq. (65) and Eq. (66) for λ = 1, in analogy with what is discussed in section IV. This means that both K diss deph (λ) and K diss asymm (λ) can be considered as arising from two independent thermal processes. The probability distribution of the total dissipated work will then be given by the convolution between the Gaussian coming from the dissipation of athermality resources with the probability distribution coming from the degradation of asymmetry. This observation also implies that the two extremal regimes studied in section V can be considered as the cornerstone of any quasistatic thermodynamic process.
The FDR are also satisfied independently by the two terms. We can then write the inequality: This equation can be read off as the fact that a coherent process dissipates less for the same unit of fluctuation. It should be noticed that, thanks to Jarzynski equality, the presence of a negative correction allows for positive higher cumulants. Therefore, the second law of thermodynamics manifests itself not as a higher dissipation on average, but rather with a fat tail for positive dissipation, that is a tendency of the system to dissipate more than the average. This means that the entropy production associated with the degradation of asymmetry is inherently different than the one associated ot the dissipation of athermality. Notice that the disordered nature of the work output from systems coupled to a coherent bath was already noticed in [49].

B. Time reversal symmetry: the Evans-Searles Fluctuation theorem
As it was pointed out in section IV, from the relation K diss (λ) = K diss (1 − λ) we can deduce the fluctuation relation p(w diss ) = p(−w diss )e βw diss , which typically referred to as the Evans-Searles relation [27,50]. This relation places a considerable constraint on the fluctuations in entropy production, with negative values exponentially suppressed. In fact, one may derive as a direct consequence of of this fluctuation theorem the following: where . + denotes an average over the positive values of entropy production Σ = βw diss and s + = Σ 2 + [51]. One may also obtain a lower bound on the likelihood of negative dissipation: where ψ(.) is the inverse of the function φ(z) = z/(1 + e −z ). These constraints are tighter than those impose.d by applying general concentration bounds that do not make use of the additional information provided by the fluctuation theorem.
The Evans-Searles relation should be compared with the weaker Crooks fluctuation theorem that relates the entropy production to a hypothetical reverse process: Here p rev (−w diss ) represents the probability induced by the time reversed driving H rev i = H 1−t . Comparing the two fluctuation relations we see that the work statistics for the forwards and reverse protocols are indistinguishable: This result has a straightforward interpretation: in fact, we are expanding the dissipated work around its minimum, i.e., the manifold of equilibrium states. This implies that the second laws in the CGF (Eq. (8)) become quadratic forms, so that the system does not differentiate between the drivingḢ t and its reverseḢ rev t = −Ḣ 1−t . If we again consider the splitting of the entropy production Σ = Σ c + Σ q , we see from the symmetry of the y-covariances in (63) that the dephased and coherent contributions actually satisfy the Evans-Searles fluctuation theorem independently: As these terms are independent random variables, we also have The fact that Eq. (74) holds true in this regime provides an interesting connection to the so-called thermodynamic uncertainty relation (TUR) [52][53][54]. The TUR imposes a trade-off between the noise-to-signal ratio S(x) = ∆x/ x of a given time-integrated current and the average entropy production Σ along a process. For small average currents x 1, it has been shown that any distribution of the form Eq. (74) leads to the TUR [55]: The TUR implies that reducing the noise-to-signal ratio associated with a given current comes at a price of increased entropy production. This is consistent with the trade-off relation (68) for currents x = {w deph diss , w asymm diss , }. Notably the TUR is saturated by a Gaussian distribution, which is satisfied for the dephased work x = w deph diss . On the other hand, any non-zero quantum contribution x = w asymm diss cannot saturate the TUR, as clearly seen in (68). This highlights the fact that coherence provides a fundamental limitation to the trade-off between reversibility (small Σ ) and constancy (small S).

VII. ENTROPY PRODUCTION FOR CONTINUOUS PROTOCOLS
Historically, the motivation to analyse step equilibration processes was to construct a framework in which effects arising from the slow driving regime could be isolated from the particular details of the relaxing dynamics [10]. In fact, a process constituted by a discrete sequence of quenches and thermalisations steps can be thought as the simplification of a continuous open-system process in which the equilibration is trivial. We will here motivate this claim and connect the CGF for a continuous protocol with the one we obtained for a discrete one.
In particular, we focus on the regime in which the state of the system can be approximated at all times as t = π t + δ t , where δ t is of order O (1/τ ), and τ is the duration of the protocol. One can expect this kind of behaviour whenever the dynamics is relaxing and the parameters are changed in a sufficiently slow manner.
The cumulant generating function for a continuous process which initially starts at equilibrium is given by [31,38]: For what follows, we denote the time evolved state as τ := U τ π 0 U † τ . Using this notation, we can give a compact expression of the dissipative CGF, which takes the form: This equation shows the close relation between the cumulants of the dissipated work and the statistical difference of the evolved state τ from the thermal state [38]. In the case of a quench the state is given by τ = π 0 , and Eq. (77) reduces to the expression for a single step of Eq. (8). In this sense, in a discrete process the dynamics, which can be assumed to be given by a generic thermalising map, is not taken into account by the CGF in any other way than in the information about the initial conditions i ≡ π i . The dependency of τ on the particular protocol is implicit in this expression. Knowing that the dissipation is path dependent, though, it is useful to highlight this dependency by rewriting Eq. (77) as: We show in Appendix I that this takes the form: where we can already see that both terms are of order O (1/τ ) in the quasistatic limit. Notice that this expression is exact and it directly connects the cumulant generating function with the particular trajectory in the parameter space taken during the protocol.
We can now pass to the slow driving limit of Eq. (79). This means that we assume the state to be given by the approximation t = π t + δ t , and that we can neglect terms of order O δ 2 . Then, Eq. (79) can be simplified to the form (Appendix I): Note that in the above derivation, the adiabatic term π t need not be thermal.
Let us now imagine that the system is embedded in a larger Hilbert space H → H ⊗ H R through coupling to a reservoir at inverse temperature β, while only the local system Hamiltonian is varied in time. We restrict attention to situations where the coupling is weak enough so that we may approximate the global state as a tensor product ρ t ⊗ π R π t ⊗ π R + δρ t ⊗ π R , where π R is a fixed equilibrium state of the reservoir. The y-covariance is then invariant under partial trace over the reservoir, and we end up with the expression (80) for the open system. We note that while neglecting the influence of correlations on the y-covariance is a rather strong assumption, this approximation is sufficient for illustrating the connection between the discrete protocols adopted in the previous sections and continuous-time dynamics. We leave a more rigorous treatment of the CGF for general open quantum systems for future work.
From here we can consider Lindblad evolutions; if the reduced dynamics of the system is given by a time dependent relaxing Lindbladianρ t = L t (ρ t ), then the correction term in the slow driving regime takes the where the cross denotes the Drazin inverse [56][57][58]. Under the assumption of quantum detailed balance [59], we show in Appendix I that the CGF (80) becomes where τ y eq (t) : represents a quantum generalisation of the integral relaxation timescale introduced in [60]. Here we denote the operatorḢ t (ν) = e νL † t [Ḣ t ] evolved in the Heisenberg picture. This timescale quantifies the time over which the fluctuations in power, as quantified by the y-covariance, decay to their equilibrium values. As an example, for a simple Lindbladian of the form L t (ρ t ) = (π t − ρ t )/Γ eq (t), the integral relaxation time reduces to a single timescale τ y eq (t) = Γ eq (t).
From this formula the expression for the average dissipated work and the work fluctuations presented in [16] can be obtained in a straightforward manner, extending the work therein to arbitrary cumulants. Furthermore, from (81) we now find a connection between the continuous protocol approach and the discrete protocol described by the CGF (24). If we identify the ratio between the integral relaxation time and total time with a uniform step size N , namely τ y eq = 1/2N , we find that the two protocols are described by the same statistics. In essence, we may therefore view the discrete protocol as indistinguishable from that of a continuous process γ described by the trivial relaxing Lindbladian L t (ρ t ) = 2N (π t − ρ t ).

VIII. DISCUSSION, CONCLUSIONS AND OUTLOOK
In this paper, we have characterised the fluctuations of work, dissipation, and entropy production in quantum quasi-isothermal processes, where the system of interest stays close to equilibrium along the thermodynamic transformation. In this regime, all cumulants of the work distribution decouple into a classical (non-coherent) and a quantum (coherent) contribution, hence extending previous considerations for average quantities [20][21][22][23]. In fact, all cumulants beyond the second one have a purely quantum origin, and they can be obtained by differentiating the quantum skew information (Eq. (31)), a measure of quantum uncertainities [24,25]. Such quantum fluctuations lead to positive skewness and excess kurtosis, witnessing a tendency of the system for extreme deviations above the average dissipation. These results shed new light on our understanding of quantum features in the work distribution [26,33,34,39,40,[61][62][63][64][65][66].
It is also important to comment on the deviations from Gaussianity in the tails of the work distribution that have been observed for classical processes [13]. These deviations can be understood here by the impossibility of expanding in terms of 1/N work values of order ∼ O (N ) (in a single step); these contributions will certainly appear in unbounded spectra, albeit being extremely unlikely. This problem does not appear for quantum systems with bounded spectra (hence there is always an N sufficiently large as compared to all possible work values during a single step). In this case we can build a strong relation between non-Gaussianity and quantumness: a work distribution becomes non-Gaussian in the quasi-isothermal regime if and only if quantum fluctuations appear. For unbounded spectra, this statement remains true for typical values of the work, while on the tails the "only if" ceases to be valid.
The thermodynamic quantum signatures reported in this article (breakdown of the classical fluctuationdissipation relation Eq. (1), non-Gaussianity of the distribution, with positive skewness and kurtosis) appear to be measurable with state-of-the-art technologies. Indeed, it is sufficient to have an experimental platform where the following three operations can be implemented: (i) projective measurement of the energy for non commuting Hamiltonians of the system {H i }; (ii) quenching of the system Hamiltonian from H i to H i+1 , and (iii) thermalisation of the system for the Hamiltonians H i+1 . Quantum signatures will then be observed whenever [H i , H i+1 ] = 0 for at least some time step i. Ion traps provide excellent controllability and have previously been used to measure quantum work distributions [67,68], and similarly for NMR systems [69]. Other promising platforms for measuring such quantum signatures in the work distribution are quantum dots and superconducting qubits [70][71][72].
From the observation that the entropy production distribution for a particular process equals the one for its time reverse, we have also proven the Evans-Searles relation [27], a stronger form of Crooks fluctuation theorem [50]. This relation enables us to set strong constrains on the fluctuations in entropy production, with negative values exponentially supressed. This result has also enabled us to make a connection with (quantum) Thermodynamic Uncertainty Relations [52][53][54][55], which in this context set a tradeoff between fluctuations and average entropy production. We have then shown that quantum coherence prevents the saturation of the TUR in the slow driving regime.
Our results have also implications for a seemingly unrelated question, namely the interconvertibility of states within the resource theory of thermodynamics [28,73]. Indeed, as it is also observed in [31], there is a close connection between work statistics and the second laws of thermodynamics of [29], both being expressed through Rényi divergences. The expansions of Rényi divergences close to thermal equilibrium developed here imply that the continuous family of second laws of [29] for the interconvertibility of diagonal states reduces to a single one (the second law of thermodynamics) close to thermal equilibrium. This is in spirit similar to the wellknown fact that in the many-copies limit such second laws converge to a single one [73], but this result holds at the single-copy level. On the other hand, the additional constraints from quantum coherence [21,30] remain untouched close to thermal equilibrium. In other words, the fact that the work distribution becomes Gaussian close to equilibrium for diagonal states translates into the many laws of [29] reducing to a single one; whereas the fact that the work distribution is non-trivial for quantum coherent states implies that the asymmetry restrictions of [21,30] do not simplify close to equilibrium. We believe our results might find other implications in the resource theory of thermodynamics; for example, from the expansions of the relative entropy and the relative entropy of variance it appears that the conversions of thermal resources close to equilibrium will always be resonant, in the sense developed in [74].
We have also discussed how the quantum signatures in the work distribution behave for macroscopic systems. The first important observation is that the quantum cor-rection Q in the quantum fluctuation dissipation relation (FDR), β 2 σ 2 w = w diss + Q, is extensive (and so are the dissipation and the variance) which means that no matter how large is the system under observation, a correction to the classical FDR in Eq. (1) will appear for protocols where [Ḣ t , H t ] = 0, hence causing the work distribution to become non-Gaussian. Note that non-commutativity is in fact ubiquitous in many-body quantum systems, where usually the (local) control does not commute with the global Hamiltonian. We also discussed how this observation relates with the central limit theorem, which implies that the rescaled work distribution will converge to a Gaussian, at least for non-interacting or locally interacting many-body systems (away from a phase transition). For such a rescaled distribution, the quantumness in the distribution is encoded in a larger variance of the distribution due to the presence of Q. We have illustrated these considerations in a Ising chain in a driven transverse field, where we have computed the average and variance of the quantum work distribution. These considerations contribute to recent efforts to characterise the work distribution of many-body systems [75][76][77][78][79][80][81][82][83].
While most of our results have been derived through a discrete model of quasi-isothermal processes, we have shown in Sec. VII that our considerations can be naturally extended to more complex continuous dynamics. It remains as an interesting future question to derive these same results by means of a quantum jump approach, or quantum trajectories, given a Lindblad master equation [84,85], a direction that we are currently exploring. Another interesting complementary question is to derive similar slow-driving expansions through linear-response theory for higher cumulants of the work distribution, hence extending previous results for average dissipation [86,87] (see also [88] for a discussion of the work moments beyond weak coupling). In all such extensions, we expect that the results reported for the discrete case should be recovered in the simplest model of an exponential relaxation to equilibrium with a well-defined time-scale, as argued in Sec. VII.
Appendix A: Work probability in the two point measurement scheme In this section we explain how one can arrive to the expression in Eq. (6) for the cumulant generating function of the work starting from the definition of the work probability in the two point measurement scheme (TPM) for discrete processes.
In order to obtain the work output after N steps, we would need to do a convolution between the probability distributions at each step, which is clearly an untreatable task. On the other hand, though, the cumulant generating function for the sum of independent variables factorises in the sum of the CGF of each variable. For this reason, it is useful to introduce the CGF for the full work as: Plugging in the definition of the probability of the work given in Eq.
(3), we then get: log Tr e −βλHn+1 e βλHn n . (A2) where we obtained the last equation by a resummation of the spectral decomposition of the exponential operators. This concludes the derivation of Eq. (6).
Appendix B: Average and fluctuations of dissipated work in the slow driving limit In this section we give the derivation of the slow driving approximation to average dissipated work and work fluctuations. The results of this appendix were already presented in [16,58], and are reproduced here due there prototypical nature, which will be encountered in most of the derivations of the paper.
We first derive Eq. (9) and Eq. (10) from the cumulant generating function, using the definition of cumulants Eq. (5). For what regards the average dissipated work we obtain: in a similar fashion, the work fluctuations can be obtained as: We can recognise in Eq. (B1) and Eq. (B2) the definition of relative entropy and relative entropy variance given in the main text.
Before passing to derive the slow driving limit of the above expression, it is useful to explain how to Taylor expand π i+1 around π i for small variations of the Hamiltonian ∆H i := (H i+1 −H i ). Using the Dyson series for the exponential operators we obtain [89,90]: (B3) One can recognise in the expansion the operator J e −βH defined in the main text (Eq. (13)). Plugging Eq. (B3) in the definition of partition function, one can also prove the equality: where we used the cyclicity of the trace. Then, it is straightforward to see that the Gibbs state π i+1 can be expanded as: where the average in the operator ∆ A := (A − Tr [ A]) comes from the expansion of the partition function in the denominator of the Gibbs state.
Using matrix analysis one can also show that [89,90]: where is definite positive, σ is a traceless perturbation, and J −1 is the inverse superoperator of the one appearing in the Dyson series, which arises from the Taylor expansion of the logarithm. Plugging π i+1 in this expression we then obtain: , which holds up to corrections O ∆H 2 i = O N −2 . In the limit N 1, if the change in H i is smooth enough, one can define an interpolating curve H t and convert all the sums in integrals. For this reason we use the substitution ∆H i → 1 NḢ t , so to rewrite Eq. (B6) as: where we used the definition of Riemann sum. This concludes the derivation of Eq. (11).
We can now pass to expand Eq. (B2). Since we want only term up to order O 1/N 2 we can ignore the second part of the sum. Then the expansion takes the particularly simple form: where to pass from the first to the second line, we used the Taylor expansion of the logarithm [89]. Taking the continuous limit completes the derivation of Eq. (12), which also holds up to order O 1/N 2 .
Appendix C: Expansion of λ-Rényi divergence In this section we prove the following expansion for the λ-Rényi divergence: where σ is a traceless perturbation, J −1 is the inverse of the operator in Eq. (13), and we defined the y-covariance as: Proof. The two main ingredients to prove Eq. (C1) are the Dyson series of the exponential [89]: and the Taylor expansion of the logarithm around the positive definite matrix [90]: Moreover, it is useful to rewrite the λ-Rényi divergence as: where we simply differentiated the trace and integrated again, and we used the fact that Tr λ σ 1−λ | λ=0 = 1. This expression highlights the dependence of the λ-Rényi divergence on the difference between and σ.
Passing to the expansion of S λ ( + εσ|| ), we first focus on the trace inside of the logarithm in Eq. (C5): The first order contribution is given by: where we used the hermiticity of J −1 and the identity J −1 [A] = −1 A, whenever [A, ] = 0. This last relation comes from the fact that the Fréchet derivative of the logarithm agrees with its usual derivative on the subspace of commutative operators.
Passing to the study of the second order contribution, it should be noticed that one can use the change of variables y → 1 − y to obtain the relation: Moreover, the presence of the trace allows us to simplify the double integration as: where passing from the first line to the second we used the substitution u = x 1 and v = x 1 − x 2 , and in the last equation the substitution x = 1 − u and y = 1 − v. Adding together the last two equalities in Eq. (C9), we obtain: Finally, we can pass to expand the logarithm which gives the result at second order: which concludes the proof of Eq. (C1).
Symmetry in the arguments. The quadratic structure of Eq. (C1) hints at the approximate symmetry of the λ-Rényi divergence. In fact, using the identity S λ ( ||σ) = −λ λ−1 S 1−λ (σ|| ), we also obtain that: where the first integral in the second line goes to zero thanks to the symmetry of the y integration bounds, and we performed the change of variables x → 1 − x in the second integral. This proves that the λ-Rényi divergence is symmetric at second order. Relative entropy. Moreover, it should be noticed that in the limit λ → 1 we regain the known expression of Eq. (B5) for the expansion of the relative entropy: where in the first line we used L'Hospital's rule for limits of the form 0/0, and we inserted in the last equation the definition of J .
Appendix D: Cumulant generating function in the slow driving limit In this section we will give the slow driving approximation to Eq. (8), and we will further discuss the remarks made in sec. IV.
Derivation and general results. Starting from the expression of the dissipative CGF and applying the expansion in Eq. (C1) we obtain: where we defined∆H i := π i+1 − π i . As explained in Appendix B, in the limit in which the variation of H i is smooth, one can use the expansion of the Gibbs state in Eq. (B4) to express∆H i , and consequently pass to the continuous limit using the definition of Riemann sum: The dissipative CGF so obtained can be given in coordinates, in the instantaneous energy eigenbasis, as: where the subscript t has been dropped for notation simplicity. This form makes particularly evident the fact that non Gaussian effects can arise solely from off diagonal terms in the driving. Indeed, defining the Wigner-Yanase-Dyson skew information: (which provides a measure of the amount of information contained in with respect to a non commuting observable L [43]), we can rewrite the cumulant generating function as: where in the last equation we explicitly carried out the x and y integration, because the variance does not depend on y. Since a Gaussian distribution has only a quadratic CGF, we see that the non Gaussian contribution can be linked in a precise manner with the lack of commutativity between H t andḢ t . This point will be further analysed in the next sections.
We can now pass to verify some properties of the CGF so obtained. First, it is straightforward to check from Eq. (D2) that both normalisation (K diss (0) = log 1 ≡ 0) and Jarzynski equality (K diss (1) ≡ 0) are preserved by the approximation. The latter can be shown by explicitly writing the integral: and by noticing the fact that the y integration bounds are antisymmetric around x = 1/2, while the y-covariance behaves as cov y t (Ḣ t ,Ḣ t ) = cov 1−y t (Ḣ t ,Ḣ t ). These results are trivial, since the original expression in Eq. (8) satisfies both conditions, but they work as a sanity check to guarantee that after the approximation we still have a probability distribution arising from a thermodynamic process.
y-Covariance and KMS condition. It is interesting to see how Eq. (D2) can be connected to linear response theory. Defining the Heisenberg picture for an operator A by A(t) := e iHt Ae −iHt , and the thermal average as A π := Tr[πA], we can connect the y-covariance with the autocorrelation function of the power operatorḢ t as: This result is analogous in spirit to linear response theory [45]. In fact, Eq. (D8) connects the work response arising from a linear perturbation with the connected two point correlation function of the power operatorḢ t . In this context, the KMS condition reads: Since we used this property to prove Eq. (D7), this expression shows the close connection between the thermality of the state, encoded by the KMS condition, and Jarzynski equality.
Entropy production rate and time reversal symmetry. The entropy production rate is connected with the breaking of time reversal symmetry between a trajectory and its inverse, as it is illustrated by fluctuation theorems [91]. We will show how this condition is encoded in the dissipative cumulant generating function.
Defining the time reversed protocolπ i := π N −i , we get the identity: where we implicitly defined the time reversed CGF. The relation just obtained can be rewritten in terms of the probability distribution of the dissipated work (which, with an abuse of notation, we denote by p(w)): Applying the inverse Fourier transform to the first and the last equation we obtain: In this way we see that Eq. (D10) is equivalent to the Crooks fluctuation relation [91]: which signals the fact that the ratio between the probability of dissipating an amount of work w and the one of getting this work back by reversing the transformation is exponentially big in w. In this way, we can interpret Eq. (D13) as the underlying explanation for the emergence of the arrow of time.
In the slow driving regime the CGF satisfies the condition: where we first used the substitution x → 1 − x, and then in the second line we applied Jarzynski equality to get rid of the first integral in dx. Comparing this relation with Eq. (D10), we can deduce that the probability distribution for the entropy production during a protocol in the slow driving regime equals the one for its time reversed. This could also be inferred from the quadratic structure of the y-covariance, since it does not distinguish betweenḢ t and −Ḣ t . In this context, the Crooks relations becomes: also known as the Evans-Searles relations, which tell us that the probability of having a negative dissipation is exponentially suppressed.

Appendix E: Computation of the cumulants
In this section we explicitly derive a formula for all the cumulants of the distribution starting from Eq. (28): First, we obtain again the expression for the average dissipated work : where we can recognise the second term in the equality as the quantum correction Q defined in Eq. (17). The work fluctuations on the other hand are given by: where we used the fact that the skew information satisfies I λ (π t ,Ḣ t ) = I 1−λ (π t ,Ḣ t ), together with the identity I 0 (π t ,Ḣ t ) = 0. From the third cumulant onwards only the Wigner-Yanase-Dyson skew information contributes to the expression of the cumulants. In fact, if we further differentiate Eq. (E3) , we can see that: Using Eq. (D5) we can give a recursive formula of the equation just obtained. First, it is useful to give the expression for the derivative of the functional: where we used the fact that log π t = −βH t − log Z t . Then, by applying iteratively applying Eq. (E5) and Eq. (E6), we can prove by induction the formula: where the index n runs on integer values starting from n = 1, and we recursively defined the operators C n by the relation: C 0 =Ḣ t and C n = [H t , C n−1 ]. From Eq. (E8) we can infer that all the even cumulants are positive. Moreover, since I 0 (π t ,Ḣ t ) = 0 and I ε (π t ,Ḣ t ) ≥ 0 for any ε ∈ (0, 1), so that from Eq. (E4) we can also deduce that κ (3) w ≥ 0, with equality iff [H t ,Ḣ t ] ≡ 0 at all times. Finally, considering the coordinate expression of the CGF, we can also prove the positivity of the higher odd cumulants. First, it is useful to rewrite Eq. (D4) as: which can be verified by expanding the hyperbolic functions in terms of exponentials. In this way we can obtain even and odd cumulants from the expansion of the hyperbolic cosine and sine, respectively. Then, it is straightforward to give the explicit formula: Since the logarithm preserves the order, the sum is negative. The additional minus in front, then, implies the positivity of κ where J sets the energy scale, h is the intensity of the transverse field and L is the number of sites. We can apply a Jordan-Wigner transformation in order to map the problem to a free fermionic model: where {c i } are fermionic annihilation operators associated with each site. After the transformation we get the free Hamiltonian: Since this representation is quadratic in the creation and annihilation operators, the Hamiltonian is block diagonal in the momentum eigenbasis. For this reason, it is useful to decompose {c j } in Fourier modes as: Substituting in the Hamiltonian one gets: where we used the shorthand notations E k (h) := J(2h − 2 cos(k)) and ∆ k := J(2 sin(k)). Choosing the basis |k, −k , ordered as { |1, 1 , |0, 0 , |1, 0 , |0, 1 }, we can rewrite each H k as: Then, thanks to the block diagonal form of the Hamiltonian, the exponential state factorises in the tensor product: where the index k runs over positive momenta only, so to account for the choice of the basis |k, −k . This rewriting also implies that it is sufficient to study the properties of the 4x4 matrix (G8) to understand the physics of the system. In particular, we can obtain the partition function as: Tr[e −βH k (h) ].
In the limit L 1, the discrete set of momenta becomes approximately continuous in the Brillouin zone. In this regime, one can rewrite the logarithm of the partition function as: We can now pass to evaluate the y-covariance. First, notice that the variation of the Hamiltonian also factorises as: where ∂ h H k (h) can be written in matrix form, in analogy with Eq. (G8), as: Thanks to the factorisation of the exponential in Eq. (G9) and the additive structure ofḢ(h), we can rewrite: cov y h (Ḣ,Ḣ) = Tr π(h) 1−y (Ḣ − Ḣ π(h) ) π(h) y (Ḣ − Ḣ π(h) ) = where the function C(k, y, h) is given by: The cumulant generating function of the dissipated work is then obtained by plugging this expression into the definition Eq. (24), which gives: In Fig. 3 we presented the dependence of cov y h on y for different temperatures. As it was noticed from Eq. (27), if the y-covariance would be constant when varying y, we would regain a Gaussian distribution. This is the case at high temperature. On the other hand, one can see that for higher and higher values of β the non-Gaussian effects become more evident, and the y-covariance starts being more and more sensitive to variations of y. Moreover, it should be noticed that this effect is more prominent for low values of h, a signal reminiscent of the zero temperature phase transition to a magnetic long range order.
In order to take the thermodynamic limit we consider the CGF of the variable −β(w−∆F ) L , which we will denote by K diss resc (λ). This rescaling makes the average dissipation finite for L → ∞. Moreover, from the definition of the CGF, it is straightforward to verify the general property: K cX (λ) = K X (cλ). Then, we can rewrite Eq. (G17) as: where in the last line we expanded in powers of x/L for L 1. Finally, notice that Eq. (G18) can be cast in a more compact form as: • in most part of the paper the Hamiltonians will carry a discrete or continuous index; we will use the letter i in the first case, and the letter t for the latter one. More precisely, we will often use a continuous description by approximating the discrete path H 0 → H 1 → ... → H N by a continuous oneH t with t ∈ (0, 1); so that H i/N = H i . Then for example, we have thatḢ i/N = lim N →∞ N (H i+1 − H i ). We will abuse notation and write H t instead ofH t so that the index t indicates a continuous description (whereas the subindex i indicates a discrete one).
• whenever an object is function of an indexed Hamiltonian alone, the index passes to the state. For example, Z i stands for Z(H i ), or π i for π(H i ); • we will denote by γ the path in the parameters space defining the protocol; • we will use the notation Var