Contributions from populations and coherences in non-equilibrium entropy production

The entropy produced when a quantum system is driven away from equilibrium can be decomposed in two parts, one related with populations and the other with quantum coherences. The latter is usually based on the so-called relative entropy of coherence, a widely used quantifier in quantum resource theories. In this paper we argue that, despite satisfying fluctuation theorems and having a clear resource-theoretic interpretation, this splitting has shortcomings. First, it predicts that at low temperatures the entropy production will always be dominated by the classical term, irrespective of the quantum nature of the process. Second, for infinitesimal quenches, the radius of convergence diverges exponentially as the temperature decreases, rendering the functions non-analytic. Motivated by this, we provide here a complementary approach, where the entropy production is split in a way such that the contributions from populations and coherences are written in terms of a thermal state of a specially dephased Hamiltonian. The physical interpretation of our proposal is discussed in detail. We also contrast the two approaches by studying work protocols in a transverse field Ising chain, and a macrospin of varying dimension.

When a system relaxes to equilibrium, in contact with a heat bath, quantum coherences are known to contribute an additional term to the entropy production [21,36,37], which quantifies the amount of irreversibility in the process. A similar effect also happens in unitary work protocols [38,39]. To be concrete, we focus on the latter and consider a scenario where a system is described by a Hamiltonian H t = H(g t ), depending on a controllable parameter g t . The system is initially prepared in thermal equilibrium at a temperature T , such that its initial state is the thermal state ρ th 0 ≡ ρ th (g 0 ) = e −βH 0 /Z 0 , where β = 1/T and Z 0 = tre −βH 0 is the partition function. At t = 0, a work protocol g t , that lasts for a total time τ, is applied to the system, driving it out of equilibrium [40,41]. Letting U denote the unitary generated by the drive, the state of the system after a time τ will be In general, ρ τ will be very different from the corresponding equilibrium state ρ th τ = e −βH τ /Z τ . This difference is captured by the entropy production (also called non-equilibrium lag in this context) [42][43][44][45], where S (ρ||σ) = trρ(ln ρ − ln σ) 0 is the quantum relative entropy. The non-equilibrium lag is directly proportional to the irreversible work [46][47][48], Σ = β W − ∆F , where W = tr H τ ρ τ − H 0 ρ th 0 is the work performed in the process and ∆F = F(g τ ) − F(g 0 ) is the change in equilibrium free energy, F(g) = tr H(g)ρ th (g) − T S (ρ th (g)) (with S (ρ) = −tr(ρ ln ρ) being the von Neumann entropy). Due to its clear thermodynamic interpretation, Σ has been widely used as a quantifier of irreversibility, both theoretically [45][46][47][48][49][50][51][52][53] and experimentally [54][55][56][57][58][59][60][61][62][63][64][65].
The entropy production Σ in Eq. (2) contains contributions of both a classical and quantum nature. This is linked with the fact that the work protocol g t can modify the Hamiltonian H(g) in two ways. On the one hand, it may alter the spacing of the energy levels; and, on the other, it may rotate the eigenvectors ( Fig. 1(a)). The latter is directly associated with quantum coherence and to the fact that [H(g t 1 ), H(g t 2 )] 0, for two different times t 1 , t 2 . It therefore has no classical counterpart, and corresponds to a fundamental feature distinguishing classical and quantum processes. In general, these two processes will become mixed, and hence identifying how each physical process contributes to Σ is in general a challenging task. In the literature, a popular choice is the splitting put forward in [20,21,36,38]: where arXiv:2102.11244v2 [quant-ph] 29 Apr 2021 Eq. (10) (c) (d) FIG. 1. Difficulties in identifying the classical and quantum contributions to the entropy production. (a) A work protocol H(g t ) can modify the Hamiltonian in two ways: altering the energy level spacings, which can be viewed as a semi-classical effect, and rotating the energy eigenbasis, which is a fully quantum property with no classical counterpart. (b) The entropy production (2) compares the final state ρ τ with the reference thermal state ρ th τ . To understand its classical and quantum contributions, the splitting (3) uses an intermediate state D Hτ (ρ τ ). Conversely, in this paper we introduce the splitting (10), which uses the stateρ th τ , in Eq. (17). (c) The contributions Γ qu and Γ cl of Eq. (3) in a minimal qubit model, as a function of βω, plotted using Eqs. (6) and (7) with θ = 1.1 (see text for details). (d) The new splitting (10) (see also Eq. (11)), which yields physically more reasonable results at low temperatures.
with D H (ρ) being the super-operator that completely dephases the state ρ in the eigenbasis of H (explicitly defined below, in Eq. (12)). The first term, Γ cl , measures the entropic distance between the populations of the actual final state ρ τ and those of the reference thermal state ρ th τ , and is generally identified with the classical contribution. The term Γ qu , in turn, is known as the relative entropy of coherence and compares the final state ρ τ with the dephased state D H τ (ρ τ ). It hence captures the contribution from coherences in the energy basis. By construction, Γ cl and Γ qu in Eq. (3) are both non-negative, which shows that coherences increase the entropy production in the process, as compared to a fully classical (incoherent) scenario. One should also clarify that, since the changes in populations and coherences are inevitably mixed, the terminology "classical" vs. "quantum" is not entirely precise, nor is there a one-to-one relationship between this and the terms "populations" and "coherences". For instance, while Γ qu depends only on the basis rotation (coherences), Γ cl depends on both the changes in energy eigenvalues, as well as the eigenbasis rotation. Notwithstanding, as we will show, in the case of infinitesimal quenches, these distinctions can be made precise.
The splitting (3), first analyzed in [20], has been studied in the context of the resource theory of thermodynamics [21], relaxation towards equilibrium [36,37],thermodynamics of quantum optical systems [66] and work protocols in the absence of a bath [9,38,39]. At the stochastic level, both Γ qu and Γ cl satisfy individual fluctuation theorems [38], which is a very desirable property. Moreover, Γ cl has a resourcetheoretic interpretation within the resource theory of athermality [67,68], while Γ qu is a natural monotone in the resource theory of coherence [69,70]. These facts make the splitting (3) a valuable tool in understanding the relative contribution of classical and quantum features to non-equilibrium processes. However, working with various models, we have observed that this splitting behaves strangely, even in some simple protocols. More specifically, we identify two main shortcomings.
The first concerns the relative magnitudes of Γ qu and Γ cl : At low temperatures, Γ cl will always be much larger than Γ qu . The reason is purely mathematical: Γ qu is a special kind of relative entropy because it can be expressed as a difference between two von Neumann entropies, as in the second equality of (5). As β → ∞, ρ th τ tends to a pure state and hence S (ρ τ ) tends to zero, while S (D H τ (ρ τ )) ∈ [0, ln d], where d is the dimension of the Hilbert space. As a consequence, Γ qu will always remain finite. The term Γ cl , on the other hand, generally diverges when the support of ρ τ is not contained in that of ρ th τ [71], meaning Γ cl will grow unbounded when β → ∞. This implies that it is impossible to construct a lowtemperature process where the quantum term dominates.
More precisely, consider again the two types of drivings depicted in Fig. 1 (a): one that alters the spacing of the energy levels (associated here to a classical process), and one which may rotate eigenvectors (associated here to a quantum process). At strictly zero temperature, a Gibbs state is invariant under the first class of protocols, and hence we may expect that any entropy production in (2) has a quantum origin. However, the opposite identification arises in the splitting (3). The reason for this apparent contradiction is rather simple: The splitting (3) is not characterising whether the driving generates quantum coherence or not; rather, given a possibly coherent process, it characterises how much the final diagonal and off-diagonal terms contribute to the total entropy production.
This issue can be neatly illustrated by a minimal qubit model. Consider a qubit which starts at H 0 = ωσ z and is suddenly quenched (U = 1) to H τ = ω(σ z cos θ + σ x sin θ) (where σ α are Pauli matrices). In this quench the energy levels remain intact and all that happens is that the eigenbasis is rotated by an angle θ. This is thus, by all accounts, a highly quantum process. The entropy production (2) for this model reads where t = tanh(βω) ∈ [0, 1]. On the other hand, the coherent contribution Γ qu in Eq. (5), reads − 1 2 ln 1 + sinh 2 (βω) sin 2 θ .
A plot of Γ qu and Γ cl = Σ − Γ qu is shown in Fig. 1(c) as a function of βω, for θ = 1.1. As can be seen, in general both quantities are comparable in magnitude. But, as the temperature goes down (β goes up), the classical contribution becomes increasingly larger and eventually dominates. Thus, at very low temperatures, most of Σ comes from the population term Γ cl and very little from coherences. The above considerations highlight the fact that splitting the total entropy production (2) in a classical and quantum contribution may be highly non-trivial, and that different splittings might provide different insights. In particular, we argue that the splitting in Eq. (3) does not appropriately distinguish coherent from non-coherent drivings (see Fig. 1), but instead characterises how populations and off-diagonal terms contribute to entropy production. In this work, we will propose a new complementary splitting that better incorporates the difference between coherent and non-coherent drivings.
A second issue with the splitting (3) concerns infinitesimal quenches. This is a very important scenario, widely studied in the context of critical systems [72][73][74][75] and quasi-isothermal processes [29,30]. The idea is to analyze the entropy production perturbatively, for a small instantaneous quench of the work parameter, from g to g + δg. The problem with Γ qu and Γ cl in this case is that, as will be shown, the parameter δg appears multiplied by a factor that increases exponentially with β. Hence, the radius of convergence of Γ qu and Γ cl , in δg, tends to zero exponentially fast as β → ∞. For Σ, no such issue arises. This is again well illustrated by the qubit example in Eqs. (6) and (7), where the quench parameter is now the angle θ. We see that Σ in (6) can be readily expanded in powers of θ, for any temperature β (or any t = tanh(βω)). The same is not true for Γ qu , however. The problem is in the third term of Eq. (7), which is a function of x = sinh 2 (βω) sin 2 θ. This quantity appears inside a logarithm, in the form ln(1+ x). However, a series expansion of ln(1 + x) only converges if |x| < 1. And since the prefactor sinh 2 (βω) grows exponentially with β, at low temperatures, extremely small values of θ are required to validate a series expansion.
More generally, one can readily show that for Σ this issue does not arise. If we use Σ = β W − ∆F , we find in the case of infinitesimal quenches that where ∆H = H(g 0 + δg) − H(g 0 ) and ∆F = F(g 0 + δg) − F(g 0 ).
A series expansion of Σ in δg therefore amounts to two things. First, an expansion of ∆H in powers of δg, which is entirely independent of β. And second, an expansion of F(g), which is an analytic and generally smooth function (except possibly at a critical point [73]). Indeed, if H(g) is linear in g, the leading order contribution to the expansion becomes [74] Σ − 1 2 showing that Σ is simply proportional to the equilibrium susceptibility, a textbook quantity used throughout equilibrium statistical mechanics.
The above results show that, despite its interesting properties (individual fluctuation theorems and resource-theoretic interpretation), the splitting (3) is not a fully satisfying splitting of the entropy production into a classical and quantum contribution (in the sense described in Fig. 1 (a)). In order to capture the difference between coherent and non-coherent drivings, in this paper we propose a different splitting, which is inspired by the recent results of [30]. We label it as The actual definitions of Λ qu and Λ cl will be given below in Sec. II and a stochastic trajectories formulation will be given in Sec. III. A comparison in the case of the minimal qubit example is also presented in Fig. 1(d). In this case, using the results of Sec. II, one finds the following elegant expression for Λ qu (to be contrasted with Eq. (7)): As seen in Fig. 1(d), Λ qu and Λ cl behave as desired: Since the process is highly coherent, Λ cl is very small; and as the temperature goes down, Λ qu grows monotonically, showing that cold processes have higher contributions from the coherences. The features discussed in Fig. 1 are not restricted to quenches. To illustrate that we show in Fig. 2 another qubit example, where the process is assumed to be cyclic, with H τ = H 0 = ωσ z , and the unitary is taken to be generated by an x-pulse with a duration τ; that is, U = e −iτ(H 0 +h x σ x ) , where h x is the pulse intensity. Fig. 2 illustrates the results for ω = 1, h x = 1.3 and two choices of τ: in the upper panels τ = 0.4 and in the lower panels τ = 1. The results show that for (3) the behavior is always roughly the same, with Γ cl always eventually dominating at low temperatures. Conversely, for the new splitting (10) a richer competition is observed. Depending on the parameters we may either have Λ qu dominating, or Λ cl , or both.
As we will show in this paper, our new splitting (10) more accurately distinguishes which part of the entropy production is generated by a commuting or non-commuting drive. This provides a complementary approach to the standard splitting in Eq. (3), which instead describes how populations and coherences in the final state contribute to entropy production. On the other hand, we also note that Λ qu and Λ cl do not share some of the nice properties of Γ qu and Γ cl . First, Λ qu cannot be directly linked with a monotone for coherence or asymmetry [70]. Second, while Λ cl always satisfies an individual fluctuation theorem, Λ qu only does so in the case of infinitesimal quenches. Different properties of each splitting are highlighted in Table I. We also show that for infinitesimal Splitting of the entropy production in a cyclic qubit model, The curves were computed using Eqs. (14), (15), (16), (19) and (18)  To illustrate the usefulness of our results, we analyze our new splitting in two quantum many-body problems. Previous works have focused on the behaviour of the statistics of work and entropy production Σ for quantum quenches [72][73][74][75][76][77], with emphasis in quantum phase transitions [73,[78][79][80][81][82][83][84][85].
Motivated by this, we analyze in Sec. IV the transverse field Ising model (TFIM), and discuss the behavior of (10) in the vicinity of the quantum critical point. This is complementary to the analysis put forth in [39], which studied Eq. (3). Then, in Sec. V, we consider a macrospin of varying size and focus on the full statistics of Λ qu and Λ cl , including their probability distributions and their first four cumulants. We finish with conclusions and future perspectives in Sec. VI.

II. SPLITTINGS OF THE ENTROPY PRODUCTION
In this section we introduce our alternative splitting of the entropy production [Eq. (10)]. We focus for now at the level of averages; the corresponding stochastic formulation will be presented in Sec. III.
Let O denote any Hermitian observable and decompose it as O = α o α Π α , where Π α are projectors onto the subspaces with eigenvalues o α . We define the dephasing operation The rationale of the splitting Eq. (3) was to introduce an intermediate step, associated with the state D H τ (ρ τ ) ( Fig. 1(b)). This represents the final state ρ τ dephased in the eigenbasis of the final Hamiltonian. If the process generates coherences, this state will differ from the actual final state ρ τ and their entropic distance will be precisely Γ qu in Eq. (5). For convenience, we introduce the non-equilibrium free energy, associated with the final Hamiltonian H τ Non-equilibrium free energies depend on two parameters, H and ρ. However, in this paper, we will henceforth only need free energies defined with respect to H τ , so we write it more simply as F(ρ). In terms of F, the entropy production (2) can be written as Similarly, one can also express Γ qu and Γ cl in terms of free energy differences. Since tr H τ D H τ (ρ τ ) = tr H τ ρ τ , one finds that which clearly add up to Σ. The splitting (3) uses D H τ (ρ τ ) as intermediate state.
Our new splitting (10) follows a similar logic, but in reverse: Instead of working with ρ τ dephased in the basis of H τ , we work with H τ dephased in the basis of ρ τ . More precisely, we definẽ which is a thermal state based only on the incoherent part of H τ , in the basis of ρ τ (as a consequence, [ρ th τ , ρ τ ] = 0). With this in mind, we now define which add up to Σ, as in Eq. (10). The first term, Λ cl , compares the two commuting states ρ τ andρ th τ and is hence associated with their population mismatch. The nonnegativity of Λ cl becomes evident by noting that it can also be written as The term Λ qu , on the other hand, compares ρ th τ ∝ e −βH τ with ρ th τ ∝ e −βD ρτ (H τ ) . Unlike Λ cl , the contribution Λ qu cannot be written as a relative entropy. In fact, written down explicitly, it reads Notwithstanding, as shown in Appendix A, it turns out that Λ qu is still non-negative, and zero if and only if [H τ , ρ τ ] = 0. Throughout this paper we will provide several additional justifications as to why the choices (18) and (19) are physically reasonable, starting in Sec. II A. But before doing so, let us briefly revisit the minimal qubit model defined above Eq. (6). The process is a quench (U = 1), so ρ τ = ρ th 0 . Hence, all we need to do in order to compute Λ qu is to dephase the final Hamiltonian H τ = ω(σ z cos θ + σ x sin θ) in the basis of ρ th 0 . Or, what is equivalent, in the basis of H 0 . The result is thus simply D ρ τ (H τ ) = ω cos(θ)σ z . Using this in (19) yields Eq. (11), which is the result plotted in Fig. 1(d) and discussed in Sec. I.

A. Infinitesimal quenches
The physics of the problem becomes particularly simpler in the case of infinitesimal quenches. We therefore now specialize the above results to this scenario. This will provide strong justifications in favor of the new splitting (10). Furthermore, in this limit the splitting (10) becomes equivalent to the one recently put forward in [30]. More precisely, in [30] the authors describe quasi-isothermal processes as a series of infinitesimal quenches, and in particular consider how Σ splits into a classical and quantum contribution. Focusing on a single infinitesimal quench, both approaches become directly comparable and, as we will show, agree with each other.
We thus analyze what happens if we take U = 1, and assume that H changes only by a small amount ∆H (i.e., we write H τ = H 0 + ∆H). Since U = 1, the state of the system remains unchanged: ρ τ = ρ th 0 . Therefore, dephasing H τ in the basis of ρ τ is equivalent to dephasing in the basis of H 0 : Let us define the dephased (incoherent) and coherent parts of the perturbation ∆H, in the initial energy basis, ∆H d = D H 0 (∆H) and ∆H c = H τ − D H 0 (H τ ). Then, following a procedure detailed in Appendix B of Ref. [30], one may show that,ρ th where . . . 0 = tr{. . . ρ th 0 } and J ρ is a super-operator defined as We see that both ρ th τ andρ th τ can be expanded essentially in a power series in β∆H. Conversely, the same is not true for the state D H τ (ρ th 0 ) entering (16) and (15). In fact, one may show that to order ∆H [86] Even though this is an expansion in ∆H, the dependence on β enters in a highly non-trivial way. This explains the nonanalytic behavior of Γ cl and Γ qu at low temperatures, discussed in Sec. I.
Plugging (23)- (24) in Eqs. (2), (20) and (21) leads, up to second order, to where we used the fact that ∆H d 0 = ∆H 0 . The interesting aspect of these results is that, within this infinitesimal quench limit, Λ cl and Λ qu are found to be related to Σ via the simple separation of the perturbation, ∆H = ∆H d + ∆H c , into a dephased and a coherent part. These results also coincide with the splitting proposed in [30].
An additional justification for the splitting (10) can be given in terms of the fluctuation-dissipation relation (FDR). As shown in Refs. [29,30], Eq. (27) can also be written as where Var 0 [∆H] = ∆H 2 0 − ∆H 2 0 , is the variance of the perturbation, and is a measure of quantum coherence, associated with the socalled Wigner-Yanase-Dyson skew information [87] For incoherent processes one recovers the usual FDR Σ = β 2 2 Var 0 [∆H] [46]. But when the process is coherent, the FDR is broken by a term −βQ. Repeating the same procedure for Λ cl and Λ qu , one readily finds that Whence, Λ cl always satisfies a standard FDR, and all violations are associated to Λ qu . This provides additional justification as to why Λ qu is referred to as a quantum contribution.
In the case of high temperatures (β → 0), one may show that Q ∼ O(β) 3 . Moreover, the state entering the variances in Eq. (33) can be replaced with the maximally mixed state I/d. As a consequence, we find that to leading order in β, Both contributions are thus found to scale as β 2 in this limit, which agrees with the observations in Fig. 1(c) and (d). However, their relative contribution will be determined by the variance of ∆H d and ∆H c in the maximally mixed state; hence, which term will be dominant will depend on the details of the process (either a commuting or a non-commuting drive). This is also expected to remain true for general drives.

III. STOCHASTIC TRAJECTORIES
We now discuss how to formulate the splittings (3) and (10) at the level of stochastic trajectories, based on a standard two-point measurement (TPM) scheme [48]. Since Σ = β( W − ∆F), the statistics of Σ can be obtained solely from measurements in the eigenbasis of the initial and final Hamiltonians. As first shown in [38], a major advantage of the original splitting (3) is that this remains true when assessing the individual contributions Γ cl and Γ qu ; that is, no additional measurements are necessary. As we will now show, the same is also true for Λ cl and Λ qu [Eq. (10)]. This means that both splittings can be assessed, at the stochastic level, with the same amount of information as a standard TPM.
Irrespective of the splitting one is interested in, the protocol may therefore be described as follows. Initially the system is in the thermal state ρ th 0 , associated with the Hamiltonian H 0 = i 0 i |i 0 i 0 |. The first measurement is performed in the basis |i 0 , which occurs with probability p 0 i = e −β 0 i /Z 0 . Conversely, the second measurement is performed at time τ, after the map (1), and in the eigenbasis of the final Hamiltonian H τ = j τ j | j τ j τ |. The bases {|i 0 } and {| j τ } are, in general, not compatible.
The conditional probability of finding the system in | j τ given that it was initially in |i 0 is | j τ |U|i 0 | 2 . The probability associated with the forward protocol |i 0 → | j τ is thus The dynamics is defined as being incoherent when | j τ |U|i 0 | 2 = δ i, j , which means U is not able to generate transitions between states of the initial and final Hamiltonians. Similarly, in the backward protocol the system starts in ρ th τ and one measures first in the basis of H τ , yielding | j τ with probability p τ j = e −β τ j /Z τ . The timereversed unitary U † is then applied, after which one measures in the basis |i 0 of H 0 . This yields the backward distribution The entropy production associated to the trajectory |i 0 → | j τ is now defined as usual: The second equality follows from the fact that | i 0 |U † | j τ | 2 = | j τ |U|i 0 | 2 . As a consequence, σ[i, j] depends only on the equilibrium populations p 0 i and p τ j , associated with the initial and final Hamiltonians. As can be readily verified, σ[i, j] = i, j σ[i, j]P F [i, j] = Σ, returns precisely Eq. (2). In addition, σ[i, j] also satisfies an integral fluctuation theorem e −σ = 1 (see Eq. (43) for more details).
A. Stochastic definitions for the splittings (3) and (10) Following [38], we now define stochastic quantities associated to Γ cl and Γ qu . In order to do that, we first write the In passing, we note that q τ j = i P F [i, j], so q τ j can also be interpreted as the marginal distribution of the final measurement. As shown in [38], we may now define Clearly , which is the stochastic analog of (3). Moreover, γ cl [i, j] = Γ cl and γ qu [i, j] = Γ qu . Similarly, we construct stochastic quantities for the new quantities Λ cl and Λ qu in Eq. (10). The central object now is the thermal stateρ th τ , defined in Eq. (17) and associated with the Hamiltonian D ρ τ (H τ ). Since the system evolves unitarily, That is, ρ τ has the same populations p 0 i as ρ th 0 , but a rotated eigenbasis. Based on this, we can now write Eq. (17) as where˜ τ i = ψ i |H τ |ψ i are the eigenvalues of the dephased Hamiltonian D ρ τ (H τ ) and F(ρ th τ ) is the same free energy as that appearing in Eq. (18). We then define from which we may compute the CGF, K σ (v) = ln e −vσ .
In addition, from the CGF we may compute any cumulant of σ as with κ 1 (σ) = Σ being the mean in Eq. (2). We may also compute the joint CGF of γ cl and γ qu , defined as K γ cl ,γ qu (v, u) = ln e −vγ cl −uγ qu . With similar manipulations, it may be written as The CGF of σ = γ cl + γ qu , Eq. (43), is recovered by setting . The reduced CGFs of γ cl and γ qu are found by setting u = 0 or v = 0, respectively: From this one may verify that γ cl and γ qu individually satisfy fluctuation theorems Note also that, except in certain particular cases, Eq. (46) cannot be written as a sum of two CGFs, which means γ cl and γ qu are statistically dependent. Similarly, we compute the joint CGF of λ cl and λ qu , defined as K λ cl ,λ qu (v, u) = ln e −vλ cl −uλ qu . It reads The reduced CGFs of λ cl and λ qu are again found by setting u = 0 an v = 0, Once again, λ cl and λ qu are, in general, statistically dependent. Eq. (51) implies that λ cl satisfies a fluctuation theorem, But the same is not true for λ qu . Notwithstanding, as we will show, this property is recovered in the limit of infinitesimal quenches.

C. Infinitesimal quenches
As before, we now specialize the above expressions to the case of infinitesimal quenches. Since U = 1, the path probability reduces to P F [i, j] = | j τ |i 0 | 2 p 0 i . Moreover, since ∆H is assumed to be small, | j τ will be close to |i 0 and τ j will be close to 0 j . For concreteness, we assume that the spectra of H 0 is non-degenerate. Standard perturbation theory then yields, to order ∆H 2 , where ∆H i j = i 0 |∆H| j 0 and E (2) j = l j |∆H jl | 2 /( 0 j − 0 l ). Note that if we split ∆H = ∆H d + ∆H c , the first non-trivial contribution of the former is ∆H j j , while that of the latter is E (2) j . Similarly, the eigenstates | j τ of the final Hamiltonian can be expanded as wheref and i . Note howf j depends only on the diagonal part of the perturbation, ∆H d . This is in line with Eq. (23). Conversely, f j , which is associated with the full probabilities p τ j , also has an additional contribution from E (2) j , which is the term associated to coherences.
We are now in the position to discuss the analyticity of the entropy production and its splittings, at the stochastic level. A series expansion of ln(1 − x) is convergent only for |x| < 1. Thus, since p 0 i /p 0 j = e −β( 0 i − 0 j ) is a well behaved function, the analyticity of σ, λ cl and λ qu are all conditioned on having | f j | < 1 and |f j | < 1, which is satisfied if β|∆H i j | 1, as intuitively expected. Thus, the quantities of our new proposed splitting (10) behave, from an analytical point of view, similarly to the full entropy production.
On the other hand, Eqs. (63)-(64) rely on |s j | < 1. Each s j in (61) is a weighted contribution from all energies 0 , with j. At low temperatures, those energies for which 0 < 0 j will yield an exponentially large contribution 1 − e −β( 0 − 0 j ) to the sum. Conversely, those with 0 > 0 j will contribute negligibly. The expansion is thus not in powers of β∆H, which is also visible from (26). Instead, it is an expansion in powers of ∆H, with coefficients that depend exponentially in β. Violating the condition |s j | < 1 is thus exponentially easier at low temperatures. These results show that the shortcomings illustrated in Sec. I, are in fact absolutely general.
To better illustrate this discussion, we revisit the minimal qubit model example treated in Sec I. Initially the system's Hamiltonian is H 0 = ωσ z , and after an instantaneous quench it becomes H 1 = ω(σ z cos θ + σ x sin θ), where we consider θ to be small. The problematic term in this case is s 1 [Eq. (61)], which is given by In comparison, we havẽ f 1 = 2βω sin 2 (θ/2)(1 + tanh βω) (68) × 1 + 2βω sin 2 (θ/2) tanh βω In Fig. 3 we plot Eqs. (67)-(69) as a function of βω, for θ = 0.1. The condition for analyticity of γ cl and γ qu in this case, |s 1 | < 1, is rapidly violated with increasing βω. For σ, λ cl and λ qu , instead, | f 1 | < 1 and |f 1 | < 1 for a much larger range of temperatures. It is also interesting to note that, at low temperatures, the excited state thermal probabilities p 0 1 ,p τ 1 and p τ 1 ∝ e −βω all tend to zero exponentially as e −βω . Conversely, q τ 1 tends to (sin θ/2) 2 . This corroborates the use of thermal states, such as (17), as intermediate states for the splitting of Σ, as it ensures that the resulting functions are analytic.
We now move on to discuss what becomes of the CGFs of Sec. III B in the infinitesimal quench regime. We start we the CGFs of σ, λ cl and λ qu in Eqs. (43), (50)- (52). Using Eqs. (62), (65) and (66), together with the path probability P F [i, j] = p 0 i | j τ |i 0 | 2 , we find to order ∆H 2 , that where These results are quite illuminating. Eq. (70) implies λ cl and λ qu become statistically independent in this limit. Moreover, since K σ (v) = K λ cl ,λ qu (v, v), we now find that This means that all cumulants of σ can be split as a sum of the cumulants of λ cl and λ qu : κ n (σ) = κ n (λ cl ) + κ n (λ qu ). For all intents and purposes, the two channels of entropy production, Λ cl and Λ qu , may thus be regarded as stemming from independent processes: Λ cl gives the entropy production associated with a quench from H 0 → D H 0 (H τ ), while Λ qu is associated with a second quench from D H 0 (H τ ) → H τ . We also note that, from Eq. (73), it can now be seen that in this limit λ qu satisfies an integral fluctuation theorem: e −λ qu = 1. In contrast, for the original splitting (3), we have (74) Once again, a series expansion of (1 − s j ) −x is convergent only if |s j | < 1. However, if |s j | < 1 is satisfied, which happens for sufficiently high temperatures, one may show that, to order ∆H 2 , we can also split K γ cl ,γ qu (v, u) = K γ cl (v) + K γ qu (u). And, what is more important, K γ cl and K γ qu coincide with K λ cl and K λ qu respectively. Whence, at sufficiently high temperatures the splittings (3) and (10) coincide, even at the stochastic level. However, this is only true for infinitesimal quenches. Otherwise, the two splittings may differ, even at high temperatures.

IV. TRANSVERSE FIELD ISING MODEL
We now turn to discuss applications of our framework. We begin with the behavior of the splitting (10) in the onedimensional transverse field Ising model (TFIM), which is a prototypical model of a quantum critical system. An analysis of (3) for the same model was recently made in [39]. Here we aim to contrast those results with our new splitting (10). We thus restrict the analysis to quench protocols, and study the problem at the level of the averages Λ cl and Λ qu (Eqs. (18) and (19)). Non-trivial unitaries and higher order statistics will be studied in Sec V, for a different model.
We begin by introducing the model and delineating the steps to compute Λ cl and Λ qu . To make the paper selfconsistent, some additional details are provided in Appendices B and C. Consider a linear chain of N spins, each described by Pauli operators σ x,y,z j and interacting via the Hamiltonian where g is the applied magnetic field and J is the coupling strength, which we henceforth set to unity (J = 1). We consider periodic boundary conditions, σ α N+1 = σ α 1 . This model presents critical points at g = ±1, where the system changes from a ferromagnetic phase, for |g| < 1, to a paramagnetic phase, for |g| > 1.
After a series of transformations (see Appendix B) this Hamiltonian can be written as where {η k } are fermionic operators and are the single-particle energies. Here, k = ±(2n + 1)π/N, with n = 0, 1, ..., N/2 − 1, denotes the system's quasi-momenta. We consider that the system initially has a transverse field g 0 and is prepared in the thermal state ρ th 0 = e −βH 0 /Z 0 . The full expression for ρ th 0 can be found in Appendix C. Due to the structure of (76), it can be decomposed as a product over individual k modes, which greatly facilitates the calculation of all thermodynamic quantities.
The system is then decoupled from the reservoir and undergoes an instantaneous quench, where the field is changed to g τ = g 0 + δg. Since the quench is instantaneous, the state of the system remains the same, but its Hamiltonian changes, from H 0 to H τ = H 0 + ∆H, where ∆H = −δg j σ z j . Full details on the computation of Λ cl and Λ qu are provided in Appendices B and C. The overall contributions of the diagonal vs. off-diagonal is described by the Bogoliubov angle cos θ k = (g 0 − cos k)/ 0 k and sin θ k = sin k/ 0 k . And the state (17), associated with the dephased final Hamiltonian, is described by the modified energies˜ τ k = 0 k + δg cos θ k . We are interested in the thermodynamic limit (N very large), where k sums can be converted to integrals and all quantities become extensive in N. In this case, we ultimately find that and Adding both contributions recovers the full entropy production Σ, which was computed in [39,73,83] and reads (80) For comparison, in Appendix D we also present the formulas for the splitting (3), which were developed in [39]. We also mention, in passing, that Eqs. (78) and (79) are not perturbative in the quench magnitude. That is, they hold for arbitrary quenches δg. The only assumption is that U = 1. For completeness, their behavior in the infinitesimal case is presented in Eqs. (C10) and (C11). Fig. 4 compares the two splittings (3) and (10) as a function of β, with fixed g 0 = 0.75 (outside criticality) and δg = 0.01. At high temperatures, one clearly sees how both splittings coincide. This corresponds to the region of parameters where Γ cl and Γ qu are analytic. But as the system is cooled, they eventually begin to differ. In particular, Γ cl tends to grow linearly with β, while Λ cl tends to zero. For the quantum components the opposite is observed: Γ qu tends to saturate while Λ qu tends to grow. Thus, at very low temperatures Λ qu becomes the dominant contribution in (10), while Γ cl becomes the dominant one in (3).
Next we turn to the behavior near criticality. In Fig. 5 we plot Γ cl , Γ qu , Λ cl and Λ qu as a function of the initial field g 0 , for different values of β (focusing on low temperatures) and fixed quench magnitude of δg = 0.01. The full entropy production Σ behaves similarly to Λ qu in Fig. 5(a); for any finite T it presents a peak at g 0 = 1, which eventually tends to a divergence as β → ∞. As is clear by comparing Figs. 5(a) and (b), the dominant contribution to the splitting (10) is Λ qu . Moreover, Λ qu is found to always grow (and eventually diverge) with β at g 0 = 1, while Λ cl sharply drops to zero. Conversely, for the splitting (3), we find in Figs. 5(c), (d) that the dominant contribution is instead that of the populations Γ cl . Crucially, we find that in this case Γ qu remains finite as β → ∞, while Γ cl diverges [39]. We also call attention to the non-monotonic dependence on β, of the quantities in Fig. 5(c). This is an artifact of the fact that Γ qu is scaled by β. The quantity Γ qu itself is monotonic, but its behavior changes from β 2 at high temperatures, to β 0 in low temperatures [39].
As highlighted in [39], the entropy production in this limit results entirely from the changes in occupations, i.e. creation/annihilation of particles, in the modes ±k, when the quench is performed. This enters in Γ cl as a population mismatch between the initial and final equilibrium states. Conversely, in the split (10), it enters in Λ qu as resulting from the rotating energy basis. On the other hand, Λ cl only quantifies the contribution to the entropy production resulting from a change in the energy levels given by˜ τ k − 0 k = δg cos θ k . In the low temperature limit, only the ground and low lying excited states are relevant, and close to the critical point g 0 = 1, the latter corresponds to creating excitations with momentum k → 0; but one can easily show that at g 0 = 1, cos θ k = | sin(k/2)|, which goes to zero when k → 0. This explains the drop in Λ cl at this point. Overall, Fig. 5 therefore shows that the critical properties of these quantities depend crucially on the type of splitting used.

V. MACROSPIN MODEL
Finally, we analyze our framework from the stochastic perspective developed in Sec. III. To emphasize the generality of our results, we also focus on non-quench scenarios (U 1). We consider a macrospin model, with d = 2S + 1 levels and spin operators S x , S y , S z [89]. We consider a scenario similar to that of Fig. 2: The initial and final Hamiltonians are taken to coincide, being given by H 0 = H τ = −h z S z . And the unitary is driven by a magnetic pulse in the x direction, so U = exp{−i(H 0 −h x S x )τ}. Since the Hamiltonian is cyclic, the eigenbases |i 0 and | j τ coincide. However, since the unitary is now non-trivial, the final state ρ τ = Uρ th 0 U † will contain coherences in the S z -basis.
A panel summarizing the results for the splitting (10) is shown in Fig. 6, where we plot the first four cumulants of λ cl (images (a)-(d)) and those of λ qu (images (f)-(i)), as a function of the Hilbert space dimension d and for different values of β. In Figs. 6(e) and (j), we also show exemplary plots of the full distributions P(λ cl ) and P(λ qu ), for fixed d = 200 and two values of β. For comparison, a similar panel, but for the quantities in (3), is shown in Fig. 7. Note also that some cumulants are scaled by either d or β, whenever a simple scaling rule could be found.
From these plots the following conclusions can be drawn. Concerning Fig. 6, all cumulants of λ cl are found to be intensive, saturating at a finite value when d → ∞. Conversely, all cumulants of λ qu are extensive, scaling proportionally to d. The cumulants of λ qu also scale with powers of β at low temperatures (Figs. 6(f)-(i)), but for higher order cumulants this scaling only becomes good at very low temperatures. For the splitting (3) the situation is reversed: now the cumulants of γ cl become extensive (and quite similar to those of λ qu ), while those of γ qu tend to saturate. The only exception is κ 1 (γ qu ), which is found to grow logarithmically with d. Notice that this dependence on the dimensionality is fundamentally different from what was found in the TFIM (Sec. IV), where all quantities were extensive in the number of particles. We also note that the statistics of λ cl (Fig. 6(e)) has significantly smaller support than that of λ qu . This is a consequence of the fact that, from its definition in Eq. (40), λ cl depends only on the initial points |i 0 , while λ qu depends on both |i 0 and | j τ .
The results in Fig. 6(e) indicate that even when d → ∞, the distribution of λ cl will never tend to a Gaussian. Conversely, for γ cl in Fig. 7(e), this is clearly the case. This is supported by a comparison of the corresponding cumulants in images (a)-(d), which are intensive for λ cl and extensive for γ cl . As for λ qu and γ qu , even though the histograms in Figs. 6(j) and 7(j) do not seem to indicate a Gaussian behavior, this is expected to eventually occur for sufficiently large d. For λ qu , the scal-ing with d is similar to γ cl , and hence the same argument as above applies. Conversely, the situation for γ qu is more delicate, since the first cumulant scales only logarithmically (and hence very slowly) with d. Extremely large sizes may thus be necessary for a Gaussian behavior to be observed.
One might expect that in the limit d → ∞ one should recover a classical spin model. This does not happen, however, as is evidenced by the fact that the coherent terms Γ qu and Λ qu do not vanish in this limit, but actually increase with d. The explanation for this rests essentially on a coarse-graining argument. Even though we take d → ∞, we continue to assume we have full access to all eigenstates of the system, as appears, for instance, in the dephasing operations involved in constructing the intermediate states.

VI. CONCLUSION
In this article, we studied how entropy production can be divided into a classical and quantum contribution, when a system is driven out of equilibrium. A popular choice in the literature is given in Eq. (3), see in particular [38]. This splitting has several interesting properties, including individual fluctuation theorems for each term [38] and a resource-theoretic interpretation [20,21,36]. However, we here noted it also has two major shortcomings. First, we showed that the classical contribution Γ cl in Eq. (3) dominates for highly coherent processes and at low temperatures, in contrast with what might be expected. We observed this undesired behaviour in all considered systems, from a simple driven qubit to a many-body Ising model at criticality, and identified the divergence of the relative entropy in (4), at low temperatures, as the underlying cause. Second, given a perturbation δg of the Hamiltonian, the radius of convergence of Γ qu and Γ cl tends to zero exponentially fast as β → ∞, making this splitting impractical to characterise the entropy production of quenched systems at low temperatures.
In order to overcome these shortcomings, we suggested a new splitting for the entropy production given in Eq. (10), which was motivated by the developments of Ref. [30] for infinitesimal quenches. The definition is valid arbitrarily outof-equilibrium. We also provided a formulation in terms of stochastic trajectories and a physical interpretation, highlighting how it can be obtained following a similar logic to the one behind (3). Indeed, both (3) and (10) can be understood by introducing intermediate states for comparison. The different choices, however, turn out to have crucial consequences, especially for highly coherent processes. Indeed, in the lowtemperature regime the quantum term Λ qu dominates in Eq. (10), but the classical one does in (3). For high temperatures and infinitesimal quenches, both splittings coincide. A comparison between the two approaches is summarized in Table I. More generally, our considerations illustrate that it is nontrivial to identify the classical and quantum contributions in entropy production for an arbitrarily out-of-equilibrium process. In analogy with the definition of work for coherent processes [26,90], the splitting of Σ in a classical and quantum term may not be unique, and will depend on the specific context into consideration. Nevertheless, there are some relevant scenarios where such a splitting is unambiguous. One is in a thermalization process described by either a Markovian master equation or as a resource-theoretic state transformation; in both cases, such a distinction seems to be very well captured by Eq. (3) [20,21,36]. On the other hand, when an equilibrium state is slightly moved out of equilibrium (e.g. by an infinitesimal quench), the splitting (10) provides a more accurate description of the quantum and classical contributions. In fact, in such a scenario, the entropy production can be decomposed into a classical and quantum contribution at all levels of the statistics, as shown in Sec. III B (see also [30,31,91]). For general out-of-equilibrium processes, however, classical and quantum contributions become inevitably mixed. Still, our results show that the splitting (10) has a more reasonable behaviour (i.e., the quantum term dominates at low temperatures and for highly coherent processes).
In a second part of the article, we applied these ideas to a transverse field Ising model, and to a macrospin undergoing finite time dynamics. For the Ising model, we found that the behavior close to criticality is fundamentally different for both splittings, with the quantum component playing a predominant role for (10) and the classical component being dominant in (3). For the macrospin model, we focused not only on the average, but on the full statistics, including the first four cumulants and the corresponding probability distributions. We have found that different cumulants scale with the Hilbert space dimension d in non-trivial ways, some being extensive, others intensive or even logarithmic.
We hope that these results help to motivate further investigations on the non-trivial way in which populations and coherences intermix in quantum thermodynamic processes. We are particularly interested in further understanding how this unfolds for many-body systems in general. In particular, the analysis of higher order cumulants for these models has been seldom explored in the literature, even for Σ itself. It would also be interesting to generalize the present results for open systems, undergoing generic interactions with a heat bath. This can be done for quasi-static processes, following the approach in [30]. Or it can be constructed in a controllable way using collisional models [92]. Finally, these ideas could also be extended to describe quantum correlations in bipartite systems. For instance, instead of studying the entropy production in a work protocol, one may analyze it in the context of heat exchange between two quantum correlated systems, which are locally thermal, as studied in [35]. energy eigenbasis, D ρ τ (H τ ) = D H τ (H τ ) = H τ , which leads to Λ qu (ρ τ ) = 0. Conversely, if we assume that Λ qu (ρ τ ) = 0 and β > 0, we must have F(ρ th τ ) − F(ρ th τ ) = 0. This implies that which means trH k τ − D ρ τ (H τ ) k = 0, ∀ k ∈ N. The case k = 0 is trivial and the case k = 1 follows directly from the definition of D ρ τ (H τ ). For the case k = 2, we use that Again, using the definition of D ρ τ (H τ ) one may verify that trD ρ τ (H τ )H c τ = 0. Therefore we are left with But since H τ − D ρ τ (H τ ) is also Hermitian, we must have H τ − D ρ τ (H τ ) = 0. Then, since D ρ τ (H τ ) = H τ , for k > 3, trH k τ − D k ρ τ (H τ ) = 0 follows trivially, and ρ τ must be incoherent in the eigenbasis of H τ , i.e., we must have [ρ τ , H τ ] = 0. ten as ρ th 0 |±k = n k =0,1, n −k =0,1 e 2β 0 k (1−n k −n −k ) 4 cosh 2 (β 0 k ) |n −k n k n −k n k |.
After the instantaneous quench, which changes the field to its final value g τ = g 0 + δg, we have the final Hamiltonian, where τ k = k (g τ ) and |n τ −k n τ k are the joint eigenstates of the post-quench fermionic operators ξ † k ξ k and ξ † −k ξ −k , which are related to the pre-quench operators {η k } according to [73] ξ k = cos(∆ k /2)η k + sin(∆ k /2)η † −k , where sin ∆ k = −δg sin(k)/ τ k 0 k . As discussed in Appendix B, Eq. (C4) shows us that the perturbation couples pairs of modes +k and −k. This is why it is more convenient to write all quantities as ρ th 0 = k>0 ρ th 0|±k and H 0 = k>0 H 0|±k , instead of a product/sum over the negative values of k.
We have all we need to compute Λ cl and Λ qu now. We just have to plug Eqs. (C2), (C5) and (C7) into (18) and (19). Because all states ρ th 0 ,ρ th τ and ρ th τ are separable in terms of ±k modes, Λ cl and Λ qu will be given as sums over k. Hence, we find Λ cl = (C9) Finally, in the limit of very large N, all k-sums can be converted to integrals and all quantities become extensive in N.
We note that Eqs. (78) and (79) do not assume that the quench is infinitesimal. All they assume is that U = 1. If, in particular, we are interested in infinitesimal quenches, then we may series expand these expressions in powers of δg, leading to Λ cl = Nβ 2 δg 2 π 0 dk 2π sech 2 (β 0 k ) cos 2 θ k , (C10) where it is clear the relation of Λ cl and Λ qu with the dephased and coherent parts of the perturbation in Eqs. (B11) and (B12). Furthermore, it is easy to check that they satisfy Eq. (33).