Quantum work statistics at strong reservoir coupling

Calculating the stochastic work done on a quantum system while strongly coupled to a reservoir is a formidable task, requiring the calculation of the full eigenspectrum of the combined system and reservoir. Here we show that this issue can be circumvented by using a polaron transformation that maps the system into a new frame where weak-coupling theory can be applied. It is shown that the work probability distribution is invariant under this transformation, allowing one to compute the full counting statistics of work at strong reservoir coupling. Crucially this polaron approach reproduces the Jarzynski fluctuation theorem, thus ensuring consistency with the laws of stochastic thermodynamics. We apply our formalism to a system driven across the Landau-Zener transition, where we identify clear signatures in the work distribution arising from a non-negligible coupling to the environment. Our results provide a new method for studying the stochastic thermodynamics of driven quantum systems beyond Markovian, weak-coupling regimes.

In macroscopic thermodynamics, the energy of the interaction between a system and a thermal environment is dwarfed by their respective bulk energies.Consequently, the interaction contribution to energy changes like heat and work can be safely ignored, and the equilibrium state of the system can be adequately described by a thermal Gibbs state at the same temperature as the environment.However, the weak coupling approximation becomes less well founded when we consider the thermodynamics of mesoscopic and quantum mechanical systems in which the interaction terms can become more comparable to the system's energy scale [1].As recent advancements in quantum control now makes this regime experimentally accessible [2][3][4], there is significant interest in exploring how non-negligible interactions and structured environments can modify the equilibrium behaviour and thermodynamic laws of strongly coupled quantum systems [5][6][7][8][9][10][11][12].
The first key deviation from standard thermodynamics lies in the fact that the equilibrium state of the system generally differs from the familiar Gibbs state when weak coupling breaks down.Studies of both exactly-solvable and approximate models of strongly coupled quantum systems provide evidence that, when left to equilibriate, the system reaches the partial trace of a global Gibbs state with respect to both system and environment [13][14][15][16][17].This motivates the introduction of the so-called mean-force free energy F * S [18][19][20][21][22], which is an effective thermodynamic potential for the system that determines its properties at equilibrium while encoding information about interactions with the environment.This generalised notion of free energy also acts as the main ingredient for extending the laws of stochastic thermodynamics to the strong coupling regime out of equilibrium.For example, one can derive an extension of the Jarzynski fluctuation theorem when a system is driven away from a non-Gibbsian equilibrium state [23]: Here W denotes the work done on the open system, .denotes an average with respect to the work distribu-tion P (W ), and ∆F * S denotes the change in the meanforce free energy of the system at inverse temperature β.While formally correct, the difficulty in verifying this theorem lies in the fact that the work distribution P (W ) is formidably difficult to compute in strongcoupling regimes since it is obtained from projecting the total system and environment into global energy eigenstates at the start and end of the process [12,23,24].Ultimately this has led to work statistics at strong coupling receiving little attention, with the exception of relatively few exactly solvable models [25][26][27][28].One significant issue is that numerically exact techniques for simulating open system dynamics, such as discrete time pathintegral methods [29][30][31], are formulated to describe the dynamics of the reduced state of the system, and therefore do not have straightforward access to the eigenspectrum of the composite system and environment required to construct the full work distribution.
In this paper we introduce an alternative approach based on the polaron framework [32][33][34] that is not hampered by these issues.The polaron method uses a unitary transformation to dress system states with vibrational modes of the environment, capturing the dominant non-Markovian effects and allowing the system-environment interaction to be included non-perturbatively.The polaron framework has been successfully adapted to study resonant energy transfer in photosynthetic complexes [35,36], exciton-phonon interactions in semiconductor quantum dots [37], and heat exchange statistics in the time-independent spin-boson model with multiple reservoirs [38,39].Within the polaron transformed frame, we are able to derive a generalised master equation (GME) that governs the evolution of the characteristic function associated with P (W ) beyond weak coupling regimes.This can be used to infer the full statistics of work and reproduces the fluctuation theorem (1) as required for thermodynamic consistency.To demonstrate our approach, we consider a Landau-Zener (LZ) driving protocol on the system Hamiltonian.We find that stronger coupling creates new peaks in the distribution as well as imparting significant broadening and renormalisation, indicating clear signatures in the work statistics arising from interactions with the reservoir that are either not captured or severely underestimated by a weak-coupling approach.
Generalised Master Equation for Quantum Work-We begin by considering a general time-dependent spinboson model, governed by the Hamiltonian, where the time-dependent system Hamiltonian is given by Here ω 0 (t) is a time-varying energy spacing and ∆ is a fixed tunneling coefficient.The reservoir Hamiltonian is given by , with b k the annihilation operator for the kth mode of frequency ω k .The final term in Eq. 2 describes the system-resevoir interaction mediated by a set of scalar coupling constants g k .This interaction is fully characterised by the spectral density J(ω) = αω −2 c ω 3 exp(−ω/ω c ), with cutoff frequency ω c and α a measure of the interaction strength.The composite state at a time t is given by the density operator ρ(t) = U (t)ρ(t i )U † (t) with initial state ρ(t i ) and timeordered unitary ti dt H(t ) .We now wish to compute the statistics of work done on the system during the driving protocol, that is we are interested in the total energetic change of the system and environment during the protocol, including any stochastic fluctuations of these quantities.Independent of the coupling strength, the work done can be defined using the two-point measurement protocol, where a projective energy measurement of the global Hamiltonian ( 2) is made at the initial and final times of the protocol, with the work done W defined as the resulting energy difference [23].The resulting distribution P (W ) can be characterised by its characteristic function (CF), given by its Fourier transform Here we define the work characteristic operator (WCO) [40], where ρ(t i ) = n Π i n ρ(t i )Π i n is the initial state dephased in the energy eigenbasis of H(t i ), with Π i n a projector onto the nth eigenstate.Without further approximations, computing the CF even numerically is challenging due to the presence of the interaction term in (2).One immediate solution would be to assume the interaction is sufficiently weak, allowing one to derive a Lindbladlike master equation for the WCO (5) [41][42][43][44].However, our aim is to probe regimes where this approximation is not appropriate, and so we require an alternative method for computing K(t, η).Our approach is to use the fact that we are free to apply a trace-preserving operation to the WCO such as K(t f , η) → K P (t f , η) = e P K(t f , η)e −P .This will always leave the CF invariant so that Φ(η) = tr{K P (t f , η)}.For our spin-boson model we choose the unitary polaron transformation [34], where . This transformation redraws the boundary between the system and environment by dressing the system states with vibrational modes of the environment.By moving into the polaron frame, we can account for a significant part of the interaction using these dressed system states, and the residual part of the interaction in this frame can be treated perturbatively.The transformed WCO is given by K ti dt H P (t ) is the evolution operator with respect to the polaron-transformed Hamiltonian Here we identify the polaron system Hamiltonian where 0 ≤ κ ≤ 1 denotes the polaron renormalisation constant that lowers the tunneling coefficient of our transformed system.The degree of renormalisation is dependent on the structure of the reservoir, and can be computed from its spectral density according to Furthermore, we have identified a new interaction term in the polaron frame given by where . This new interaction can now be treated in a perturbative manner quantified by the parameter [34] Taking a weak-coupling approximation in this new frame will allow us to determine the work statistics of driven systems that are strongly coupled to the reservoir with respect to the original frame.
Let us assume that the system and environment have initially thermalised to a global Gibbs state with respect to the Hamiltonian in the original frame, ρ(t i ) = π eq (t i ) := e −βH(ti) / tr e −βH(ti) with β = 1/k B T an inverse temperature.As our first approximation we suppose that the interaction strength in the polaron frame is weak enough relative to the reservoir temperature so that β 2 g 2 1.We can then neglect any correlations in the initial state and contributions of the interaction to the initial and final energy measurements, so that the polaron transformed WCO simplifies to × π eq P S (t i ) ⊗ π eq R U † P (t f ), (11) where π eq R is a Gibbs state of the bare reservoir with Hamiltonian H R while π eq P S (t i ) is a Gibbs state of the polaron system with the initial Hamiltonian H P S (t i ).We emphasise that this approximation is valid in strong coupling regimes, unlike the typical weak coupling limit in which ρ(t i ) π eq S (t i ) ⊗ π eq R [24].To obtain the work statistics we derive a master equation for the WCO of the system degrees of freedom in the polaron frame, which is obtained by taking a partial trace over the reservoir so that K P S (t f , η) = Tr R K P (t f , η) .The derivation of an equation for the evolution of K P S (t f , η) then follows similar steps used to obtain an adiabatic master equation for the polaron state itself, and details are provided in Appendix A. This rests on a number of assumptions such as (i) the Born-Markov approximation, stating that the coupling in the polaron frame is much smaller than the inverse of the decay time of the reservoir correlation function and (ii) the adiabatic approximation, which ensures that the system Hamiltonian does not change too rapidly in time relative to its energy gaps and bath correlation time.Note that applying the Born-Markov approximation in the polaron frame is not equivalent to applying it in the original frame, and retains the dominant non-Markovian effects in the dynamics.Under these assumptions the general form of equation for the polaron WCO is then (12) where we define the work superoperator W(t, η)(.) = ∂ t e iηH P S (t) e −iηH P S (t) (.), (13) and L t (.) = −i[H P S (t) + H P LS (t), (.)] + D t (.) is an adiabatic Lindblad generator with H P LS (t) the polaron lamb shift and D t (.) the dissipator, which are both derived and defined in Appendix B. Physically, the dissipator drives transitions between the instantaneous energy eigenstates of the polaron system, defined as |ε with θ t = arctan(∆κ/ω 0 (t)) and energy eigenvalues ε ± (t) = ± 1 2 ω(t), where ω(t) = ω 0 (t) 2 + κ 2 ∆ 2 is the transition frequency.This contrasts with the standard weak-coupling master equation (WCME) for timedependent driving [45], whose dissipator causes transitions between eigenstates of the original system Hamiltonian H S (t) (see Appendix B).
A key property of this generator is that it has an instantaneous thermal fixed point so that L t [π eq P S (t)] = 0; π eq P S (t) = e −βH P S (t) tr P S e −βH P S (t) .(14) Taking the trace of the solution to (12) for a given protocol with an initial condition set to K P S (t i , η) = π eq P S (t i ) finally allows us to calculate the work CF.As a main consistency check, our formalism can be used to recover the quantum fluctuation theorem in the strong coupling regime [23].In particular, one may combine the fixed point property (14) with our chosen initial conditions to show that K P S (t f , iβ) = e −βH P S (t f ) /tr P S e −βH P S (ti) is a solution to (12).This immediately implies the quantum Jarzysnki equality e −βW = e −β∆F P S , ( where ∆F P S = F P S (t f ) − F P S (t i ) is the change in the polaron system free energy Here we may interpret this quantity as an approximation to the original system's mean-force free energy [1,22,23], defined by which is the difference between the total free energy F SR (t) = −β −1 ln tr e −βH(t) of the composite state and the free energy of a bare reservoir, F R = −β −1 ln tr R e −βH R .It then follows from our approximation β 2 g 2 1 that the polaron free energy is equal to this mean force free energy to the same order in coupling; F P S (t) F * S (t) and hence we have recovered the Jarzynski equality (1) obtained within the original frame.From the validity of the fluctuation theorem, one can therefore conclude that the work distribution predicted by our equation ( 12) is consistent with the laws of stochastic thermodynamics [12].
It is important to highlight the parameter regimes for which our model is valid; to do this, we shall consider the specific case of a two level system driven over a Landau-Zener crossing.In this case, the time dependent driving in the Hamiltonian may be written as ω 0 (t) = νt, where ν is the rate at which the driving occurs.First we consider the Born-Markov approximation in the polaron frame, which amounts to assuming that the environment can respond instantaneously to the dynamics of the system, where the dynamical timescale of the reservoir is proportional to the cutoff frequency ω c .From (10) this requires [34] So long as we have a sufficiently large cutoff frequency ∆ ω c , this condition can be satisfied from low to high temperatures, and from weak to strong coupling as defined by the original frame (small to large α).This provides us with a much wider parameter regime to probe the work statistics than would be possible with the WCME.Second, to ensure that the adiabatic approximation holds for the time-dependent protocol, we need the rate of change of the polaron's eigenbasis to be small relative to the energy gaps of its Hamiltonian.This is achieved when where t i and t f are the initial and final times respectively [46].The maximum is found at the intermediate avoided crossing where t = 0.If we compare this condition with the Born-Markov condition (18), we see that this can be maintained at any temperature so long as we choose a large enough cutoff frequency ω c .As a final condition we also require that the rate of change of the polaron energy eigenbasis must be small relative to the reservoir dynamics [45,[47][48][49], which means ν 2ω c ; this is guaranteed if the previous conditions are satisfied.To confirm that polaron theory gives a faithful description of the reduced system, we have benchmarked the system population dynamics using the numerically exact TEMPO method [30] in Appendix C and find excellent agreement with the adiabatic polaron dynamics at strong coupling, in contrast to the WCME.Therefore, we conclude that in the Landau-Zener model our master equation for the WCO Eq. ( 12) is accurate so long as ∆ ω c and ν 2∆ 2 κ 2 .Work Statistics-We now numerically solve the master equation (12) for the Landau-Zener model.In all following calculations, we choose ν = 0.1∆ 2 and start and end the protocol far from the avoided crossing −t 0 = t f = 100, this ensures that the adiabatic condition is satisfied.To ensure the validity of the PME we work in the scaling limit ∆/ω c 1 [37,50], taking ω c = 10∆.To calculate the work distributions we sample the CF via the solution to the generalised PME (12) over a range of countingfield values at intervals δη up to some cut-off η max , before performing a numerical inverse Fourier transformation on this data to obtain the work probability density function.A discrete work probability distribution is then found by numerically integrating the probability density over fixed intervals δW .
Before considering the work statistics of the full open quantum system, it is worth briefly discussing the expected behaviour for a closed system.If we start the system in the thermal state, such that the initial thermal occupation of the ground state is p − = 1/(1 − e −βω ), then we expect there to be three possible contributions to the work distribution: if the system's trajectory remains adiabatic there is no net change in the energy of the system, and therefore no work is done.This occurs with probability P (W = 0) = 1 − P LZ , where P LZ = 1 − exp −π∆ 2 /2ν is the LZ transition probability.If a LZ transition occurs, as shown in Fig. 1(i), the system can move from its ground to excited state.The work done on the system is given by the net energy change W = ∆E = (ω 0 (t f ) − ω 0 (t i ))/2, and occurs with probability P (W = ∆E) = p − P LZ .As the system is initially in thermal equilibrium, the converse process may also occur, where a LZ transition occurs between the excited and ground states, extracting W = −∆E of work.This occurs with probability P (W = −∆E) = (1 − p − )P LZ , and is expected to be heavily suppressed if the LZ protocol is initialised far from the avoided crossing.Now let's consider the more interesting case of the open quantum system.In Fig. 1 we show polaron and weak-coupling predictions for the quantum work distribution for a choice of α well into the strong coupling regime.The most immediate difference in the open system distribution is that it is continuous in contrast to the discrete distribution of the closed system.Since the bath is composed of a continuum of modes, it may exchange energy with the system at any of the continuous range of transition frequencies that the system moves through during the protocol.
In analogy to the closed system, both the polaron and weak-coupling distributions show a peak at W = 0 and W = ∆E corresponding to the adiabatic and nonadiabatic system trajectories (again, see Fig. 1(i) for a schematic of the latter process [51]).Interestingly, nonadiabatic transitions are over two orders of magnitude more likely to occur in the polaron formalism versus the weak-coupling approach.This is a direct consequence of the renormalisation of the energy gap to ∆κ, with κ given by Eq. ( 8), meaning that the system dynamics must be slower to maintain the same degree of adiabaticity.Unlike the closed system, both the weak-coupling and polaron probability distributions show a much richer structure due to phonon mediated processes.The simplest phonon process is shown schematically in Fig. 1(ii), where a phonon excites the system over the avoided crossing, which then continues to evolve adiabatically in the excited state.This requires W = ∆E − ∆ ≈ 9ω(t f ) of work in weak coupling theory, which is reduced to W = ∆E − ∆κ ≈ 9.5ω(t f ) for the polaron approach due to the renormalisation of the avoided crossing.
The peaks in the distribution either side of the origin, which are due to a heat exchange event at the avoided crossing in combination with a nonadiabatic transition (Fig. 1(iii)), can be interpreted similarly and are brought closer together in the PME due to renormalisation.While renormalisation tends to shift such peaks, we also observe that strong coupling tends to increase the likelihood of larger amounts of work being done on the system.One explanation for this is routed in the fact that, as mentioned, non-adiabatic transitions are less suppressed in the PME, as can be seen from its larger adiabatic parameter (19) when compared with the WCME (for which κ = 1).From a thermodynamic perspective, this leads to a greater probability of excitation and thus an increased likelihood of dissipation, which here corresponds to larger amounts of work being done to the system.As a final observation, the polaron approach predicts another type of process not apparent in the weak coupling theory and given schematically by Fig 1(iv), whereby a LZ transition excites the system, which emits a phonon, before a second LZ transition reexcites the system.The total work done in this situation is W = ∆E + ∆κ ≈ 10.5ω(t f ).
We can also observe how the increased likelihood of stochastic dissipation of the system in the PME impacts both the average and fluctuations in work done.tions of the coupling and inverse temperature, and compared to the corresponding predictions from the standard weak-coupling model.Since there is no change in the free-energy of the system-bath from the beginning to end of the protocol, any work done is dissipative.We see that for this range of coupling strengths and temperatures, the weak-coupling theory typically underestimates the average work and its variance, with this discrepancy only becoming negligible at small coupling and low temperature.
To summarise our findings, we see clear differences between predictions of the WCME and PME at the level of full work statistics.The stronger coupling manifests in the work distribution by shifting the positions of various peaks associated with both adiabatic and non-adiabatic transitions due to a renormalisation of the system Hamiltonian.Furthermore, higher coupling tends to drive the system further from equilibrium due to the increased level of non-adiabaticity in the PME.This results in an increased likelihood of stochastic work done on the system and leads to significant increases in the average dissipation and work fluctuations when compared to predictions from the weak coupling approach.

Conclusion-
We have studied the full-counting statistics of work done in the time-dependent spin-boson model beyond the weak coupling regime and have derived a generalised master equation for the corresponding WCO (5).This has been achieved by exploiting the unitary invariance of the work characteristic function, from which it was possible to apply the polaron transformation to rotate the WCO into a new frame where we could revert to standard approximations from weak-coupling theory.This resulted in the master equation (12) that can be solved to give the full work probability distribution in parameter regimes where the standard Markovian model breaks down, such as strong coupling.We illustrated our general approach by solving the dissipative Landau-Zener model, where the polaron framework indicated significant increases in both average dissipation and work fluctuations in contrast to predictions from the weak coupling master equation.This serves to highlight the limitations of the usual weak-coupling approximation, as even modest increases in the coupling can have a noticeable impact on the stochastic thermodynamics of the system.
Our results lay the groundwork for further investigations into how correlations, non-Markovianity [26] and other strong coupling effects can influence the full statistics of work in driven open quantum systems.The polaron framework could be applied to study higher order work fluctuations in thermal machines beyond weak coupling such as quantum Otto cycles [52] and Carnot engines [9].It could also provide a route towards free energy estimation in strongly coupled systems via the fluctuation theorem (15).In principle, one could apply quantum jump unravelling techniques [53,54] with respect to the polaron renormalised Hamiltonian to record the work statistics, which in turn could be used to establish estimates of the mean force free energy change using standard methods [55,56].
We conclude by noting a number of ways in which our approach could be generalised.Firstly, the polaron transformation that we have used is restricted to reservoirs with super-Ohmic spectral densities.However, this could be straightforwardly generalised using variational polaron theory [57], which allows treatment of both Ohmic and sub-Ohmic spectral densities, as has been used recently to study heat transport in non-equilibrium steady states [58].Furthermore, the unitary invariance of the work characteristic function allows for other transformations such as the reaction coordinate mapping [59].This approach could provide a means for computing work statistics in other strongly coupled systems, such as linearly coupled oscillators and fermionic environments. Acknowledgements ∞ 0 ds e i(s−η)ωnm(t) B ab (s − η) .

(A45)
The integrals over the bath correlation function may be rewritten as Fourier transforms (denoted by operation F(.)).
For example, the last integral of (A45) can be expressed as where we have defined the counting-field shifted non-Hermitian part of the spectral density matrix: In the second equality we used the property of the correlation function B ab (τ ) = B ba (−τ ) * .Using the same procedure we find where we have defined S a,± (t) ≡ S aa (∓ω(t)) and S a,0 ≡ S aa (0).