A convenient Keldysh contour for thermodynamically consistent perturbative and semiclassical expansions

The work fluctuation theorem (FT) is a symmetry connecting the moment generating functions (MGFs) of the work extracted in a given process and in its time-reversed counterpart. We show that, equivalently, the FT for work in isolated quantum systems can be expressed as an invariance property of a modified Keldysh contour. Modified contours can be used as starting points of perturbative and path integral approaches to quantum thermodynamics, as recently pointed out in the literature. After reviewing the derivation of the contour-based perturbation theory, we use the symmetry of the modified contour to show that the theory satisfies the FT at every order. Furthermore, we extend textbook diagrammatic techniques to the computation of work MGFs, showing that the contributions of the different Feynman diagrams can be added to obtain a general expression of the work statistics in terms of a sum of independent rescaled Poisson processes. In this context, the FT takes the form of a detailed balance condition linking every Feynman diagram with its time-reversed variant. In the second part, we study path integral approaches to the calculation of the MGF, and discuss how the arbitrariness in the choice of the contour impacts the final form of the path integral action. In particular, we show how using a symmetrized contour makes it possible to easily generalize the Keldysh rotation in the context of work statistics, a procedure paving the way to a semiclassical expansion of the work MGF. Furthermore, we use our results to discuss a generalization of the detailed balance conditions at the level of the quantum trajectories.


Introduction
Green's functions (GFs) are commonly used to tackle problems of quantum transport in areas like molecular transport, electronics, superconductivity, nanojunctions and nuclear physics [1][2][3][4][5].They can also be used to quantify the linear response to external perturbations, a fundamental tool to analyse experiments, for instance in magnetic resonance, Raman spectroscopy and crystallography [6][7][8].Path integral techniques on the other hand, are convenient to study the relation between the quantum and the classical (stochastic) dynamics in open quantum systems [9][10][11][12][13].They are also very suitable for studying the physics of phase transitions, instantons and critical phenomena in general [14,15].In this seemingly heterogeneous landscape, the theoretical framework introduced by Keldysh [16][17][18] represents a unified way to extend both, the usual Feynman approach to path integration and the equilibrium theory of GFs, to nonequilibrium systems [19,20] and has recently found applications in quantum and stochastic thermodynamics [21][22][23][24].The idea behind the standard Schwinger-Keldysh technique is to simplify calculations by introducing an extended version of the time domain, called the Schwinger-Keldysh contour, that acts as a new ordered domain for the time variables of the theory.Different extensions of this formalism have proved to be suitable for computing charge and energy statistics using non-equilibrium GFs [23,[25][26][27][28][29][30][31][32][33].
In the present paper, inspired by the extension of the contour idea to quantum thermodynamics in the papers of Funo and Quan [22] and Fei and Quan [24], we focus on the issue of thermodynamic consistency of perturbative and path integral approaches based on modified contours.Our first step is to introduce a modified Keldysh contour that can be used for calculating the work MGF in the evolution of an isolated driven quantum system.We then make use of the symmetries of the theory, in particular the FT [34][35][36][37][38][39], proving that it can be seen as a geometrical symmetry of the modified contour.This is crucial to the derivation of our other results.First, using the modified contour as a starting point for a perturbative expansion of the MGF, we show that its symmetry can be used to prove that the expansion is thermodynamically consistent, that is, it satisfies the FT at every perturbative order.
Then, following textbook diagrammatic techniques, we introduce an approach to the perturbative expansion of the MGF based on Feynman diagrams.The structure of the diagrammatic theory is similar to the one of standard perturbation theory [20,40], but while the architecture is the same, and we can picture it in terms of the topology of the diagrams, the building blocks, i.e., the propagators and the contour that acts as a domain of integration for the vertex variables, are different.As an application of this technique, we compute the work statistics when a small non-linear perturbation of a quadratic Hamiltonian is switched on between two energy measurements performed at two different times, showing that the probability distribution can be expressed as a linear combination of Poisson processes.We then analyze the thermodynamic consistency at the level of single Feynman diagrams and prove that the FT can also be interpreted as a detailed balance condition relating each diagram to its time-reversed counterpart.
We proceed in the second part of the paper by generalizing the Feynman path integral technique [41] to the modified contour.We use a similar strategy as in Ref. [22], but our approach is more suitable for defining a generalization of the Keldysh rotation (used to obtain a semiclassical expansion of the path integral formulation of dissipative systems [19,[42][43][44][45]) in the context of quantum thermodynamics.We do so by introducing a symmetrization of the modified contour.After performing the symmetrization, we obtain a semiclassical expansion of the MGF where we explicitly compute the zeroth (classical) order and the first quantum correction to work fluctuations.Another benefit of the symmetrized contour approach is that it enables us to express the fluctuating work as a function of the endpoints of the quantum trajectories.This can be used to discuss the concept of detailed balance at the quantum trajectory level.
This paper is organized as follows.In Sec. 2 we introduce the Keldysh contour and its extension to the modified contour to express the MGF.In Sec. 3 we investigate the work statistics and show how to obtain the FT as a symmetry of the modified contour.In Sec. 4 we introduce the perturbation theory in terms of Green's functions and use the perturbation expansion to obtain the MGF of work in weakly perturbed quantum systems.Finally, in Sec. 5 we apply path integral techniques on the modified contour to compute the MGF.By symmetrizing the contour, we derive the MGF of work, compare it with similar approaches in the literature, and discuss the semiclassical limit in the path integral expression of the fluctuating work.

The contour in a general two-point measurement
The Schwinger-Keldysh contour makes it possible to simplify and express in an elegant way the theory of nonequilibrium many-body quantum systems.The same idea can be applied to the context of counting statistics in a double measurement scheme.In this section we briefly discuss the idea behind the contour and its generalization to counting statistics, before focusing our attention on more specific problems, like the work MGF.
Many physical properties of open quantum systems are encoded in correlation functions between pairs of operators [46].Here we will use them as a starting point to understand the intuition behind the Keldysh contour (taking a similar route as more complete textbook derivations, e.g., Ref. [20]).Given two observables A 1 , A 2 with time arguments in the Heisenberg picture given by t f , t ′ with t f ≥ t ′ we can introduce the correlation function The quantity above can be written in the Schrödinger picture after writing out explicitly the evolution operators contained in A (H) where with T and T representing the time-ordering and anti-time-ordering operator, respectively.If we expand the evolution operators in Eq. ( 2) we have to deal with three ordered products and the calculations can become cumbersome.It is however possible to simplify the expansion by introducing some "bookkeeping" time arguments A 1 → A 1 (t f ), A 2 → A 2 (t ′ ), and by letting the time-ordering operator act on A 1 (t f ), A 2 (t ′ ).In this way, the last term of Eq. (2) reduces to It is natural to ask if we can further simplify Eq. ( 2) by including also U † (t f , 0) in the time-ordered product in Eq. (4).A possible issue arises from the fact that the time arguments in U † (t f , 0) are the same as U(t f , t ′ ) and U(t ′ , 0), however, the position of these operators in the r.h.s. of Eq. ( 2) are evidently different.An elegant way to circumvent this problem is to introduce a new ordering operator T K , which orders the operators according to the position of their time arguments on a new domain γ K , called the Schwinger-Keldysh contour [19,20,40,47].This contour comprises two branches: 1.The forward branch γ − , that goes from the initial time 0 to a final time t f .The ordering operator T K behaves like the usual time-ordering on this branch, so that, for instance, we can rewrite Eq. (4) as 2. The backward branch γ + , that comes after the forward branch and goes back from the final time to 0 (see Fig. 1A).Since this branch of the contour goes backward in time, T K is naturally sorting the operators from left to right with their time arguments increasing.This means that an anti-time-ordered product can be written as an ordered product over γ + .For instance we have If we put the arguments of T K in Eqs. ( 5) and ( 6) under a single ordering, the operators with time argument on γ + will be placed to the left of those with time argument on γ − (since γ + follows γ − on the contour).In this way, U † (t f , 0) in Eq. ( 2) will be correctly placed to the left of A 1 and we finally obtain Note that the Hamiltonian H(t) is now a function with argument in the contour t ∈ γ K that is equal in each branch to the physical Hamiltonian H(t Using an extension of the Keldysh idea, we can write the MGF in a double measurement process in a simple and compact way.In this framework a generic observable Λ is measured at time 0 and time t f while the system evolves under the action of a time-dependent Hamiltonian between the two measurements.The generating function for the statistics of the difference of the outcomes of the first and second measurement is given by [37] which is a special case of Eq. ( 1) with A 1 = e λΛ(t f ) and A 2 = e −λΛ(0) , as well as t ′ = 0. Note that Eq. ( 8) assumes that the initial measurement operator Λ(0) commutes with the initial state density matrix ρ 0 , which is something we will always consider to be true.After expressing the MGF (8) in the form (7) and introducing the dummy integrations e λΛ(t dτΛ(0) , we obtain It is easy to see that the sum of the three integrals can be expressed as a single integration along a different contour γ C that consists in a modification of the Schwinger-Keldysh contour, in which two vertical branches are added at time 0 and t f (see Fig. 1B).Accordingly, the ordering operator also has to be extended to this new contour.This allows us to write where now the Hamiltonian H C is defined as where we labelled the vertical tracks of the contour drawn in Fig. 1B by γ ↑,↓ .As a last step, we point out that the dependence on the initial state ρ 0 can also be cast as an integration over an additional track of the contour.Indeed, if we define an effective Hamiltonian H M by expressing ρ 0 as log dzH M /Z(0) and express Eq. ( 10) as in Ref. [24] where γ is the contour defined in Fig. 1C, built by adding an additional track γ M along the imaginary axis, T γ denotes the time-ordering operator along this contour, and the extended Hamiltonian H γ (z) is defined as Moreover, we have defined If the initial state is a Gibbs state, we have H M = H(0) and β = 1/(k B T ) assumes the role of the inverse of the physical temperature T .Note that a similar modified contour can be derived for computing the characteristic function M (iλ, t f ), but in this case the vertical branches of length λ in Fig. 1C are replaced by horizontal ones [24].We will discuss this point in more detail in the following sections (see Fig. 5).

Work fluctuation theorem as a symmetry of the contour
If the measurement operator Λ(t) is the system Hamiltonian itself, the difference between the outcomes corresponds to the difference of the final and initial energies, i.e., the work performed on the system by an external source.The MGF of the work statistics is obtained from Eq. ( 8) by choosing Λ(0) = H(0) and Let us suppose now that the system is initialized in the Gibbs state, so that H M = H(0).In this setting it is possible to prove a fluctuation theorem (FT), a symmetry connecting the generating function M W and the generating function M rev W associated with the time-reversed work extraction process [35,48].We will now show that the FT can be seen as a symmetry of the contour γ.
We first present the standard derivation of the fluctuation theorem.To compute M rev W , we first introduce the time-reversed process, in which the system is initialized in a Gibbs state corresponding to the final Hamiltonian, i.e., ρ rev (0) = e −β H(t f ) /Z(t f ), and then undergoes a time-reversed evolution.The latter is a combination of a time-reversal operation Ξ, which is an anti-unitary operator satisfying Ξi = −iΞ, and a forward evolution with the time-reversed driving protocol, i.e., H rev (s) = H(t f − s), where 0 dsH rev (s) .Note that we suppose that the Hamiltonian does not depend on a magnetic field B, as such a case would lead to an additional minus sign, H rev (s, B) = H(t f − s, −B).Using Eq. ( 15) and ΞU rev (t f , 0)Ξ = U † (t f , 0), we can deduce the following identity where we introduced the shorthand notation U ≡ U(t f , 0).After introducing the free energy difference as Z(t f )/Z(0) = e −β∆F , the above result leads to the fluctuation symmetry To formulate this symmetry in the contour framework, we first note that the entire dependence of Eq. ( 12) on the counting field and the final time is established through the modified contour γ in Fig. 1C.With this in mind, every step of the derivation ( 16) can be seen as a transformation of the contour itself as shown in Fig. 2. By making the λ dependence explicit, this geometrical symmetry corresponds to γ λ = γ rev −λ−β , where γ rev is the contour associated to the time-reversed generating function (panel D in Fig. 2).The weighted ratio between the generating functions of the forward and time-reversed processes becomes which is equivalent to Eq. ( 17).The contour-based proof of the FT is a structural symmetry of any theory based on the contour γ of Fig. 1C.Note that the idea of using a modification of the contour to encode symmetries that are intrinsic in the quantum theory has already proved to be a valid tool in non-equilibrium physics.For instance, outside the double energy measurement framework, it has been used to show that equilibrium theories share a fundamental symmetry at the level of the Schwinger-Keldysh action [49][50][51][52].

Thermodynamically consistent perturbation theory
The modified contour γ can be used to build perturbative expansions of work generating functions in weakly perturbed quantum systems.Following the standard approach of perturbation theory, we decompose the contour Hamiltonian as H γ (z) = H 0 (z) + χ(z)H 1 (z) where H 0 (z) and H 1 (z) are the unperturbed Hamiltonian and the perturbation, which are both defined on the contour γ, and χ(z) is a switching function.In complete analogy with Eq. ( 13), H 1 (z) has different physical interpretations depending on the position of z on the contour γ: it represents an external perturbation of the physical Hamiltonian for z = t ∈ γ K , while it is a correction to the measurement operator or to the initial state when we choose z ∈ γ ↑,↓ or z ∈ γ M , respectively.
We will consider the work extraction in two scenarios.In the first scenario, the final Hamiltonian is equal to the initial one, i.e., χ(0) = χ(t f ) = 0. We will refer to this protocol as switching on/off, since the perturbation H 1 (z) is turned on after the first measurement and switched off before the second one.This is the setting in which Bochkov and Kuzovlev derived the first integral fluctuation theorem (Eq.( 8) in Ref. [53]) and the switching function in this case reads In the second scenario, the final Hamiltonian is arbitrary, i.e., χ(0) = 0 but χ(t f ) ̸ = 0.This is the general setting considered by Jarzynski [36] and the associated contour Hamiltonian is The results about the perturbation theory and the diagrammatic expansion of the MGF discussed in the next subsections are true for both the protocols (19) and (20), while at the end of Sec.4.2 we will continue our calculation by choosing specifically the protocol (19) and leave some comments about the generalization to the scenario (20) in App.D.

Interaction picture and non-interacting Green's functions
The customary approach to time-dependent perturbation theory begins with the introduction of the interaction picture [54].To generalize this concept to the modified contour we use the following relation where we denoted with z γ the integration between the initial point of the modified contour and the point z and defined with the contour interaction Hamiltonian defined as H1 (z) = U γ 0 (0, z)H 1 (z)U γ 0 (z, 0) (see App.A for details).The l.h.s. of Eq. ( 21) evaluated in −iħ hβ gives the numerator of (12), thus, after some manipulations we obtain where we used that H 0 (z) = H M on the γ M branch of the contour.Equation ( 23) can also be written as an explicit function of H 1 (z) The relevance of the equation above is due to the fact that when H 1 (z) is a small perturbation of H 0 (z), we can expand the second exponential in the numerator and obtain a perturbative series for the generating function.Before doing so let us focus first on the fundamental building block of perturbation theory which in systems of coupled bosons and fermions is given by the Green's function where c is a bosonic/fermionic annihilation operator and the angular brackets denote an average over the initial-state density matrix.Equation ( 25) is a straightforward generalization of the Keldysh contour GFs [19], but we stress that z, z ′ are defined on the contour in Fig. 1C instead of the standard Keldysh contour.Using the evolution in the contour interaction picture (see App. A), the dynamical equation for the Green's function reads where the boundary conditions are given by the Kubo-Martin-Schwinger relations G(z, 0) = ±G(z, −iħ hβ) for any z ∈ γ, where the upper or lower sign stands for bosons or fermions, respectively.
Let us now replace H 0 = ħ hωc † c.In this case Eq. ( 26) can be easily solved to obtain the non-interacting GFs (see Apps.B and B.1) where the upper (lower) sign again refers to bosons (fermions), n = 1 ± n, and n = n b, f (ω) = (e ħ hβω ∓ 1) −1 denotes the Bose (Fermi) distribution function and ω is the associated frequency.In Eq. ( 27) we also introduced the contour step function Θ γ (z − z ′ ), that is equal to 1 if z is placed after z ′ in the contour γ and equal to 0 otherwise.As in standard Schwinger-Keldysh theory, the contour Green's function does not have an evident physical meaning, unlike its components which are obtained by restricting the domain of z and z ′ to selected branches of the contour (we use ± subscripts to denote time variables on the branches γ ± ), for instance are called greater (>) and lesser (<) components, the dependence by λ arises from the vertical displacement in the complex plane between the γ − and γ + tracks.Other notable components can be obtained by choosing different contour arguments (see App. B.1), e.g., the time-ordered and anti-time-ordered components defined, respectively, as b, f (t + + iħ hλ, t ′ + + iħ hλ), which have the property of being λ-independent.

Perturbative expansion and diagrams
Within the perturbative framework, we can expand the second integral in the numerator of Eq. ( 24) in terms of correlation functions of arbitrary order, and obtain an expression of the form (see also Ref. [24]) At each order n we have a n-point correlation function of the Hamiltonian H 1 .If we assume that H 1 is a linear combination of products of fermionic and bosonic creation and annihilation operators, the correlation functions can be decomposed into non-interacting GFs following Wick's theorem [20,40].It is important to underline that the expansion of the correlation functions in terms of non-interacting GFs is the same as in standard perturbation theory because Wick's theorem is a consequence of the commuting (or anti-commuting) algebra of the operators c, c † and as such holds independently of the contour of integration.The discussion above can be reformulated in diagrammatic terms since the Feynman diagrams entering in the calculation of Eq. ( 29) are the same as in standard perturbation theory, the only difference being the contour of integration of the vertex variables z 1 , . . ., z n and the non-interacting GFs themselves, which are given by Eq. ( 27).To illustrate this, let us consider as an example the second-order perturbation theory of the Hamiltonian H 1 = ħ hω χ (a + a † )c † c where c and a represent, respectively, the annihilation operator of a fermionic and a bosonic mode.In this case the contribution of the second order (n = 2) to the numerator of Eq. ( 29) Figure 3: An example of how different Keldysh components appears in the diagram depicted on the left.The solid and wiggly lines denote the fermionic and bosonic non-interacting GFs, respectively.will correspond to two connected Feynman diagrams, the one in Fig. 3A and the "dumbbell" diagram (see Fig. 6 in App.C).Focusing on the former, which we call Assuming a switching on/off protocol as in Eq. ( 19) and dividing the integration over γ in an integration over γ − and γ + (the general rule to divide an integration over the contour in many integrals over the branches has been introduced by Langreth [55]), we rewrite the above integral in terms of the components of the GFs as In complete analogy to the connection between Eq. ( 30) and Fig. 3A, the four contributions to Eq. ( 31) can be expressed in terms of Feynman diagrams with a "charge" ± to take into account the position of the vertex variables on the contour, see Fig. 3B.
The time-ordered and the anti-time-ordered components are independent of λ while G <,> include a factor of e ±ħ hωλ which carries the frequency ħ hω [see Eq. ( 28)].This simple dependence allows us to introduce the energy transferred in a "charged" diagram d as E d = ħ h i s d i ω i where i runs over all the propagators and s d i = 0, −1, 1 depending on whether the propagator . Each diagram will contribute to the generating function a factor Γ d (t f )e λE d , where Γ d (t f ) is a prefactor that can be found by computing the contribution of the diagram (in the example of Eq. ( 31) by computing the associated double integral over the time variables).
It is convenient to introduce the cumulant generating function, given by the logarithm of Eq. ( 12), The diagrammatic expansion of each of the two logarithms contains only connected diagrams, as a consequence of the linked cluster theorem (see for instance Ref. [20], chapter 11).The second logarithm contains exactly the same diagrams of the first one, but has λ = 0.As a result we obtain the following general formula for the cumulant generating function where the summation variable d runs over all the connected "charged" diagrams.We stress again that the result ( 33) is true when assuming that the switching function is non-zero only in the horizontal branches, that is Eq. ( 19).The cumulant generating function ( 33) can be interpreted as a sum of independent rescaled Poisson processes with average Γ d (t f ) and jumps in energy given by E d .Therefore, every diagram represents a channel through which a quantized amount of energy E d can be exchanged with a Poissonian rate given by Γ d (t f ).The universality of Eq. ( 33) can be seen as a consequence of the Poisson law of rare events, since within the perturbative approach the rate of each energy exchange process is small.We expect Eq. ( 33) to break down when an infinite number of diagrams is resummed, in analogy to what happens in the case of charge statistics [26,56].
In App.C we explore several applications of the result (33), by computing the rates Γ d (t f ) explicitly for some specific models.In the case of the Holstein coupling H 1 = ħ hω χ (a + a † )c † c our results can also be verified by direct calculation of the cumulant generating function using the matrix form of the time evolution operator (see App. E).
We are now in a position to come back to our initial claim that the fluctuation theorem is satisfied at every order of the perturbative expansion.In Sec. 3 we showed that the contours used to compute M W (λ, t f ) and M rev W (−λ − β, t f ) are the same.Since the propagators (see Eq. ( 27)) and the integration domain of the vertex variables (see Eq. ( 29)) are fully determined by the contour, the two perturbative expansions with the same contour are identical.This ensures, as announced, that for the switching-on/off scenario (19) the symmetry ) is preserved at every perturbative order.

Time-reversed diagrams and detailed balance
We can investigate the effects of the FT at the level of the single Feynman diagrams.For this sake, we have to connect the diagrams appearing in the perturbative expansion of the moment generating function M W (λ, t f ), with the ones appearing in the expansion of its time-reversed counterpart M rev W (λ, t f ).The two MGFs share the same Hamiltonian, but the shape of the contour is different: while the vertex coordinates of the diagrams associated to M W (λ, t f ) live on γ in Fig. 2 A, in the case of M rev W they are defined on γ rev in Fig. 2 D. It is thus clear that the definition of the time-reversed diagrams goes through a mapping between the coordinates of the contours γ and γ rev .To understand the nature of this mapping, let us consider a "charged" diagram d appearing in the decomposition of the MGF in the scenario (20) (see e.g. one of the four contributions to Eq. ( 31)).We introduce its time-reversed as the diagram d appearing in the expansion of M rev W (λ, t f ) in which the sign of all the "charges" is the same as in d.Note that the positions of γ + and γ − are inverted in γ rev if compared to γ (see also [31]), thus for any greater GF G (0)> b, f (t, t ′ ) appearing in the formula for d, there will be a lesser GF G (0)< b, f (t, t ′ ) appearing in d and vice versa.
Since from Eq. ( 28) we have we can conclude that the ratio between the contribution of a "charged" diagram and its time-reversed contains a factor ±e βħ hω e −2λħ hω for every lesser GF and a factor ±e −βħ hω e +2λħ hω for every greater GF contained in the original diagram, where ω is the frequency of the fermionic/bosonic mode associated with the propagator.After introducing E d as in Sec.4.2 and setting λ = 0, we find that the weights Γ d (t f ) and Γ rev d (t f ) of a given diagram and its time-reversed are related by where the sign contribution ± in Eq. ( 34) can be neglected by assuming that the interaction Hamiltonian contains an even number of fermionic fields.Equation ( 35) can be seen as a diagrammatic version of the detailed balance conditions [57][58][59][60].Note that the condition (35) is stronger than the FT for the cumulant generating function, ), since it holds at a more detailed level, i.e., the level of single transitions.The FT at the level of the cumulant generating function can now be easily recovered by using Eqs.(35) and (33).To do so, we write the time-reversed generating function explicitly where E rev d is the energy transferred in the time reversed charged diagrams, that satisfies Using the detailed balance relation in (35), we replace 36) which gives To write the second term as a function of the forward rates Γ d (t f ), we note that since the second addend on the right hand side of Eq. ( 32) is equal to both log Z(0) and the λindependent term in Eq. ( 33), we have d,conn Γ d (t f ) = log Z(0).In the time-reversed generating function, Z(0) should be replaced by Z(t f ) (see Sec. 3), however, the two quantities are equal due to the assumption (19) (switching on/off scenario).Therefore, we conclude that d,conn Γ rev Replacing this into Eq.( 37) we conclude that

Path integration using the modified contour
Another useful approach to studying the generating function is to express Eq. ( 12) in terms of path integrals [22,23,61].This approach is an extension of the usual Feynman path integral approach on the Keldysh contour [41,62,63].We will see that our modified contour is particularly suitable for work statistics and for describing the semiclassical limit of the MGF by considering an expansion for small ħ h of the generating function [43].Let us consider the Hamiltonian of a single particle in an external potential, where P is the momentum operator, m is mass and V [α(t), X ] is a single particle potential in the particle position X , which depends parametrically on an external driving parameter α(t).
Plugging the Hamiltonian (38) into Eq.( 12) and performing a Trotter decomposition of the contour-ordered exponentials, we obtain the path integral form of the moment generating function, where p(z) and x(z) are the momentum and position fields defined on the modified contour γ, Z(0) is the partition function corresponding to the Hamiltonian at time 0 and S is the classical action Here we present the calculation for a measurement operator corresponding to the total energy, Λ = P 2 2m + V (X ), but the formalism works similarly for other choices of the measurement operator.We split the fields on the contour by defining their components, in a similar way to what was done for the GFs.Defining z = t + iτ we have and analogously for p(z).Eliminating the momentum operator p by integrating over Dp(z), one obtains the Lagrangian representation of the path integral (see App. F) where D ′ x(z) is the new measure of integration following the elimination of the momenta, and the Lagrangian action reads where L γ is the Lagrangian on the modified contour that according to the portion of the contour of interest can assume different forms as where we have This indicates that the Lagrangian on the vertical branches is the negative of the classical energy.The generating function in terms of the action in Eq. ( 40) can be used to study the characterization of the work fluctuations at the path integral level and its semiclassical limit.

Symmetrization of the contour and Keldysh rotation
Following the quantum-classical correspondence principle [64], the generating function in the semiclassical limit should reproduce its classical analog at first non-zero order in ħ h.This means that the path integral form of the generating function (42), in which fields are defined on the contour γ, should reproduce its stochastic path integral counterpart in a suitable limit [65][66][67].
The first obstacle that we have to overcome when trying to find such a correspondence is the difference between the domains of integration of the stochastic path integral, that is [0, t f ] and the Keldysh contour [6,68,69].In the path integral representation of the dynamics, this obstacle is simply overcome by the fact that the forward (γ − ) and the backward (γ + ) branches are equal, so the integration on the Keldysh contour can be seen as an integration over the segment [0, t f ] of the difference between the forward and backward actions [19].Applying this reasoning to the contour in Fig. 1C, we look for a transformation of the contour that makes it possible to divide it into two equal halves.This symmetrization is carried out by assigning half of the vertical lines of γ ↑ , γ ↓ and γ M to an upper branch named γ ⊕ and the other half of the lines to a lower branch γ ⊖ (see Fig. 4 for a detailed explanation of this procedure).In the symmetrized contour in Fig. 4C, an argument z ∈ γ ⊕ can be mapped to γ ⊖ simply by complex conjugation.This allows us to write the action in Eq. ( 43) as where in the last equality we separated the contribution of the real and imaginary parts of the differential dz.Interestingly, the contour in Fig. 4 has a great similarity with the the symmetric contour used by Aron et al. [49] to study the symmetries of the Schwinger-Keldysh action in equilibrium and non-equilibrium systems.We can now perform an analog of the Keldysh rotation [42][43][44]70] and introduce the classical and quantum fields as Inserting the above expressions into the action in Eq. ( 46), we can rewrite the generating function as a path integral over γ ⊖ .Since the integrand of the action in Eq. ( 46) depends on the branch (Im z nullifies on γ (⊖,−) while Re z nullifies on γ (⊖,↑) , γ (⊖,↓) and γ (⊖,M ) ) it is convenient to separate the contributions of the horizontal and vertical branches, obtaining where D ′ x cl/q (z) = D ′ x cl (z)D ′ x q (z) and we introduced the functions The γ ↑,↓,M branches of the contour in Fig. 1C can be naturally divided into two parts in this representation, which we denote by γ (⊖,↑) ,γ (⊕,↑) , by γ (⊖,↓) , γ (⊕,↓) and by γ (⊖,M ) , γ (⊕,M ) respectively.The branch γ (⊖,M ) goes from z 0 to z = 0, while the branch γ (⊖,↓) goes from z = 0 to z = −i ħ hλ 2 .
Notice that in Eq. ( 46) the domain of the fields of the path integral is γ ⊖ , with starting and ending points given by z = iβħ h/2 and z = t f .Since in these points the forward and backward fields are equal, the boundary conditions for Eq. ( 48) are x q (t f ) = x q (iβħ h/2) = 0.It is now natural, taking inspiration from the classical case in which the fluctuating work can be defined as a function of the endpoints of the stochastic trajectories, to introduce a quantum energy function, defined at time t f , as and its analogue at the initial time t = 0, where γ (⊖,↑) is replaced by γ (⊖,↓) .The difference between the initial and final energy functions gives a characterization of the fluctuating work at the trajectory level, and we can write the MGF as where M(z) and Σ(z) are given in Eq. ( 49) and The specific choice of the symmetrization shown in Fig. 4 is essential to obtain a representation of the MGF in which W (λ, t f ) depends only on the initial and final points of the trajectories.If we had chosen a different contour instead of considering the contour in Fig. 4C and its lower half γ ⊖ as the domain of integration in the last of Eqs.(46), W (λ, t f ) would have acquired a different functional dependence on the fields.A particularly interesting choice is the one used by Funo and Quan [22] that we will summarize in the next section.

Connection with the results of Funo and Quan
To compare our approach with the one presented in Ref. [22], we assume as a first step that the counting field is purely imaginary and replace λ → −iλ.This amounts to replacing the contour γ with the flat one in Fig. 5B.Such a contour allows a calculation of the characteristic function of work, which using path integrals can be written as where S 1 , S 2 and S M are the the actions of the forward, backward and γ M branches in the contour of Fig. 5B, respectively, The fields on the path integral (53) have boundary conditions x(ħ hλ , and x M (0) = y(0).These conditions are a direct consequence of the continuity of the field x(z) in the modified contour, since x, y, x M are its components in the contour in Fig. 5B (analogously to the case of the contour in Fig. 5A in which they are given by Eq. ( 41)).Note that differently from the contour in Fig. 5B, the Lagrangian in the forward and backward branches are now different.This lays behind the definition of the work functional as a difference between the forward and backward action in terms of the forward path x (for details we refer to Ref. [22]) which after some manipulations makes it possible to express the work functional at the level of quantum trajectories: where the integrand of the first integral acts as a quantum generalization of the instantaneous power.This result is different from the representation (52) in which the work is expressed as a function of the endpoints of the trajectories only.

Semiclassical limit
To analyze the semiclassical limit it is convenient to represent the path integral in the Hamiltonian convention, retaining the integral over the p fields until the end of the calculation.In this case, we can repeat the procedure in Sec.5.1 and obtain the same result as in Eq. ( 48) but with the path integral now including the momentum variables Dp x/cl .The functions M and Σ are replaced by (see App. G.1 for details) where K = 2p q (d x cl /dz) − 2x q (d p cl /dz).Since we expect the final energy to be a function of x cl (t f ), p cl (t f ) in the classical limit, we will assume that the contribution of x q and p q is small on the vertical branches.This intuition can be verified by keeping the nonlinear contributions in x q and p q , and showing that they scale as an higher power of ħ h (see App. G and [43]).To first order in p q and x q , Eq. ( 57) leads to . Therefore, at this order, the path integrals over x q and p q become trivial.Since the function K contains the first order of the quantum fields, the path integral on γ (⊖,↑) is written as Hence, the only admissible classical paths in this limit are the ones in which x cl (z), p cl (z) in the vertical branches are always equal to x cl (t f ), p cl (t f ).With this in mind, the path integral over the classical fields can be easily carried out, Combining this result with the one we obtain doing the analogous calculation on γ (⊖,↓) and γ (⊖,M ) , we finally conclude that that is the classical version of Eq. ( 52).This result is independent of the counting field λ and coincides with work fluctuations in classical driven isolated systems [36,71].We obtained Eq. ( 60) as a limit of Eq. ( 52), but we could have obtained it as the classical limit of Eq. ( 56) as well.In the latter case it can be shown (see [22]) that the term inside the time integral in Eq. ( 56) reduces to the instantaneous power output generated by the time-dependent driving and Eq. ( 60) is recovered after carrying out the time integration.
The first quantum correction to Eq. ( 60) can be obtained by keeping the quadratic terms in x q and p q in the expansion of Σ h .In this case we are left with a Gaussian integral over the quantum fields, which we computed explicitly in App.G.3.If we introduce where W 1 contributes to the variance of work, that corresponds to the second derivative of the generating function with respect to λ.Moreover, the above results become exact for timedependent harmonic potentials (see App. H).In this case, we recover the results of Ref. [72].

Detailed balance and connections with the Wigner function
One of the main advantages of the representation of the work ( 51), ( 52) arising from the contour in Fig. 4C is that while the information on the dynamics is all contained in the forward branch γ (⊖,−) , the λ dependence is isolated in the vertical branches.This separation between the dynamical contribution to the action and the "thermodynamical" one allows us to solve, at least formally, both contributions separately.This is not possible if we represent the work as in Eq. ( 56), since the fields appearing in the definition of W λ are the same appearing in the dynamical action (i.e., the fields in the forward branch).We can exploit this advantage to study a generalization of the concept of detailed balance at the level of quantum trajectories, in a way that is conceptually similar to what we did in Sec.4.3 in the case of the diagrammatic approach to perturbation theory.As a first step, we introduce the formal solution to the path integral on the vertical branch γ (⊖,↑) , that is the last part of Eq. ( 48) where the boundary conditions of the path integral are x q↑ (0) = 0, x q↑ (−ħ hλ/2) = y f q , and x cl↑ (−ħ hλ/2) = y f cl , and we made the dependence on α(t f ) explicit for convenience.Note that, in a similar way to Eq. ( 41), x q↑ (τ) and x cl↑ (τ) are defined as the quantum and classical fields on the branch γ (⊖,↑) , by replacing z = t + iτ.Since the contribution of the other vertical branch γ (⊖,↓),(⊖,M ) is functionally the same, apart from replacing λ → −λ − β, we can write Eq. ( 48) as with U representing the propagator from the initial values of the fields ( y i q , y i cl ) to the final values ( y f q , y f cl ), that can be obtained by integrating the contribution of the forward branch in Eq. ( 48) with the appropriate boundary conditions (see also Eqs. ( 62), ( 143), (144) of App.I for details).Let us define the integrand of Eq.( 48) for λ = 0 as It has a similar structure as the joint probability distribution of a classical system with a two-dimensional configuration space spanned by ( y q , y cl ), prepared in a Gibbs state with Hamiltonian E −β [ y i q , y i cl , α(0)], and following a stochastic evolution that brings the system in ( y f q , y f cl ).However, since U is in general a complex number, K lacks a probabilistic interpretation.Despite this, we can give an interpretation to Eq. ( 64) in terms of Wigner functions [42,43,[73][74][75].Indeed it is possible to prove that where W β is the Wigner representation of the initial state, Π is the Weyl transform of the propagator (see App. I) and C is a normalization constant.It is known that W β and Π on the r.h.s. of Eq. ( 65) reduce, respectively, to a classical Gibbs state distribution in the phase space, and to a classical propagator that imposes the deterministic equations of motion in the phase space [73].We can use both the integrand of the l.h.s. and the integrand of the r.h.s. in Eq. ( 65) to discuss the concept of detailed balance.Let us introduce K rev as the integrand of Eq. ( 48) for λ = 0, but in the symmetrized version of the time-reversed contour γ rev in Fig. 2D.In the time-reversed contour with λ = 0 the only vertical branch is placed at Rez = t f , and represents the initial preparation of the time-reversed trajectories.With this in mind, it is easy to prove that where e −β∆F = Z(t f )/Z(0) and the minus sign in front of the quantum variables at the denominator comes from the fact that the forward and backward branches are inverted in the time-reversed contour.As discussed above, K is not necessarily real and has no operational meaning in terms of measurements.On the contrary, if we choose y i q = y f q = 0 it represents the joint probability distribution P( y f cl , y i cl ) in a two-point measurement process of the position operator, with initial and final measured values given by y i cl and y f cl (see App. I).In the classical limit, this reduces to that is, a detailed balance condition for the exchanges of the potential energy of the system.This is expected since in the initial measurement there is no information on the initial momentum of the particle.To obtain the classical form of the detailed balance in the phase space representation, we have to rely on the r.h.s. of Eq. ( 65) instead.Indeed, even if the Wigner function is non-positive in general, it is in the classical limit, where we obtain where W 0 is given by Eq. ( 60), in analogy with the classical results [71], and W rev β denotes the Wigner function associated to the initial state of the time-reversed trajectory.The last equality also follows from the equivalence of the Wigner function and classical Gibbs state in the classical limit [73].The minus signs in front of the final momentum has been added for the sake of completeness, in our calculations the Hamiltonian is always quadratic in p and the change of sign becomes irrelevant.

Conclusions
We used the modified Keldysh contour as a versatile tool to investigate the perturbative approach to MGFs and their semiclassical limit.The symmetry property of the extended contour proved instrumental in showing the consistency of the perturbative expansion, ensuring that the fluctuation theorem is satisfied at every order in the expansion of the MGF.We showed that the modified contour technique is particularly suitable for diagrammatic approaches, and we used it to derive the work distribution in a switch on/off scenario, showing the work to be distributed as a linear combination of rescaled Poisson processes.Each of those processes represent a channel in which a discrete package of energy is exchanged between the system and the experimental driving apparatus.To investigate the semiclassical limit of the MGF, we expressed it in terms of a Feynman path integral and discussed how the form of the action depends on the choice of the modified contour.By symmetrizing the contour, we used the Keldysh rotation to write the action in terms of the classical and quantum fields, in a generalization of the Keldysh rotation approach.Our technique allows for a discussion of the detailed balance conditions, both at the level of the diagrammatic and the path integral approaches.In the former case, we show how the fluctuation theorem can be seen as resulting from a stronger detailed balance symmetry at the level of the single diagrams, while in the latter case we proposed a way to generalize the detailed balance to quantum trajectories and showed how this approach allows us to make contact with the Wigner function and with classical phase space approaches to thermodynamics.A interesting future perspective would be to generalize our approach to many-body open quantum systems [3,20] and assess the advantages of our new contour.This should be particularly relevant for preserving thermodynamic consistency while calculating work and heat counting statistics using known approximations schemes such as GW or random phase (RPA) [46,76,77].To conclude, we expect that generalizing our formalism to the case in which the system is in contact with many baths at different temperatures will present a very interesting challenge, since in this extended scenario the temperature is not unique, nor is the counting field, considering that one is typically interested in measuring many different thermodynamic fluxes (e.g. the different heat flows in each one of the reservoirs).We expect that studying in detail the many baths scenario could help in bridging the non-equilibrium GF formalism and the modified contour formalism with other general approaches to thermodynamic consistency in stochastic thermodyanmics [59].

A Interaction picture on the modified contour
The standard approach to time-dependent perturbation theory is based on the interaction picture [54].We start by proving Eq. ( 21).The derivative of U γ (z, 0) in respect to z reads The preceding equation can be shown by doing an infinitesimal displacement of the generator and expanding to first order Let us now compute the derivative of the right hand side of Eq. ( 21) and verify that is the same as Eq. ( 70).
and the last term is equal to the r.h.s. of ( 69) after regrouping H 0 (z) and χ(z)H 1 (z).Note now that the numerator of Eq. ( 12) in terms of U γ (z, 0) is simply given by Tr T γ e and we can use Eq.
After expanding Ũγ,1 as a contour ordered exponential and dividing by the proper normalization, we obtain Eq. ( 23) of the main text.

B Green's functions and higher order correlation functions
A generic two-point correlation function on the contour in the Schrodinger picture is written as To ease the calculations we rewrite the two-point correlation function as To take the the derivatives with respect to z 1 and z 2 , we use the equations bellow which come from ( 70) Using above equations we can take the derivative of ( 75) to obtain the equations of motion for the two-point correlation function More general correlation functions can be obtained by including more operators in the trace (74).For instance, we introduce the n-operator correlation function as In this manuscript we are mainly interested in the case in which the operators A i are products of bosonic or fermionic creation and annihilation operators.However, the equations of motion (77) can be solved only in specific cases, for instance when H γ is quadratic and A 1,2 are single creation/annihilation operators.At the contrary, when H γ is not quadratic, the equations ( 77) are not closed in G A 1 ,A 2 and depend on higher order correlation functions, giving rise to the Martin-Schwinger hierarchy [20].

B.1 Calculation of the non-interacting Green's functions
Using the results of App.B it is easy to find that for a bosonic operator a governed by the Hamiltonian H 0 = ωa † a, we have where the δ-term comes from the fact that, if z ′ > z, the ordering operator switches the position of a and a † .Substituting the definition of the Green's function into the equation above, we have for a bosonic Green's function, Table 1: The values of the contour propagator G (0) (z 1 , z 2 ) for selected values of the vertex time arguments.
We can derive the same equation for a fermionic operator c as follows, which is formally identical to Eq. ( 80).The solution of Eq. ( 80) for a generic GF reads The constant A in the equation above can be determined by means of the boundary conditions [20]: if we replace any of the arguments on the contour with the earliest and latest instants on the contour, the two results should be the same (for bosons) or differ by a sign (for fermions).
For a two point GF we thus obtain which leads to Eq. ( 27) of the main text.When λ = 0 and t = t ′ , Eqs. (28) represent the components of the non-interacting density matrix [40].Choosing instead both time arguments on γ − or γ + we obtain the time ordered and anti-time ordered GF where Θ is the Heaviside step function on the real axis.These functions coincide with the conventional Green's functions.A summary of the main GF components according to the position of the arguments on the contour is given in Table 1.It is possible to consider more general components by choosing the arguments in the new branches of the contour, e.g., a mixed non-interacting GF G (0) b, f (t, iτ) for t ∈ γ ± and iτ ∈ γ ↑ .These new GFs play no role in discussing the work statistics in the switching on/off scenario (19), but are relevant when considering the more general assumption (20).We will discuss them more in detail in App.D.

C Examples
Let us consider a simple Holstein coupling model, in which a fermionic energy level is shifted by an amount that depends on the position of a quantum harmonic oscillator.The Hamiltonian of this model is given by The Feynman diagram contributing to the second-order expansion of Eq. ( 29) for this potential is represented in Fig. 3A.The only other connected diagram contributing at second order is appears in every second order diagram, 2 comes from the freedom of exchanging the two vertices in fig 3, −1 from the single + in Fig. 3, i 3 from the definition of many body Green's function (see Eqs. ( 25) and ( 30)) and −1 is the factor prescribed by the Wick theorem (see for instance [20]).The energy jumps for these diagrams d = 3, 4 are E 3 = −ħ hω b and E 4 = ħ hω b , respectively, as we can verify directly by replacing the components in Eq. (86), By direct comparison with Eq. ( 33), we find with The contributions of the diagrams d = 7, 8 represented in Fig. 6 are the same as the ones computed for d = 3, 4 a part from the prefactor (−n 2 f + n f ) being exchanged with n 2 f .Summing all the contributions, we obtain In App.E we compute the generating function of the work exactly, and check that the expansion up to the second order in χ indeed leads back to Eq. ( 89).Another simple application is the anharmonic oscillator with Hamiltonian As usual, we consider the system as initialized in a Gibbs state of the unperturbed Hamiltonian, switch on the anharmonic term at time zero and perform energy measurements at times 0 and t, before switching on and after switching off the perturbation.The vertex in Eq. ( 90) is given by two inward or outward lines, that at second order produces a single diagram (a loop with two bosonic propagators between time t 1 and t 2 ).The contribution of such a diagram gives the cumulant generating function at second order In the notation of Eq. ( 33), this cumulant generating function consists of two rescaled Poissonian energy jumps with energies equal to E d = ±2ħ hω b and rates given by

D The role of the other GF components on the contour
When writing the contour GF in its component we also need to consider situations in which one of the arguments is in the vertical track γ ↑ .We can define a non-interacting GF with both arguments on the vertical track, Another possibility is represented by the case in which one argument is on the horizontal branches and one on γ ↑ .We can identify four different cases G We explicitly write this four new components below This definitions pave the way to an application of the theory to protocols beyond the switching on/off paradigm in which χ(z) is defined as in Eq. (20).For a generic integral of a function in the modified contour we have where f ↑ , f ± are the components of the function f in the modified contour.Notice that the contribution of the downward and Matsubara tracks is absent, since we are implicitly assuming χ(z) = 0 on γ M ,↓ as prescribed by the protocols (19) and (20).The same decomposition of the integral can be done for the equations arising from the perturbative expansion, like Eq. ( 30).
To make the calculations simpler, let us carry them on in the case in which the switching function is given by this case correspond to a system in which the Hamiltonian ω χ H 1 is quenched and then immediately measured (without leaving the system enough time to evolve).We can compute the generating function for the case of the Hamiltonian (85) discussed in the main text.The diagrams involved are the same (see A panels of Fig. 3 and 6) but we integrate over the vertical branch, and the only Green's function appearing in the calculations is G (0)↑↑ b, f .The integrals in Eq. ( 30) reduce to Where we denoted with d = 1 the contribution of the diagram in Fig. 3A with the value of the vertex variables to be chosen in γ ↑ .The contribution of the dumbbell diagram is instead given by The cumulant generating function at second order can be obtained by summing the two contributions in ( 96) and ( 97).

E Work statistics in the dispersive coupling case: non-perturbative approach
The dispersive coupling Hamiltonian is common to several models in open quantum systems, where a and c denote a bosonic and fermionic annihilation operator, respectively.The fermionic operators in the Hamiltonian (98) can be expressed as Pauli matrices acting over a two dimensional Hilbert space, by defining σ + = c † , σ z = 2c † c − 1.The anticommutation relations {σ − , σ + } = {c, c + } = 1 are preserved and we obtain The Hamiltonian commutes with σ z , so we can write the time evolution operator as The top-left element of the matrix (100) is the same as from a driven harmonic oscillator.We can define H b = ħ hω b (a † a + 1/2) and the exponential above can be written in terms of displacement operators [78] e where θ (t) is a phase factor and is the displacement operator of argument δ(t).The time dependent parameter δ(t) is connected to the classical solution for the Hamilton equations for the variable α(t) = x(t)+ip(t) 2 , with α(0) = 0.This variable evolve as α(t) = αe −iωt + δ(t), where δ(t) = so that it is clear using Eq. ( 100) that we have to evaluate the tilted displacement D[δ(t f ), λ] = e λH 0 D[δ(t f )]e −λH 0 .Since e λH 0 e −λH 0 = 1 we can bring the two matrices at the exponent in the displacement operator and obtain Using Eq. ( 104)) we find We are interested in a weak perturbation expansion of equation above, that is δ(t f ) << ω b , so we can Taylor expand D[−δ(t f )]D[δ(t f ), λ] keeping only the terms with the same number of a and a † operators, the others averaging to 0: The final expression for M (λ, t) thus reads computing the logarithm to find the expansion for the cumulant generating function and noting that 2 sin 2 ( 2 ) = 1 − cos ω b t f we obtain and confirm the result of the Eq.(89).

F Summation over the momentum variables in the path integral
To write the path integral in terms of the position coordinates we consider the quantity S = i/ħ hS in the exponent of the Eq.(39).From the main text we have where for ease of notation we exclude the parameter α from the argument of V .Depending on the branch we are considering the integration of the momentum variables leads to different results.Adopting discrete notation, the exponent of the path integral reads (112) The terms above give the contribution to exponent in the path integral and can be written respectively like i T, −i T, T, −T with T the kinetic energy.After going back to the contour variable z this kinetic contribution sums with the potential energy and gives In total we can write the exponent of the path integral as where L is defined on the contour as

G Semiclassical limit G.1 Hamiltonian action in the modified contour
The same manipulations done to obtain (46) can also be done before the integration of the momentum variables, at the level of the action (40): where we used that in the vertical branches, in which d Im z ̸ = 0, we have d dz * = − d dz , while in the horizotontal branches, in which d Re z ̸ = 0, we have d dz * = d dz .
To have an expression of the action in terms of the quantum and classical components of the fields we apply the linear transformation (47) and substitute From it, using again Eq. ( 47), we obtain an expression only in terms of the quantum and classical components If we replace the equation above in Eq. ( 116) we obtain the action in terms of the quantum and classical variables i

G.2 Expansion of the action with respect to the fields
In this section, we perform the expansion of the action (119) under the assumptions that the quantum fields admit an expansion in powers of ħ h.Now we can expand both the integrals over the quantum variables x q and p q up to the second order thus for the horizontal track we will have and for the vertical branches After integrating by parts Eq. ( 119) and replacing the expansions (120) and (121) we obtain Notice that the boundary terms of the integration by parts vanish as an effect of the boundary conditions x q (t f ) = x q (iħ hβ/2) = 0. Let us decompose the integral over the vertical lines ( γ ⊖ d Im z) in two separate contributions given by the vertical tracks γ (⊖,↑) and γ (⊖,↓),(⊖,M ) where we write the components of the quantum and classical fields as Note that in (41), we identified the components of the branches γ M and γ ↓ by x M (τ) and x ↓ (τ), respectively.However, here on the symmetrized contour, for ease of calculations, we have denoted the fields on the branches γ (⊖,↓) and γ (⊖,M ) by x cl/q↓ (τ).In this way we will have

G.3 Gaussian integration of the fields
In the Eq. ( 124) the action is quadratic in the fields x q , p q .They can be eliminated through a Gaussian integration.If define i ħ h S as the exponent of the path integral after such Gaussian integration, we have Notice that a byproduct of the path integration above is a change of the normalization factor proportional to 1/ To make further manipulations possible we will assume that the quantities V ′′ [x cl↓ (τ)], V ′′ [x cl↑ (τ)] can be respectively replaced with in the kinetic terms of the action, and the normalization.The idea behind this replacement is that the vertical tracks are short in the semiclassical limit, so that it is possible to replace the values of the fields in γ ↑ , γ ↓ with the values of the fields at their boundaries.Let us focus now on the integral in the momentum variables.After defining m(t where for now we restrict our analysis to γ ↑ .To integral over p cl↑ (0) reads Similarly for the other vertical branch we will have A further integration over p cl↓ (βħ h/2) yields Thus we can write the contribution to the path integral due to vertical lines as ] is exact.Therefore using (130) and doing the path integral over x cl we will have D x cl D x q Dp q D p cl e To go further we need to integrate over p(t f ), x(t f ), p cl (0) and x cl (0), However, the relation between them is given by the path integral over the horizontal line that results in the equation of motion for a time dependent harmonic oscillator as x(t f ) = A(t f )x(0) + B(t f )p(0)/M and p(t f ) = M Ȧ(t f )x(0) + Ḃ(t f )p(0).Therefore we can write the above relation as where the term sinh βħ h/2 comes form the normalization of the generating function.Therefore the integration over p 0 and x 0 will result in 2π/ det(C) which can be written as M W (λ, t f ) = sinh(βħ h/2) cosh(ħ hω 1 λ/2) cosh(ħ hω 0 (λ + β)/2) (139) We notice that the last term under the square root is the Wronskian and for our given boundary conditions is 1.Also for Q * we have Thus we write (139) as This allows us to recover the results of ref [72].

I Classical limit in detailed balance conditions
In Sec.5.4 we introduced a notion of detailed balance at the level of the quantum trajectories.The strategy of Sec.5.4 is to divide the path integration in Eq. ( 48) in three separate propagators, that read where B( y i ), B ′ ( y f ), B ′′ ( y i , y f ) are shortcuts to denote the boundary conditions of the path integrals.For Eq. ( 142) the boundary conditions are x q↑ (−ħ hλ/2) = y f q , x q↑ (0) = 0, x cl↑ (−ħ hλ/2) = y f cl , for Eq.(143) we have 1 x q↓ (ħ hβ/2) = 0, x q↓ (−ħ hλ/2) = y i q , x cl↓ (−ħ hλ/2) = y i cl while for Eq.(144) we have x q (0), x cl (0) = y i q , y i cl , x q (t f ), x cl (t f ) = y f q , y f cl .The definition of the time-reversed trajectory is arbitrary, but similarly to what we did in Sec.4.3, for every trajectory defined in the Keldysh contour γ K , we can consider the trajectory that attains the same values but on a reversed contour in which the forward and backward branches are exchanged.This trajectory starts in the final points ( y f cl , − y f q ) and ends in the initial points ( y i cl , − y i q ), where the sign of the quantum components has been changed due to the inversion of the forward and backward branches.We have U( y f q , y f cl , y i q , y i cl ) = U rev (− y i q , y i cl , − y f q , y f cl ) from which Eq. ( 66) follows.Let us now focus on the classical expansion of the quantity E −β , that is the integral in the branch γ (⊖,M ) .The action in this case is given by Σ, that can be expanded for small x q↓2 as Σ = m ẋ2 cl↓ + ẋ2 q↓ + 2V (α, x cl↓ ) + V ′′ (α, x cl↓ )x 2 q↓ + O(ħ h 2 ). (145) Note that respect to Eq. ( 49) there is a change of sign in the kinetic term, due to the use of x cl↓ (τ) instead of γ cl (z), since we have d dz = −i d dτ in the vertical branches.Since the branch γ (⊖,↓) is short for ħ h → 0, the value of x cl↓ will not variate too much from the value attained at the boundary with γ (⊖,−) (that is y i cl ) so we can define x cl↓ (s) = y i cl + δ(s) and obtain Σ =m δ2 (s) + ẋ2 q↓ (s) + 2V (α(0), y i cl ) + 2V ′ (α(0), y i cl )δ(s) + V ′′ (α(0), We are mainly interested in the case in which the boundaries in the path integration in Eq. ( 143) with λ = 0 are given by x q↓ (ħ hβ/2) = x q↓ (0) = 0, a case in which the contribution of the integral over the quantum variables is negligible.The path integration in x cl↓ reduces to a path integration in δ(s) with boundary conditions δ(0) = 0, while the initial value of δ is free, we call it δ: dτ m δ2 (s)+2V (α(0), y i cl )+2V ′ (α(0), y i cl )δ(s)+V ′′ (α(0), y i cl )δ 2 (s) where we introduced Ω = V ′′ (α(0), y i cl ) m and ∆ = V ′ (α(0), y i cl ) V ′′ (α(0), y i cl ) .After doing the change of variable δ → δ + ∆ and solving the path integral (remember that changing the variable also changes the boundaries of the path integration) we obtain as a result e −β[V (α(0), y i cl )− where the two contributions in the last exponential come from the zeroth and second order expansions of the hyperbolic cosine in respect to ħ h.Now note that the terms proportional to ∆ 2 cancels out so we are left with the integral over δ where the second term in the integrand in the first line can be neglected since it is of the next order in ħ h.This proves the result (67).To prove the last result it is sufficient to use the properties of the Wigner function for small ħ h, we can focus on Eq. ( 65).
d y i q d y f q K( y f q , y f cl , y i q , y i cl ) = d y i q d y f q d p i d p f dη i dη f K(η f , y f cl , η i , y i cl )e ip i ( y i q −η i ) e ip f ( y f q −η f ) .
Let us focus on the dependence by η i , y i q , p i , after replacing the definition of K in the equation above we have a contribution of the form Z(0) −1 dη i d y i q e ip i ( y i q −η i ) U[η f , y f cl , t f ; η i , y i cl , 0]e −βE −β [η i , y i cl ,α(0)] .
Remembering the definition of E −β we can also write e −iη i p i Z(0) −1 e −βE −β [η i , y i cl ,α(0)] = e −iη i p i y i cl − η i ρ 0 y i cl After performing the integral in η i , the quantity above becomes (up to irrelevant prefactors) the Wigner function associated to the initial state, that we denoted with W β .In a similar way, we can show that the integral over y i q and the integrals over y f q transform the propagator U in the quantum phase space propagator.

Figure 1 :
Figure 1: Three different time integration contours: A) The Keldysh contour γ K , in which γ − and γ + denote, respectively, the forward and the backward branches.B) The contour γ C , obtained by adding two vertical tracks γ ↑,↓ to the Keldysh contour.C) The contour γ, realized by augmenting the previous one with the track γ M encoding the initial condition.

Figure 2 :
Figure 2: A graphical representation of the proof of Eq. (16).The correspondence between exponentials of time-dependent Hamiltonians and contour tracks allow us to establish a handy parallelism.Multiplying and dividing by exp[−β H(t)] is equivalent to adding and subtracting a new track (A → B) in the picture.Exchanging −λ with β + λ can be expressed as an inversion of the blue/green stars with the blue/green circles (B → C).The cyclic property of trace tells us that the starting point of the contour is arbitrary, so that if we mark the corners of the contour with 1, 2, 3, 4, we can draw an equivalent contour starting from point 3 and running clockwise, instead of starting from point 1 (C → D).Note that from point C to D we also shifted the contour upwards.The last step assumes time-reversal symmetry of the Hamiltonian and shows that the contour D is indeed the one associated with the backward evolution of the original one (D → A).

Figure 4 :i ħ hβ 2 + i ħ hλ 2 ,
Figure 4: The equivalence between the contour of Fig. 1C and a new symmetric contour.A) The last interval [−i ħ hβ 2 + i ħ hλ 2 , −iħ hβ], denoted with γ in the picture, can be removed and attached to the initial part of the contour using the cyclic property of trace.The result of this operation is represented in the next panel.B) We can translate along the imaginary axis by ħ hλ 2 using the invariance of the generating function under such a translation, thus obtaining the contour in the next panel.C) The final symmetric version of the contour.In this picture, we have used z 0 = i ħ hβ 2 andz 1 = t f + i ħ hλ 2 .The γ ↑,↓,M branches of the contour in Fig.1Ccan be naturally divided into two parts in this representation, which we denote by γ (⊖,↑) ,γ (⊕,↑) , by γ (⊖,↓) , γ (⊕,↓) and by γ (⊖,M ) , γ (⊕,M ) respectively.The branch γ (⊖,M ) goes from z 0 to z = 0, while the branch γ (⊖,↓) goes from z = 0 to z = −i ħ hλ 2 .

Figure 5 :
Figure 5: Transforming the modified contour to the asymmetric one by changing λ → −iλ.A) the modified contour.B) Asymmetric contour resulting from the transformation.One can use this contour to define the work functional in terms of the forward paths [22].

Figure 6 :
Figure 6: Dumbbell diagrams appearing in the calculations of the work statistics for the Hamiltonian (85).

χω χ ω b ( 1
− e −iωt ).For counting the work statistics in the sudden quench scenario we have the following generating function