The quantum-classical correspondence principle for work distributions

For closed quantum systems driven away from equilibrium, work is often defined in terms of projective measurements of initial and final energies. This definition leads to statistical distributions of work that satisfy nonequilibrium work and fluctuation relations. While this two-point measurement definition of quantum work can be justified heuristically by appeal to the first law of thermodynamics, its relationship to the classical definition of work has not been carefully examined. In this paper we employ semiclassical methods, combined with numerical simulations of a driven quartic oscillator, to study the correspondence between classical and quantal definitions of work in systems with one degree of freedom. We find that a semiclassical work distribution, built from classical trajectories that connect the initial and final energies, provides an excellent approximation to the quantum work distribution when the trajectories are assigned suitable phases and are allowed to interfere. Neglecting the interferences between trajectories reduces the distribution to that of the corresponding classical process. Hence, in the semiclassical limit, the quantum work distribution converges to the classical distribution, decorated by a quantum interference pattern. We also derive the form of the quantum work distribution at the boundary between classically allowed and forbidden regions, where this distribution tunnels into the forbidden region. Our results clarify how the correspondence principle applies in the context of quantum and classical work distributions, and contribute to the understanding of work and nonequilibrium work relations in the quantum regime.

For closed quantum systems driven away from equilibrium, work is often defined in terms of projective measurements of initial and final energies. This definition leads to statistical distributions of work that satisfy nonequilibrium work and fluctuation relations. While this two-point measurement definition of quantum work can be justified heuristically by appeal to the first law of thermodynamics, its relationship to the classical definition of work has not been carefully examined. In this paper we employ semiclassical methods, combined with numerical simulations of a driven quartic oscillator, to study the correspondence between classical and quantal definitions of work in systems with one degree of freedom. We find that a semiclassical work distribution, built from classical trajectories that connect the initial and final energies, provides an excellent approximation to the quantum work distribution when the trajectories are assigned suitable phases and are allowed to interfere. Neglecting the interferences between trajectories reduces the distribution to that of the corresponding classical process. Hence, in the semiclassical limit, the quantum work distribution converges to the classical distribution, decorated by a quantum interference pattern. We also derive the form of the quantum work distribution at the boundary between classically allowed and forbidden regions, where this distribution tunnels into the forbidden region. Our results clarify how the correspondence principle applies in the context of quantum and classical work distributions, and contribute to the understanding of work and nonequilibrium work relations in the quantum regime.

I. INTRODUCTION
Work is a familiar concept in elementary mechanics and a central one in thermodynamics. In recent years, interest in the nonequilibrium thermodynamics of small systems [1][2][3][4] has motivated careful examinations of how to define quantum work . In this context one often considers a process in which a quantum system evolves under the Schrödinger equation as its Hamiltonian is varied in time -for instance a quantum particle in a piston undergoing compression or expansion [48]. It is typically assumed that the system is initialized in thermal equilibrium, and the question becomes: how do we appropriately define the work performed on the system during a single realization of this process?
One answer involves two projective measurements of the system's energy, at the start of the process (t = 0) and at the end (t = τ ). If the system Hamiltonian is varied fromĤ(0) =Ĥ A toĤ(τ ) =Ĥ B , then the measurement outcomes will be eigenvalues of these operators. The work performed during the process is then defined to be the difference between these two values, e.g.
if the measurements produce the m'th and n'th eigenvalues of the initial and final Hamiltonians. In this definition, work is inherently stochastic, with two sources of randomness: the statistical randomness associated with sampling an initial energy E A m from the canonical equilibrium distribution (Eq. 21), and the purely quantal randomness associated with the "collapse" of the final wavefunction |ψ τ upon making a projective measurement of the final energy (Eq. 22) [49].
In an analogous classical process, a system is prepared in thermal equilibrium, then it evolves under Hamilton's equations as the Hamiltonian function is varied from H A (z) to H B (z), where z denotes a point in the system's phase space. Since the system is thermally isolated, it is natural to define the work as the change in its internal energy: for a realization during which the system evolves from z 0 to z τ . As in the quantum case W is a stochastic variable, but here the randomness arises solely from the sampling of the initial microstate z 0 from an equilibrium distribution.
Much of the recent interest in quantum work has been stimulated by the discovery and experimental verification of classical nonequilibrium work relations [50][51][52][53][54][55], which have provided insights into the second law of thermodynamics, particularly as it applies to small systems where fluctuations are important [2]. For thermally isolated systems, the quantal counterparts of these relations follow directly from the definition of work given by Eq. 1 [6,7,9,15]. Moreover, experimental tests of quantum nonequilibrium work relations have recently been proposed [56][57][58][59] and implemented [60,61], providing direct verification of the validity of these results.
Despite the evident similarity between Eqs. 1 and 2, and despite the relevance of Eq. 1 to nonequilibrium work relations, the definition of quantum work given by Eq. 1 might seem ad hoc. Indeed, given the absence of a broadly agreed-upon "textbook" definition of quantum work, one might reasonably suspect that Eq. 1 has been introduced precisely because it leads to quantum nonequilibrium work relations. Furthermore, the wavefunction collapse that occurs upon measuring the final energy has no classical counterpart. Due to this collapse, the quantum work has no independent physical reality until the measurement is performed at t = τ . By contrast the classical work can be viewed as a well-defined function of time, namely the net change in the value of the time-dependent Hamiltonian along the trajectory z t . Our aim in this paper is to use the tools of semiclassical mechanics, together with numerical simulations, to investigate the relationship between Eqs. 1 and 2, and specifically to clarify how the correspondence principle applies in this context.
We restrict ourselves to systems with one degree of freedom, for which there exist explicit semiclassical approximations of energy eigenstates. We focus on the transition probability P Q (n|m) from the m'th eigenstate ofĤ A to the n'th eigenstate ofĤ B , and on its classical analogue, P C (n|m), defined in Sec. II. In Sec. III we investigate both quantities numerically for the example of a forced quartic oscillator. For high-lying initial energies P Q oscillates rapidly with n, whereas P C is a smooth function over a finite range n min < n < n max (Fig. 2). After integrating over the oscillations the two quantities are nearly identical (Fig. 3), which provides some justification for viewing Eq. 1 as the quantal counterpart of Eq. 2. In Sec. IV we derive P SC (n|m), a semiclassical approximation for P Q expressed as a sum over classical trajectories. Each trajectory in this sum carries a quantal phase, giving rise to coherent interferences between trajectories. When these interferences are neglected, P SC agrees with P C ; when they are included, P SC accurately captures the oscillations in P Q (Fig. 5). Thus the oscillations in P Q can be understood as a quantum interference pattern superimposed on a classical background. This picture breaks down around n min and n max , where P Q exhibits tails that tunnel into classically forbidden regions. In Sec. V we derive a semiclassical approximation, expressed in terms of the Airy function, that accurately describes these tails (Fig. 9).
The analysis in Secs. IV and V relies on theoretical tools that are familiar in the field of semiclassical mechanics. To the best of our knowledge, these tools have not previously been applied to study the relationship between classical and quantum work distributions, although similar analyses have been performed in the context of molecular scattering theory [62][63][64] and laser-pulsed atoms [65,66], among other examples. In particular, our calculations in Sec. IV B closely parallel those of Schwieters and Delos [66].

II. CLASSICAL AND QUANTUM TRANSITION PROBABILITIES
Here we introduce notation and specify the problem we plan to study. In Sec. II A we define a discretized classical work distribution, Eq. 17, that can be compared directly to the quantum work distribution, Eq. 20. These distributions involve classical and quantum transition probabilities, P C (n|m) and P Q (n|m), which will be the central objects of study throughout the rest of the paper.

A. Classical setup
Consider a system with one degree of freedom, described by a Hamiltonian where z = (q, p) denotes a point in phase space, and λ is an externally controlled parameter. We assume that the energy shells (level surfaces) of the Hamiltonian form simple, closed curves in phase space. We will consider the evolution of this system under the time-dependent Hamiltonian H(z; λ t ), where λ is varied from λ 0 = A to λ τ = B.
For compact notation, we define H A (z) ≡ H(z; λ 0 ) and H B (z) ≡ H(z; λ τ ). Suppose that prior to the start of the process, the system has come to equilibrium with a thermal reservoir at temperature T and the reservoir has been removed. Therefore at t = 0 the microstate of the system will be treated as a random sample from a canonical distribution corresponding to the Hamiltonian H A . From t = 0 to t = τ , the system is described by a trajectory z t evolving under Hamilton's equations, and the work performed on the system is given by the difference between its initial and final energies, as per Eq. 2: For an ensemble of realizations of this process, the distribution of values of work performed on the system is whereP C A (E 0 ) is the probability distribution of initial energies, sampled from equilibrium, andP C (E τ |E 0 ) is the conditional probability distribution to end with a final energy E τ , given an initial energy E 0 . Let us consider these factors separately.
P C A is simply the classical equilibrium energy distribution: where are the classical partition function and density of states, β = (k B T ) −1 , and h is Planck's constant. We have included the factors h −1 (which cancel in Eq. 6) for later convenience. Let us also define which is the phase space volume enclosed by the energy shell E of the Hamiltonian H, equivalently the integral of p dq around this energy shell. We then have The conditional probability distributionP C (E τ |E 0 ) is given by the expression where z τ (z 0 ) denotes the final conditions of a trajectory that evolves from initial conditions z 0 . It is useful to imagine an ensemble of initial conditions sampled microcanonically from the energy shell E 0 of the Hamiltonian H A (see Fig. 1(a)). The swarm of trajectories that evolves from these initial conditions defines a time-dependent, closed curve A t in phase space. At t = τ , the points of intersection between A τ and the energy shell E τ of the Hamiltonian H B (black dots in Fig. 1(b)) represent the final conditions of those trajectories that end with energy E τ , having begun with energy E 0 . If we consider two nearby energy shells E τ and E τ + dE, thenP C (E τ |E 0 ) dE τ is the fraction of trajectories whose final energies fall within this energy interval. Let us define the m'th energy interval of the Hamiltonian H A to be the range of energy values E satisfying mh ≤ Ω(E, A) < (m + 1)h (11) where m = 0, 1, 2 · · · . Furthermore, let us use a mid-point rule to assign a particular energy value E A m to each interval: Analogous definitions apply to the n'th energy interval of the Hamiltonian H B , and the corresponding energy E B n . From Eqs. 9 and 11 the width of the m'th energy interval of H A is For a smooth function of energy f (E), Eq. 13 leads to where the integral is taken over the m'th interval. Readers will recognize Eq. 12 as the semiclassical quantization condition p dq = [m + (1/2)]h, which provides an excellent approximation for the m'th eigenvalue of the quantum HamiltonianĤ A (Eq. 18 below) when m 1. For convenience, we will later use the notation E A m to denote this eigenvalue, without making a distinction between the exact eigenvalue and its semiclassical approximation.
The probability to obtain initial conditions z 0 within the m'th interval of H A , when sampling from equilibrium, is using Eq. 14. Similarly, the conditional probability to end in the n'th energy interval, given a representative initial energy in the m'th interval, is We will refer to this as the classical transition probability. The interpretation of P C (n|m) is straightforward. Imagine a swarm of trajectories evolving from initial conditions sampled from a microcanonical ensemble with energy E A m (see Fig. 1(a)). At the end of the process, t = τ , the fraction of these trajectories that fall into the energy window [E B n , E B n+1 ] is equal to P C (n|m) [66]. The classical work distribution can now be rewritten as a sum over initial and final energy intervals: Finally, we point out that the approximations appearing in Eqs. 13 -17 assume that the width of the m'th energy interval, δE m , is very small on a classical scale. Formally, this assumption can be justified by taking the semiclassical limit, → 0, with all classical quantities held fixed. In this limit the density of states increases, as does the quantum number at a given energy. In practice the semiclassical limit is often simply identified with the condition m 1.

B. Quantum setup
In the quantum version of this problem, H(z; λ) is replaced by the Hermitian operator A wavefunction ψ(q, t) = q|ψ t evolves under the Schrödinger equation asĤ(λ t ) is varied fromĤ(λ 0 ) =Ĥ A tô H(λ τ ) =Ĥ B . The eigenstates and eigenvalues ofĤ A are specified as follows: with similar notation forĤ B . Evolution in time is represented by the unitary operatorÛ t satisfyingĤ(λ t )Û t = i ∂Û t /∂t andÛ 0 =Î, whereÎ is the identity operator. Using Eq. 1 and assuming thermal equilibration at t = 0, the work distribution is [15] P Q A (m) is the probability of obtaining the m'th eigenstate ofĤ A when making the initial energy measurement: The quantum transition probability P Q (n|m) is the conditional probability to obtain the n'th eigenstate ofĤ B upon making the final measurement, given the m'th eigenstate ofĤ A at the initial measurement: The classical and quantum work distributions given by Eqs. 17 and 20 can now be compared directly. We note that hence P C A (m) ≈ P Q A (m) (see Eqs. 15 and 21). Thus what remains is to clarify the relationship between the classical and quantum transition probabilities, P C (n|m) and P Q (n|m). In the following section, we compare these transition probabilities in a model system for which both the classical and quantum dynamics are simulated numerically. We will find that P C and P Q are manifestly different (Fig. 2), but these differences mostly vanish after appropriate smoothing (Fig. 3). In Sections IV and V we develop a semiclassical theory to explain these features.

III. NUMERICAL CASE STUDY: THE FORCED QUARTIC OSCILLATOR
We begin with a model system, the forced quartic oscillator: We set M = 1/2 and = h/2π = 1, and we vary the work parameter at a constant rate, λ t = λ 0 +vt, from λ 0 = 1 = A to λ τ = 5 = B, taking v = 50 and therefore τ = 0.08. A number of groups have previously compared quantum and classical work distributions for a driven harmonic oscillator, for which analytical solutions are available. Specifically, Deffner et al [67,68] have studied an oscillator with a time-varying stiffness; Talkner et al [69] have obtained the work distribution for an oscillator driven by timedependent perturbations that are linear in q and p; Campisi [70] has studied changes in the Boltzmann entropy for forced quantum and classical oscillators; Ford et al [71] have investigated a harmonic oscillator with a cyclically driven stiffness; and Talkner et al have considered the sudden quench of a two-dimensional oscillator [26]. While the harmonic oscillator has the advantage of being analytically tractable, it is somewhat special in that its quantum dynamics can be reduced to its classical dynamics [72]. Here we have chosen the more "generic" quartic oscillator, which requires numerical simulations.
To evaluate P C (n|m) we simulated 10 4 Hamiltonian trajectories evolving under H(z; λ t ). Initial conditions were sampled microcanonically, and final conditions were binned according to the energy intervals of H B (z), as defined in Sec. II A. For our initial microcanonical ensemble, we chose m = 150, corresponding to E 0 = E A m = 1749. 23. The results are shown by the dashed line in Fig. 2. The transition probability is a smooth function of n from n min = 105 to n max = 215, and zero outside this range.
These sharp cutoffs reflect the fact that our initial conditions were sampled from a microcanonical distribution: n min and n max correspond to the minimal and maximal final energies that can be reached by trajectories launched with initial energies E 0 = E A m , as illustrated in Fig. 7 of Sec. V. The characteristic U-shape of the classical work distribution within the allowed region is related to the fact that the curves A and B shown in Fig. 1(b) become tangent to one another at E = E min and E = E max (again, see Fig. 7). For an exactly solvable model in which the same U-shape appears, see Fig. 1 of Ref. [70]. For the quantum system, we expand a wavefunction evolving underĤ(λ t ) =p 2 /2M + λ tq 4 as follows: HereĤ|φ n = E n |φ n and the c n 's are expansion coefficients. The time-dependent Schrödinger equation i ψ =Ĥψ then produces the set of coupled ordinary differential equations, We integrated these equations numerically, from initial conditions c n (0) = δ mn , to obtain In Fig. 2 we plot the quantum transition probability as a function of the final quantum number n (solid line). Although a correspondence between P Q (n|m) and P C (n|m) is visually evident, the quantum and classical cases differ in two distinct ways: (i) the quantum probability oscillates rapidly with n, and (ii) P Q (n|m) shows tails that "tunnel" into the classically forbidden regions n < n min and n > n max . Both features trace their origin to the wave nature of the quantum system. We also see in Fig. 2 that P Q (n|m) = 0 when n − m is an odd. This result follows from the fact that φ n q 4 φ k = 0 when the parities of n and k differ, hence only transitions between states of the same parity occur under Eq. 26. FIG. 3. Accumulated transition probabilities for the forced quartic oscillator. The jagged black curve represents the quantum case, n k=0 P Q (k|m), while the smooth red curve represents the classical case, n k=0 P C (k|m).
In Fig. 3 we plot the accumulated transition probabilities n k=0 P Q (k|m) and n k=0 P C (k|m), thereby smoothing out the rapid oscillations in the quantum transition probability. The close agreement observed in Fig. 3 suggests that Eq. 1 is indeed the appropriate quantum counterpart of Eq. 2, though distinctly non-classical features are visible in the interference and tunneling effects in Fig. 2.
In Sections IV and V we investigate these issues analytically. We will find that both the agreement seen in Fig. 3 and the quantum effects evident in Fig. 2 can be understood quantitatively, within a semiclassical interpretation in which quantum dynamics are approximated by classical trajectories bearing time-dependent phases.

IV. SEMICLASSICAL THEORY
In Sec. IV A we rewrite the classical transition probability P C (n|m) (Eq. 16) as a sum over trajectories that begin and end with energies E A m and E B n , respectively (Eq. 39). In Sec. IV B we use time-dependent WKB theory to derive a semiclassical approximation P SC (n|m) for the quantum transition probability P Q (n|m). Our main result, Eq. 49, is expressed as a sum over the same trajectories that contribute to P C (n|m), only now each of these trajectories carries a phase, resulting in interference effects between different trajectories. When these interferences are ignored, we recover the classical transition probability (Eq. 50). When they are included, the interferences give rise to the oscillations observed in the quantum transition probability. In particular, interferences between symmetry-related trajectories account for the fastest oscillations observed in Fig. 2, as we discuss in Sec. IV C.

A. Classical transition probabilities
We begin with the conditional probability distributionP C (E|E 0 ) (Eq. 10), where we have dropped the subscript τ from the final energy, for convenience. As discussed in Sec. II, we imagine a swarm of trajectories evolving in time under H(z; λ t ), with initial conditions sampled from the microcanonical phase space distribution where E 0 is a parameter of the distribution. The quantityP C (E|E 0 ) dE is the fraction of these trajectories that end with a final energy in the infintesimal range (E, E + dE). As in Fig. 1, let A t denote the time-dependent curve that evolves from the initial energy shell H A = E 0 , and let B denote an energy shell of the final Hamiltonian, H B = E. The evolving curve A t can be described by a multivalued, time-dependent momentum field p At b (q), where the index b labels the branches of A t at the coordinate value q. E.g. in Fig. 1(b) the momentum field for A τ has two branches at q = q 1 . Similarly, the multivalued momentum field p B b (q) describes the fixed energy shell H B = E, whose two branches are For convenience, we will henceforth use the notation A (without a subscript) to indicate the surface A τ . The points of intersection between A and B, highlighted by dots in Fig. 1(b), represent the trajectories that contribute toP C (E|E 0 ). Our swarm of trajectories at time t is described not only by the evolving surface A t , but also by a probability density on that surface. If we project this density onto the q-axis, the projected density ρ At (q) is a sum over the branches of the surface A t : The fixed microcanonical ensemble corresponding to H B = E can similarly be projected onto a density with ∂H B /∂p = p/M evaluated at p = p B b (q). Let l = 1, 2, · · · K be an index labeling the intersection points of A and B (e.g. K = 4 in Fig. 1), and let b l denote the branch corresponding to the l th intersection point: With some abuse of notation, we have written b l rather than b A l and b B l , though in general the l'th intersection point may occur at differently indexed branches of the two surfaces.
To evaluateP C (E|E 0 ), we begin with the contributionP C l (E|E 0 ) from a single intersection point, l. Fig. 4 depicts this intersection, together with an intersection point involving an infinitesimally displaced energy shell H B = E + δE. The highlighted line segment connecting these two points represents a set of trajectories with final energies between E and E + δE, hence Since we focus here on the contribution from a single intersection between A and B, we suppress the subscripts b l indicating the branches of these surfaces. By construction, the ratio δE/δq is the rate of change of H B with respect to q along the curve A: Since p B (q) specifies a surface of constant H B , we have Combining Eqs. (31), (33), (34), and (36) we get Summing over all intersection points gives us If we now set E 0 = E A m and E = E B n , then Eq. 16 gives us This is the main result of this subsection. We will use the term classically allowed region to denote the range E min ≤ E ≤ E max within which the function P C (E|E 0 ) does not vanish, representing all final energies that can be attained by trajectories with initial energy E 0 . In the coarse-grained function P C (n|m), the classically allowed region is bracketed by the values n min and n max , denoting the energy intervals of H B containing E min and E max . As illustrated in Fig. 7 below, the values E min and E max correspond, respectively, to the largest energy shell of H B that fits entirely within the closed curve A, and the smallest energy shell of H B that surrounds the entire curve A. At these two energies, the curves A and B become tangent to one another, ∂p A b l /∂q = ∂p B b l /∂q, leading to divergences in Eq. (38), which are reflected in the two sharp peaks in P C (n|m) in Fig. 2. Similar divergences appear in a more familiar context, namely the density ρ B (q) that describes the projection of a fixed microcanonical ensemble H B = E onto the coordinate axis (see Eq. (31)). In this case the divergences occur at classical turning points at which the velocity ∂H B /∂p vanishes.
In the example introduced in Sec. III, the reflection symmetry of the quartic potential implies that trajectories come in symmetry-related pairs: if (q t , p t ) is a solution of Hamilton's equations under the time-dependent Hamiltonian H(q, p; λ t ), then so is (−q t , −p t ). This means that for every intersection point (q l , p l ) between A and B, there will be another intersection point at (−q l , −p l ), and the two will contribute equally to P C (n|m).

B. Semiclassical transition probabilities
Let us now evaluateP Q (n|m) in the semiclassical limit. Appendix A provides a brief introduction to time-dependent WKB theory, and for further details we refer to the reviews by Delos [73] and Littlejohn [74]. Using Eq. A5, the wavefunction ψ(q, t) = q|Û t |φ A m and the n'th eigenstate ofĤ B can be written as where the actions S At b (q) and S B b (q) generate the momentum fields which appears as Eq. 37b in Ref. [66]. Using the stationary phase approximation to evaluate the integral, we obtain a sum of contributions from points q l satisfying the condition which is equivalent to Hence φ B n |ψ τ reduces to a sum of contributions arising from the K intersection points of the surfaces A and B, representing trajectories with initial and final energies E A m and E B n , just as in the classical case (Eq. 32). To determine the contribution from the l'th intersection point, we perform a quadratic expansion around q = q l : Performing the stationary phase integral then gives us: where ρ A,B b l are evaluated at q = q l ; ∆µ l = µ A b l − µ B b l ; and σ l = sign(κ l ) = ±1. In Ref. [66], where similar calculations are performed, a quantity equivalent to our πσ l /4 is identified as a Maslov-like phase (Ref. [66], p. 1037).
Eq. (47) expresses φ B n |ψ τ as a sum over trajectories, each contributing an amplitude and a phase: The semiclassical transition probability P SC (n|m) ≈ φ B n ψ τ 2 is then given by The general form of Eq. 49 is not surprising. The path integral formulation of quantum mechanics tells us that solutions of the time-dependent Schrödinger equation can be represented in terms of sums over classical paths decorated by phases exp(iS/ ), and when → 0 the dominant contributions to the sums come from those paths that are solutions of the classical equations of motion [75]. We expect the semiclassical transition probability function P SC (n|m) to provide a good approximation to the quantum transition probability P Q (n|m), for large quantum numbers. For the forced quartic oscillator, Fig. 5 compares these functions, taking m = 150. P Q was computed as described in Sec. III, and P SC was evaluated directly from classical simulations of trajectories evolving from an initial microcanonical ensemble. (The evolution of the multivalued action S At was obtained by integrating p dq − H dt along each trajectory, see e.g. Sec. 3.5 of Ref. [74].) We observe that the agreement is excellent, except in the vicinity of the boundaries between the classically allowed and forbidden regions. We will examine these boundaries in more detail in Sec. V. Now let us simplify Eq. 49 by simply ignoring the cross terms in the double sum. This is the so-called diagonal approximation [76][77][78], and it leads to which is identical to the expression for the classical transition probability, P C (n|m) (Eq. 39). The agreement between P Q and P SC = | a l e iθ l | 2 seen in Fig. 5 suggests that the rapid oscillations of the quantum probability distribution can be understood in terms of phase interference between different trajectories. When these interferences are ignored, we recover the classical probability distribution, as we see in Eq. 50. Hence P SC acts as a bridge to between the quantum and classical transition probabilities, and thus between the quantum and classical work distributions (Eqs. 20, 17). The first approximation in Eq. 51 involves time-dependent WKB theory together with a stationary phase evaluation of integrals, and the second is the diagonal approximation in which interferences are ignored.

C. Interferences and symmetries
Let us now examine the interferences in more detail. Recall from Sec. IV A that the trajectories contributing to Eq. 49 come in symmetry-related pairs. This allows us to write: a l e iθ l + a l+J e iθ l+J (52) where J = K/2 is an integer, and the l'th and (l + J)'th trajectories are related by symmetry: Here we have indexed the points (q l , p l ) from l = 1 to K, according to the order in which they appear as we proceed clockwise around the energy shell B, first rightward along the upper branch and then leftward along the lower branch. Eq. 53 implies Using the fact that the points l and l + J are located directly opposite one another on the closed curves A and B (Eq. 53), we obtain and hence (In Eq. 55 we have invoked the semiclassical quantization condition, Eq. 12, together with Liouville's theorem, which guarantees that At p dq remains constant with time. In Eq. 56 we have used the fact that the Maslov index is incremented by +2 as one proceeds around the closed curve A or B, see e.g. Fig. 1  a l e iθ l 1 + e i(m−n)π .
When m − n is even, the symmetry-related trajectories interfere constructively, but when m and n have opposite parities they interfere destructively and we get P SC (n|m) = 0. This provides a semiclassical explanation for the observation, noted in Sec. III, that the quantum transition probability vanishes when m − n is odd: the most rapid oscillations in P Q (n|m) arise from interference between symmetry-related trajectories. The same interpretation was given by Miller to explain angular momentum selection rules in the scattering matrix describing an atom impinging on a homonuclear diatomic molecule (see Sec. III C of Ref. [64]). The interference effects in P SC , due to cross terms in Eq. 49, include contributions not only from pairs of trajectories that are related by symmetry (Eq. 53), but also from those that are not related by symmetry: where s denotes a sum of symmetry-related pairs of trajectories, (q k , p k ) = −(q l , p l ), and n is a sum over nonsymmetry-related pairs, (q k , p k ) = −(q l , p l ). To illustrate these contributions separately, we constructed the quantity P SC sym by omitting the sum n on the right side of Eq. 59, and we similarly constructed P SC nosym by omitting s . The results are shown in Fig. 6. P SC sym displays a regular oscillation, as the symmetry-related trajectories are either exactly in phase or exactly out of phase, according to the parity of m − n. The oscillations in P SC nosym are less regular, as the phases θ l − θ k are not necessarily integer multiples of π. When all cross terms are included, both sets of oscillations combine to give the pattern in Fig. 5. For instance, the gap observed near n = 200 in Fig. 5 reflects the corresponding minimum in P SC nosym seen in Fig. 6(b).

V. AIRY TAILS
In Fig. 5 we see that that our semiclassical approximation fails around the boundaries of the classical allowed region at n min and n max . Although the semiclassical transition probability derived in Sec. IV vanishes identically outside the classically allowed region (since only classical trajectories contribute to P SC ), the quantum transition probability tunnels into the classically forbidden region, as seen in Fig. 5. In this section we construct a semiclassical approximation that describes this behavior. Our central results are given by Eqs. 64 and 69, and illustrated in Fig. 9.
Similar tunneling behavior is observed in the wavefunctions of energy eigenstates, which exhibit tails that reach into classically forbidden regions [79], as well as in semiclassical treatments of the S matrix in molecular collisions [63]; in both cases the tails are approximated by Airy functions. The expressions that we obtain in this section will also be expressed in terms of Airy functions.
The failure of P SC (n|m) in the boundary regions originates in the fact that the curves A and B become tangent to one another at E = E min and E = E max . In the example we have studied numerically the curves are tangent simultaneously at two points, due to the symmetry of the quartic potential, as shown in Fig. 7(b). However, it is useful first to discuss the generic case, where there exists only a single point of tangency, depicted in Fig. 7(a). For specificity we will discuss the upper boundary of the classically allowed region at E max , but similar comments apply to the lower boundary at E min .
Let us consider the function that appears inside the exponent in Eq. 42, and its derivative We have suppressed the subscripts b and b , and we have explicitly indicated that S B and p B depend on the final energy shell E. In Sec. IV we evaluated the integral dq exp(i∆S/ ) by performing a quadratic expansion of ∆S around a single intersection point q l between the curves A and B, where ∆p vanishes (Fig. 4). As the energy E approaches E max from below, two points of intersection between A and B coalesce at a point (q c , p c ), indicated in Fig. 7(a). Fig. 8 illustrates the behavior of ∆p(q, E) and ∆S(q, E) in the vicinity of q c and E max . In particular, ∆S has two stationary points when E < E max and none when E > E max . The stationary phase integration performed in Sec. IV (Eqs. 43 -47) implicitly assumed well-isolated stationary points, but this assumption breaks down near E = E max . In this situation we evaluate the integral dq exp(i∆S/ ) by expanding ∆S to cubic rather than quadratic order. Leaving the details of the calculation to Appendix B, here we present the result: where and Ai(·) denotes the Airy function. Combining this result with Eq. 42, we obtain the following expression, which is valid when E B n ≈ E max : where ρ A,B c = ρ A,B (q c ), and the subscript 1 emphasizes that we have assumed only a single point of tangency between A and B. At the lower boundary of the classically allowed region, we obtain the same result but with E max replaced by E min in Eqs. 63 and 64.
The Airy function decays rapidly for positive values of its argument, and oscillates for negative values: for real ζ 1. One can establish by inspection that ν > 0 > k when E > E max and ν, k > 0 when E < E min , therefore the argument of the Airy function in P SC,tail 1 is positive outside the classically allowed region. As a result the transition probability decays monotonically, penetrating the classically forbidden region over a characteristic skin depth |k 1/3 2/3 ν −1 |.
Let us now consider the situation corresponding to the example of Sec. III, in which there are simultaneously two points of tangency between the curves A and B at E = E max , and these two points are related by symmetry ( Fig. 7(b)). Using the subscripts c and d to distinguish these two points, we have: Summing the contributions from these two points gives us where ∆S c,d are evaluated at E = E max . Following arguments similar to those leading to Eq. 58, we obtain Eqs. 67 and 68 lead to a result identical to Eq. 64, apart from a factor that captures the interference between the two parity-related trajectories: with the subscript indicating two points of tangency between A and B. Fig. 9 compares the quantum transition probability P Q with our semiclassical approximations given by P SC (Eq. 49) and P SC,tail 2 (Eq. 69), in the central and boundary regions, respectively. The agreement is excellent, illustrating that the entire quantum work distribution can be understood in terms of contributions from individual classical trajectories, and the interferences between them.

VI. DISCUSSION AND CONCLUSION
The correspondence principle broadly states that classical mechanics is recovered from quantum mechanics in the limit of large quantum numbers. In textbook discussions this principle is often illustrated in terms of static properties, typically by comparing the coordinate space probability distribution of a harmonic oscillator eigenstate with the classical, microcanonical distribution at the same energy [79,80]. Here, by contrast, we have considered the correspondence principle in a dynamic setting, comparing classical and quantum transition probabilities under a time-dependent Hamiltonian. In this setting, time-dependent WKB theory provides a set of tools for using classical trajectories to construct approximate solutions of quantum dynamics.
We have applied these tools to study the relationship between the definitions of quantum and classical work given by Eqs. 1 and 2. Focusing on the transition probabilities P Q (n|m) and P C (n|m) for systems with one degree of freedom, we have obtained three main conclusions, and have illustrated them with simulations of a driven quartic oscillator: • In the classically allowed region, the quantum transition probability can be approximated in terms of interfering classical trajectories (Eq. 49, Fig. 5).
• When the interferences between these trajectories are ignored, the classical transition probability is obtained (Eq. 50, Fig. 3).
• The tunneling of the quantum transition probability into the classically forbidden region is accurately described by a decaying Airy function, whose arguments are expressed in terms of classical quantities (Eqs. 64, 69, Fig. 9).
The second conclusion, and Fig. 3 in particular, clarifies the sense in which classical mechanics is "recovered" in the limit of large quantum numbers (in this context), whereas the first and third conclusions describe inherently quantal effects -interference, tunneling -that can nevertheless be understood and approximated in terms of classical trajectories. In view of similar analyses in atomic and molecular contexts (see e.g. Refs. [62][63][64][65][66]) these conclusions are not surprising. However, given recent interest in quantum nonequilibrium work relations, our results provide a timely investigation of the correspondence principle as it applies to the work performed on a system driven away from equilibrium. As we have shown, semiclassical mechanics provides a bridge between classical and quantal work distributions, and our conclusions provide some justification for the two-point measurement definition of quantum work (Eq. 1). It will be interesting to study whether our results generalize to systems with more than one degree of freedom. For N -particle systems with N 1, we generically expect that the classical transition probability P C (n|m) will be bell-shaped rather than U-shaped, with n min and n max found deep in the tails of P C (n|m). In this situation the quantum tunneling into the forbidden regions will be negligible, since both P C (n|m) and P Q (n|m) will be very small at n min and n max . An obstacle to studying N -particle systems analytically is the lack of semiclassical representations of energy eigenstates, analogous to Eq. 40b, for generic systems with multiple degrees of freedom. Progress might be made by studying two limiting cases: integrable and fully chaotic systems [81]. For systems whose classical dynamics are integrable, semiclassical expressions for energy eigenstates can be constructed using classical action-angle variables (see e.g. Eq. 11 of Ref. [82] or Eq. 3.5 of Ref. [83] or Eq. 37 of Ref. [84]). At the other extreme, for fully chaotic systems, Berry's conjecture [85] suggests that energy eigenstates can be treated as Gaussian random functions of the configuration, q. In either case, the expression for the energy eigenstate might be combined with time-dependent WKB theory as a first step toward generalizing Eq. 42 and subsequent results to multi-dimensional systems. 2010363). SR gratefully acknowledges the support of the Israel Science Foundation (grant 924/11). This work was partially supported by the COST Action MP1209.

Appendix A: Brief review of time-dependent WKB theory
In time-dependent WKB theory, approximate solutions of the time-dependent Schrödinger equation are expressed in terms of classical constructions in phase space. In the brief summary that follows, we restrict ourselves to a single degree of freedom, and a Hamiltonian of the form H(q, p, t) = p 2 /2M + V (q, t).
A wavefunction ψ(q, t) is said to be in WKB form when written as the product of a slowly varying amplitude and a rapidly oscillating phase, or else as a sum of such terms, as in Eq. (A5) below. Locally, Eq. (A1) describes a wave train ψ ∝ exp(ikq), with wave number k = −1 ∂S/∂q. This implies a local momentum and the Born rule implies a probability density These functions inherit their time-dependence from ψ(q, t). A standard calculation, in which Eq. (A1) is substituted into the Schrödinger equation and is treated as a small parameter [74], produces the equations of motion where terms of order 2 have been neglected. Eq. A4a is the Hamilton-Jacobi equation. Eq. A4b is the continuity equation for evolution under a velocity field p(q, t)/M . To interpret these equations, let the function p(q, 0) describe a curve A 0 in phase space, at an initial time t = 0. If the wavefunction ψ is represented as a sum of terms of the form given by Eq. A1, then p(q, 0) is multivalued, and the curve A 0 has multiple branches, but for the moment we focus on a single branch. A 0 is a Lagrangian manifold, determined by its generating function, S(q, 0), via Eq. A2. Now imagine that this curve carries a probability density whose projection onto the coordinate axis is ρ(q, 0). It is convenient to picture a swarm of particles, sprinkled along A 0 so that the fraction of particles between q and q + dq is given by ρ(q, 0) dq. From these initial conditions, each particle in the swarm evolves under Hamilton's equations, giving rise to a time-dependent curve A t and density ρ(q, t). The Hamilton-Jacobi equation governs the evolution of A t , through its generating function S(q, t), and the continuity equation governs the evolution of the projected density ρ(q, t).
With these considerations in mind, approximate solutions of the time-dependent Schrödinger equation are written in the form The sum is taken over the branches of a Lagrangian manifold A t whose generating function S(q, t) satisfies the Hamilton-Jacobi equation (Eq. A4a), and the projected probability density for each branch satisfies the continuity equation (Eq. A4b). The quantities µ b are Maslov indices. These determine the relative quantum phases of the various branches of A t , and they differ from one another by integer values. (By convention, the overall phase of ψ is adjusted so that the µ b 's themselves are integers.) See Ref. [74] for a more detailed discussion of Maslov indices, as well as a careful treatment of caustic points, where two branches of A t meet and the manifold becomes "vertical" in the (q, p)-plane.
At the semiclassical level of approximation, Eq. (A5) connects the evolution of a quantum wavefunction to that of a swarm of classical particles "surfing" on a Lagrangian manifold. This perspective provides a natural starting point for a comparison between quantum and classical work distributions.
As an important example of a WKB wavefunction, consider a manifold A 0 that is a single energy shell E of a fixed Hamiltonian H(q, p), and consider a microcanonical distribution of initial conditions on this manifold. Under time evolution, neither the manifold nor the probability distribution changes, as each trajectory merely goes round the energy shell. The solution of Eq. A4a in this case is hence by Eq. A5 the wavefunction simply acquires a phase e −iEt/ . Thus in the semiclassical limit, energy eigenstates correspond to microcanonical ensembles. Indeed, Eq. A5 in this case leads to the WKB approximation for energy eigenstates that is familiar from undergraduate textbooks on quantum mechanics, where it is typically derived in a different manner. Moreover, a proper treatment of the Maslov indices leads to the quantization condition p(q) dq = [m + (1/2)]h.
with ∆p = p A − p B = ∆S (Eq. 61). Note that since A and B are tangent at q c when E = E c . The integral I(E) is now given by Eq. (B5), with β 0 , β 1 and β 3 obtained from α 0 , α 1 , α 2 and α 3 via Eq. (B4). This result simplifies if we consider how the β's behave in the vicinity of E = E c as → 0. Let us define Now look at the argument of the Airy function in Eq. (B5), and consider a fixed range of ζ values, ζ − < 0 < ζ + , chosen so that the asymptotic approximations for the Airy function (Eq. 65) are accurate at both ζ − and ζ + . Noting that ζ(E c , ) = 0 and expanding ζ(E, ) to first order in (E − E c ), we see that as → 0, the fixed range [ζ − , ζ + ] translates to a range of energies whose width scales like 2/3 . Thus we will be interested in energies that have the following scaling relation: where = E − E c . To leading order in , and expressing quantities in terms of rather than E, we have Now recall that β 1 = α 1 − α 2 2 /2α 3 . In the vicinity of = 0 (where α 1 = α 2 = 0), we have α 1 ∼ ∼ 2/3 and α 2 2 /2α 3 ∼ 2 ∼ 4/3 , hence the latter term can be ignored. We thus write Combining results and discarding terms that vanish as → 0, Eq. (B5) finally gives us which is equivalent to Eq. 62.