Theory of Quantum Work in Metallic Grains

We generalize Anderson's orthogonality determinant formula to describe the statistics of work performed on generic disordered, non-interacting fermionic nanograins during quantum quenches. The energy absorbed increases linearly with time, while its variance exhibits a superdiffusive behavior due to Pauli's exclusion principle. The probability of adiabatic evolution decays as a stretched exponential. In slowly driven systems, work statistics exhibits universal features, and can be understood in terms of fermion diffusion in energy space, generated by Landau-Zener transitions. This diffusion is very well captured by a Markovian symmetrical exclusion process, with the diffusion constant identified as the energy absorption rate. The energy absorption rate shows an anomalous frequency dependence at small energies, reflecting the symmetry class of the underlying Hamiltonian. Our predictions can be experimentally verified by calorimetric measurements performed on nanoscale circuits.


I. INTRODUCTION
While energy transfer, heat, and work are fundamental concepts in thermodynamics and statistical physics, it is far from trivial to extend their concepts to generic non-equilibrium quantum systems [1,2], where energy becomes a fluctuating statistical quantity even in pure quantum states, and energy transfer can only be understood in terms of a precise measurement protocol. Recent experimental developments allow, however, the investigation of these notions in quantum systems ranging from individual molecules subject to mechanical forces [3][4][5], through nuclear spins in a magnetic field [6], to mesoscopic grains [7,8], and allow even to extract the full distribution of the energy transferred.
In this work, we study the properties and universal aspects of quantum work produced in course of time dependent quantum quenches in generic, non-interacting, but disordered fermionic many-body nanosystems. These systems provide an ideal platform to study quantum thermodynamics and the interplay of quantum time evolution, disorder, and quantum statistics.
Generic noninteracting nanosystems obey random matrix statistics [26,27], and can be described in terms of a HamiltonianĤ with theâ i denoting fermionic annihilation operators, and H(t) an N × N random matrix. We assume that the system is in its M -particle ground state at time t = 0, and then work is performed by changing external gate voltages or by applying time dependent magnetic fields [28] (see Fig. 1). We describe this situation by performing a quench (ramp) at a constant pace, v ≡ Tr(d t H 2 ) 1/2 , related to the frequency ω of external parameters, and investigate the distribution of the internal energy injected, with the two averages referring to quantum and random matrix averages, respectively. Though work and heat cannot be quite disentangled in course of the quench process, throughout this paper we follow the standard practice [2] and refer to this internal energy change -directly related to calorimetric measurements in driven nanosystems [29,30] -as work. Following Ref. [31], we apply a quench protocol, H(t) = H 1 cos λ(t)+ H 2 sin λ(t), withλ set to constant, and H 1,2 independent N ×N matrices, drawn from a Gaussian random matrix ensemble [26,32], P (H) ∼ e −β N 4 TrH 2 . Here we focus on Gaussian orthogonal (β = 1) and Gaussian unitary (β = 2) ensembles, corresponding to integer spin time reversal invariant systems, and systems with broken time reversal symmetry, respectively.
We first construct a determinant formula for the generating function of work in non-interacting fermionic systems. The determinant formula presented can be considered as the dynamical analogue of Anderson's determinant formula for the orthogonality catastrophe [11]. However, while Anderson restricts himself to the probability of staying in the ground state after a local sudden quench (placing a scatterer at the origin), the determinant formula used here describes the response to any dynamical quench in a non-interacting fermionic many-body system, and describes all possible transitions and the corresponding overlaps.
We find that energy is absorbed by the system via hard core particle diffusion in energy space, and that the statistical properties of work depend crucially on the speed v of the quench as well as on underlying symmetries. For slow, almost adiabatic changes, in particular, P t (W ) displays a universal structure, with surprisingly large, superdiffusive work fluctuations, δW 2 ∼ W 3/2 . Energy absorption at small frequencies, ω ∼ v, reflects the symmetry (universality class) of the Hamiltonian and is predicted to scale as, W ∼ ω 1+β/2 , with the random matrix parameter β = 1, 2 and 4 corresponding to orthogonal, unitary, and symplectic symmetries.
The paper is organized as follows. We discuss the time dependent many-body wave function and present a determinant formula for the characteristic function of work in Sec. II. We study the full distribution of work in Sec. III. In Sec. IV we turn to the average work and show that it can be understood in terms of a diffusion in energy space. To shed more light on the dynamics, in Sec. V we show that the most important features of the work statistics can be captured by a classical symmetric exclusion process, and in Sec. VI we also present a simple mean field description allowing the derivation of analyti-cal results. We summarize our main findings in Sec. VII. Technical details of the calculations, as well as additional numerical results for the Gaussian unitary ensemble are relegated to appendices.

II. MANY-BODY WAVE FUNCTION
To compute the distribution (2), we first need to determine the many-body wave function of our electron system after the quench, |Ψ(t) . To do that, we exploit the fact that our Hamiltonian is non-interacting, and -the initial state of the system being a Slater determinant -the state |Ψ(t) can also be expressed as a Slater determinant at any time in terms of the time-evolved single particle wave functions. Alternatively, |Ψ(t) can be written as with the operatorsĉ † m,t creating single particles in the states ϕ m (t);ĉ † m,t |0 ≡ |ϕ m (t) . Here the vectors ϕ m (t) satisfy the single particle Schrödinger equation with the boundary conditions ϕ m i (0) = δ m i . To solve Eq. (4), we use the adiabatic approach. We introduce the instantaneous eigenvectors of H(t), η m t , satisfying H(t) η m t = ε m (t)η m t , and corresponding fermionic creation operators, |η m t ≡b † m,t |0 . We then expand the ϕ m (t)'s in the instantaneous basis as and determine the expansion coefficients by solving the corresponding equation of motion, with A kl = −i η k t · ∂ t η l t the Berry connection [33] (for details of the numerical calculation see Appendix B).
Sinceĉ † m,t = k α m k (t)b † k,t , knowledge of the coefficients α m k (t) allows us to express the many-body state |Ψ(t) in the instantaneous basis. Observing furthermore that the many-body Hamiltonian assumes a particularly simple form in the instantaneous basis,Ĥ(t) = m ε m (t)b † m,tb m,t , we can rewrite the expectation value Ψ(t)| e iuĤ(t) |Ψ(t) as We can now evaluate the expectation values using Wick's theorem and eliminate the summation over the labels Transferring furthermore the permutation operators from the lower index of the α m k 's to their upper index and reordering sums and products yields with P running over all permutations of the M occupied states, and the M × M matrix g t (u) incorporating all information on the overlap between initial and final single particle states, This leads immediately to the dererminant formula, Eq. (10) is one of the important results of this work: it establishes a connection between the work statistics of the many-body system and the time evolution of individual single particle states [34,35]. In the following, we shall use this formula to study the properties of quantum quenches in generic fermionic nano-systems described by random matrix theory.

III. QUANTUM STATISTICS OF WORK
To determine the full distribution P t (W ), we utilized Eqs. (6) and (10). We generated random matrices H 1 and H 2 , determined the expansion coefficients α m k (t) and the determinant det g t (u) numerically, averaged over the random matrix ensemble, and finally determined P t (W ) by taking the Fourier transform of Eq. (10). In our numerical simulations, we focus on the half-filled case, M = N/2, though the results are independent of this assumption and carry over to any filling M/N . Our results are summarized in Figs. 2 and 3.
The function P t (W ) can be disentangled into an adiabatic (ground state) and a regular part as These functions depend not only on time, but also on the velocity v, by which we drag the Hamiltonian through the random matrix manifold, on the position of the Fermi level, and on the symmetry class, too. However, once expressed in terms of appropriate dimensionless quantities, ∼ W ≡ W/∆,t ≡ t ∆, andṽ ≡ v/∆, measured in units of the average single particle level spacing ∆ at the Fermi energy, they become universal functions in the limit N → ∞, P reg → P reg ( ∼ W ;t) and P GS → P GS (t), displayed in Figs. 2 and 3 for the orthogonal ensemble.
They depend implicitly onṽ but, supposedly, they are independent of all microscopic details, and depend just on the symmetry of the underlying Hamiltonian (see also Appendix C for data for the Gaussian unitary ensemble). Remarkably, these functions simplify even further in the slow quench limit,ṽ < ∼ 1, and become functions only of the average work, ∼ W , rather than time and velocity. This is demonstrated for the probability of adiabatic transitions in Fig. 3 within the orthogonal ensemble.
The functions P reg ( ∼ W ;t) show a rather slow evolution towards a Gaussian distribution, as we increase ∼ W . As closer numerical investigation reveals, the work fluctuations follow a superdiffusive scaling, δ The behavior of P GS (t) is also somewhat unexpected: the probability of an adiabatic transition is found to decay as a stretched exponential, P GS (t) ∼t 1/4 e −C √t . As we discuss and demonstrate below, both are particular features of a symmetrical exclusion process which governs the dynamics in energy space.

IV. ENERGY SPACE DIFFUSION AND
AVERAGE WORK Fig. 4 shows the average of the occupation of each level f k (t) ≡ n k,t RM . [37] The occupation profile exhibits a clear diffusive character, and is very precisely described by a diffusively broadened Fermi-sea, where ∆k ≡ k − M is the distance of level k from the Fermi energy, and D denotes a dimensionless diffusion constant in energy space. The diffusive ∼ √ t broadening of the Fermi surface immediately implies a linear internal energy absorption, ∼ t, as clearly demonstrated in the inset of Fig. 4. The diffusion constant D can thus be interpreted as an overall dimensionless energy absorption rate. The constant D depends on the universality class of the system as well as on the velocity,ṽ. We can distinguish two distinct regimes: in the "fast limit",ṽ ≫ 1, electronhole transitions between remote levels dominate energy absorption, and we find D(ṽ) ∼ṽ 2 ∼ ω 2 , as expected in metals. However, in the "slow limit",ṽ < ∼ 1, nearest neighbor transitions mediated by Landau-Zener transitions dominate. The statistics of these latter have been studied thoroughly [31,38,42], and it has been shown that the parameters of Landau-Zener transitions display universal distributions (see also Appendix A for details). A simple calculation making use of this universal statistics then yields an energy absorption D(ṽ) ∼ṽ 1+β/2 ∼ |ω| 1+β/2 forṽ < ∼ 1, as indeed confirmed by our detailed numerical simulations (see Appendix A).

FIG. 2.
Time evolution of work statistics for the Gaussian orthogonal ensemble in the limit,ṽ < ∼ 1. Circles: quantum results for the regular part, Preg( ∼ W ;t), for dimensionless velocityṽ = 0.8 and N = 20 levels. In the three panels, we used different dimensionless quench timest = 4.88, 6.25, and 162.5 (from (a) to (c)), corresponding to dimensionless average work ∼ W = 1.5, 5, and 50, respectively [36]. Solid lines: SEP simulations. The distribution becomes more and more Gaussian-like as ∼ W increases.

V. MARKOVIAN SIMULATION AND SYMMETRICAL EXCLUSION PROCESS
In the slow limit,ṽ < ∼ 1, dominated by Landau-Zener transitions, most features of the work statistics can be understood in terms of a simple, classical model, the symmetrical exclusion process (SEP) in energy space. In this approach, we consider the occupation of each level as a classical statistical variable, taking values n k,t = 1 and 0, and think of Landau-Zener transitions as random, Markovian events, transferring particles between neighboring levels with some probability p LZ . The probabilities P {n k } can then be obtained by performing Monte Carlo averaging. First we initialize the fermions on the M lowest levels, then we apply a Markov process in the space of level occupations. During this process, the fermion configuration can only change by nearest neighbor transitions and each such transition happens with rateD(ṽ). Note that this model uses a single parameter extracted from the full quantum simulation, the diffusion constant  [43]. We then determine the work statistics as , and the random matrix average performed only on the final eigenenergies, ε k (t).
As shown in Figs. 2, 3, and 4, apart from the very short time behavior where quantum mechanics rules, the simple SEP approach reproduces the results of our fully quantum mechanical computations with amazing accuracy for slow quenches. This result now allows us to use SEP computations to obtain predictions for the work distribution function in the regime of large injected work, ∼ W > ∼ 20, inaccessible to our quantum mechanical simulations. These results are shown in the rightmost panel of Fig. 2.

VI. MEAN FIELD DESCRIPTION
The SEP model thus provides an accurate description of energy absorption for the most interesting, universal regime,ṽ < ∼ 1, but provides limited analytical understanding. We can construct, however, a simple mean field theory of work distribution, by assuming that the occupations n k are classical binary variables with expectation values f k (t) = n k , and that they are independent -apart from overall particle number conservation. This simple mean field theory incorporates three important ingredients: Pauli principle, particle number conservation, and the diffusive character of Fermi surface broadening. Below we demonstrate that this approach can be used to obtain analytical assymptotic estimates for the probability of adiabaticity and for the variance of work.

A. Probability of adiabaticity
Within the mean field approach, we consider each occupation number n k as a binary probability variable, having values n k = 0 and 1. In the simplest classical approximation neglecting all correlations between different levels, at time t we can assign the probabilities f k (t) and 1-f k (t) to n k = 1 and n k = 0, respectively, Here we set the expectation value of n k , f k (t), in accordance with the diffusively broadened Fermi sea, Eq. (12). To incorporate correlations at the lowest order, we supplement Eq. (15) with a global constraint expressing particle number conservation, k n k = M . For simplicity, here we focus on the half-filled case M ≡ N/2, though calculations can easily be extended to any value of M . We enforce the constraint by inserting a Kronecker δ function into the joint probability distribution, Here N t is a time dependent normalization factor, which can be estimated by first carrying out the summation over {n k }, and then applying the saddle point approximation as N t ≈ π −π dλ 2π j>0 cos 2 (λ/2) + sin 2 (λ/2)erf 2 j 4 Dt The probability of staying in the ground state then reads Here, again, an integral approximation has been made by assuming Dt ≈ ∼ W ≫ 1 and the saddle point approximation has been carried out, leading to the asymptotic estimate, yielding an excellent fit for ∼ W > 1, as shown in Fig. 3.

B. Variance of work
For ∼ W ≫ 1, we can estimate the variance of the work by neglecting the fluctuations of the individual energy levels [50], ε k (t) → ∆k, while keeping track of the fluctuations of the occupation numbers. For a given realization of H(t), this yields the approximate expression the average signs denoting here just quantum averages. By separating the 'diagonal' contributions, the average δ ∼ W 2 RM can be rewritten as with δn k,t ≡n k,t − n k,t denoting the deviation of the occupation number. Since then k,t behave as binary variables, averages in the first term can be expressed as δn 2 The correlators appearing in this equation can be expressed in terms of the amplitudes α m (t) as δn k,t δn k ′ ,t = − Notice that this correction is negative, indicating that the level occupations are anticorrelated, as dictated by particle number conservation. By neglecting for a moment these correlations and replacing sums by integrals, we obtain the estimate reproducing the superdiffusive scaling δ ∼ W 2 ∼ ∼ W 3/2 . We note that while neglecting the correlations explains the observed behavior of the work's variance, it does not reproduce the correct prefactor. Indeed, a more careful mean field calculation along the lines of the previous subsection shows [41] that correlations (the conservation of fermion number) cannot be entirely neglected, and they also give a similar but smaller ∼t 3/2 contribution to the variance -without changing the overall scaling.

VII. CONCLUSIONS
We have derived and used a determinant formula to compute the distribution of work, P t (W ), in course of a quantum quench in generic, disordered, fermionic nanosystems, and have shown that it displays a large degree of universality, especially in the slow quench limit. Even if a complete experimental characterization of P t (W ) may still be challenging, many of our curious findings such as the diffusive broadening of the Fermi energy, the ∼ t 3/2 scaling of the variance of the energy absorbed, or the low frequency ∼ ω 1+β/2 absorption rate, could be readily verified experimentally.
Our results also demonstrate that quantum work statistics in these systems is, after all, to a large extent classical. In the most interesting slow quench limit, we demonstrate a close connection to the symmetrical exclusion process (SEP), a classical diffusion of hard core particles in energy space. Quantum mechanics enters here through level collisions, giving rise to Landau-Zener transitions, the exclusion process, mirroring the Pauli principle and, finally, level statistics, reflecting the symmetry of the underlying quantum mechanical system.
Though diffusion in energy space and its relation to thermalization has been studied by several authors in different contexts [38][39][40], earlier works have not focused on the impact of quantum statistics on energy diffusion. In particular, here we find that energy absorption in this many-body system corresponds to diffusion of strongly interacting particles. This very strong interaction, captured at the classical level by the exclusion process, implies that while the total energy absorbed shows a clear drift in time, W ∼ t, its spread is non-diffusive but rather increases superdiffusively as δW 2 ∼ t 3/2 during a quench.
Our investigations, focusing on T = 0 temperature non-interacting systems, represent clearly only the first step. Inclusion of interactions, generalizations to the symplectic universality class as well as to finite temperatures and open systems are all exciting open questions for future research. and investigate the separation of neighboring levels, ∆(λ) ≡ ε k+1 (λ) − ε k (λ). This function displays welldefined minima, ∆ min , at positions λ 0 , close to which the level separation can be well-described in terms of an effective two-level Hamiltonian as The number of level crossings occurring between two selected neighboring levels is directly proportional to the distance covered within the random matrix manifold, N cross ∼ s tot ∼ √ N dλ. Therefore, we can define the probability distribution of ∆ min as dN cross = ρ(∆ min ) ds d∆ min .
(A3) The data trace out a universal curve, scaling as D ∼ṽ 1+β/2 in the slow quench limit,ṽ < ∼ 1. For fast quenches,ṽ ≫ 1, the statistical aspects of level repulsion lose their role and no longer govern the many-body dynamics, leading to D ∼ṽ 2 for both ensembles. Here the density function ρ(∆ min ) depends implicitly on the size of the random matrices, N , as well as on the universality class considered [46]. The dependence on N , however, appears just in a trivial way in the large N limit, through the level spacing ∆ at the energy of the Landau-Zener transitions. This dependence can be scaled out, yielding a universal distribution function, with ∆ min = ∆ min /∆ denoting the dimensionless Landau-Zener gap. Similarly, the distribution of the slopes γ is also universal, provided we measure it in its natural unit, The numerically computed statistics of the distributions ̺( ∆ min ) and̺(γ) are displayed for the Gaussian orthogonal ensemble (GOE, β = 1) and Gaussian unitary ensemble (GUE, β = 2) in Fig. 6. The numerically extracted distributions fit perfectly the analytical predictions of Ref. 44. In particular, the distribution of the slopes has a peaked structure, i.e., there exist a typical value ofγ, characterizing the transitions. The distribution̺( ∆ min ) is, on the other hand, more peculiar. In particular,̺( ∆ min ) ∼ ∆ β−1 min , and scales as in the orthogonal and unitary universality classes, investigated here in detail [47]. This curious behavior implies that avoided level crossings with very small gap are abundant in the orthogonal class, which leads to a breakdown of adiabatic perturbation theory for β = 1.
The probability p LZ of a Landau-Zener transition between two neighboring levels depends on the velocity at which the transition is approached, and is known analytically [48,49], Given the exponential sensitivity to ∆ min , small Landau-Zener gap transitions dominate in the limit of small velocities,ṽ ≪ 1. This leads to the estimate The universal scaling of the many-body diffusion constant, as extracted from our simulations, is plotted in Fig. 5. For time reversal symmetry breaking (unitary) systems, the diffusion constant scales both at small and high velocities asṽ 2 . In the orthogonal case, however, a clear crossover is demonstrated between a Landau-Zener dominated small velocity regime with D ∼ṽ 3/2 , and a high velocity fast quench regime with D ∼ṽ 2 .
Appendix B: Time evolution of single particle states To solve Eq. (6) in the main text, it is practical to make a simple gauge transformation, and eliminate the phase dependence generated by the instantaneous eigenenergy, ε k (t), with the dynamical phase defined as Importantly, dynamical phases cancel in the expression of the occupation numbers, f k (t), as well in the expression of the work and its generating function. Therefore, in all these equations, we can just replace α k (t) byα k (t), obeying the simpler differential equation, In practical calculations, it is useful to avoid numerical differentiation, and calculate the Berry connection as In our numerics, we made use of the fourth order Runge-Kutta method to solve the single particle time-dependent Schrödinger equation. Dynamical phases in Eq. (B1) were determined numerically using the Simpson-formula. Since the phases of the instantaneous states often make jumps, we enforced the choice A kk = −i η k t · ∂ t η k t = 0 by requiring that the overlaps of two consecutive eigenstates remain close to 1, η k t |η k t+δt ≈ 1.

Appendix C: Fingerprints of level repulsion and deviations from universality
As discussed in the main text, the distribution of the work is universal, i.e., it is independent from the system size N as well as the Fermi energy; it depends solely on the dimensionless time, the dimensionless velocity, and the symmetry class of the random Hamiltonians. This is illustrated in Fig. 7, where we have fixed the dimensionless quench velocity toṽ = 0.8, and the dimensionless times such that they correspond to an injected dimensionless work ∼ W = 1.5, 3, and 5 in both symmetry classes.
For relatively small injected internal energies, W < ∼ 5 ∆, the discreteness of the energy levels becomes apparent, and fingerprints of level repulsion can be observed. First, P ( ∼ W ) vanishes continuously at zero (within SEP as ∼ W β ), since the first empty level is repelled from the last occupied level. Moreover, additional wiggles appear in P ( ∼ W ), reflecting the positions of first, second and third neighbors. These features are more pronounced in the Gaussian unitary ensemble, where level repulsions are stronger, and are expected to become even more pronounced for the symplectic ensemble, not studied here.
The observed distributions are, however, independent of the matrix size, N , as stated above.
Remarkably, even forṽ = 0.8, SEP simulations (continuous lines in Fig. 7) provide an almost perfect description of the full quantum results. For larger velocities, however, deviations occur. This is demonstrated in Fig. 8. Clearly, forṽ > ∼ 1 the role of Landau-Zener transitions is reduced, the SEP description loses its validity, and distributions depend explicitly on the velocityṽ.