Statistics of the work distribution for a quenched Fermi gas

The local quench of a Fermi gas, giving rise to the Fermi edge singularity and the Anderson orthogonality catastrophe, is a rare example of an analytically tractable out of equilibrium problem in condensed matter. It describes the universal physics which occurs when a localized scattering potential is suddenly introduced in a Fermi sea leading to a brutal disturbance of the quantum state. It has recently been proposed that the effect could be efficiently simulated in a controlled manner using the tunability of ultra-cold atoms. In this work, we analyze the quench problem in a gas of trapped ultra-cold fermions from a thermodynamic perspective using the full statistics of the so called work distribution. The statistics of work are shown to provide an accurate insight into the fundamental physics of the process.


Introduction
In the past decade ultra-cold quantum gases have emerged as ideal candidates for clean and controllable simulation of condensed matter physics [1]. Ultra-cold quantum gases are now created in a variety of configurations in laboratories worldwide. In particular, both Bosonic and Fermionic atoms can be trapped and manipulated on optical lattice potentials [2]. The lack of thermal phonons coupled with the tunability of the interactions by means of Feschbach resonances [3] has allowed for the detailed study of a multitude of phase diagrams, the most celebrated example is perhaps the Bose-Hubbard model [4,5,6]. Equilibrium properties aside, over the past number of years there has been a surge in interest in the out of equilibrium behavior of closed quantum systems following a quench of a Hamiltonian parameter. Fundamentally, this is due to a series of spectacular experiments in ultra-cold atoms whereby the high degree of isolation and long coherence times permits the study of dynamics over long timescales [5,7]. These experiments have raised a number of important theoretical issues such as the relationship between thermalisation and integrability and the universality of defect generation following evolution across a critical point [8]. Given the controllability of ultra-cold atomic systems and the current interest in quench dynamics, it is a natural question to ask if there are any out of equilibrium condensed matter physics problems which could be simulated. Unfortunately, due to their intrinsic complexity, there are very few examples of exactly solvable problems in the out-of-equilibrium domain. A worthy exception is the phenomenon of orthogonality catastrophe [9,10,11] and the Fermi-edge singularity. This problem was first pointed out by P. W. Anderson over 40 years ago when he showed that the overlap of two manybody wave-functions, which describe deformed and undeformed Fermi seas, vanishes in the thermodynamic limit [9]. The corresponding 'quench' problem, was investigated a few years later with the prediction of a universal absorption-edge singularity in the X-ray spectrum of metals, the 'Fermi-edge' singularity [12]. The universal physics of the Anderson orthogonality catastrophe and the Fermi-edge singularity was recently explored by Goold et al. in the context of ultra-cold quantum gases [13]. In this work it was suggested that the physics maybe simulated in a controlled fashion by the appropriate embedding of a single probe qubit. The approach was further formalized in [14] where the authors solved the dynamical problem in the inhomogeneous system by means of a linked cluster expansion. This qubit probe approach was further suggested as a mechanism to probe the physics of the orthogonality catastrophe in [15,16]. Interestingly enough, a connection was made by Heyl and Kehrein [17] between the absorption and emission spectrum in the original X-ray experiments and the so called quantum work distribution and corresponding fluctuation relations in classical and quantum statistical mechanics [18,19]. Treating a quench problem in manybody physics as a thermodynamic transformation and analyzing the statistics of work done has recently shown to be a useful approach to understand the intrinsic out of equilibrium dynamics in manybody systems [20,21,22,23,24]. The approach is based on the study of the moments of a quantity known as the quantum work distribution [25,26] which have been found to encode both thermodynamic and universal features of the model in question. In this work this relationship will be explored in detail in the context of a locally quenched trapped Fermi gas. In particular, in section 2 the general problem of a system Hamiltonian depending on a work-parameter will be introduced, adopting the description based on the grand canonical ensemble [27]. The formalism provided by the work distribution and its characteristic function will be discussed, and the concept of irreversible work will be explored. In Section 3 the focus will be moved to a Fermi gas in equilibrium with a harmonic trap, being suddenly perturbed by a spatially structureless perturbation. In section 4 the vacuum persistence amplitude and the linked cluster expansion will be used to reduce the calculation of the characteristic function of work to the sum of connected Feynman diagrams. The relation of these diagrams to the characteristic function of work will be covered and an analytic approximation holding at low temperature will be presented. In section 5 the first three cumulants of the work distribution will be computed and their link with thermodynamics will be discussed. Finally, in section 6 the irreversible work will be calculated using a perturbative and a numerical approaches, while Section 7 will provide some further comments and conclusions.

Non-equilibrium quantum thermodynamics
Consider a system evolving according to the HamiltonianĤ(η), which depends on some externally tunable parameter η. The system is brought into weak contact with a heat reservoir, at inverse thermal energy β, and allowed to equilibrate for each chosen value of η, before the system-reservoir coupling is turned off. To each η there corresponds a well defined Gibbs statê being the grand-canonical partition function andN the particle number operator. Suppose that η is some explicitly time-dependent degree of freedom, which interacts directly with the system and not with the reservoir, i.e., a work parameter η = η t . Then, for an initial value η 0 , at the time t = 0, the state of the system isρ(η 0 ). The initial Hamiltonian and particle number operators have, respectively, the spectral decompositionsĤ where |n is the n th simultaneous eigenstate ofĤ(η 0 ) andN, with eigenvalues E n (η 0 ) and N n . Next, some 'work' is performed on the system taking the work parameter from η 0 to a final value η τ , at a later time t = τ . The final Hamiltonian, connected by the protocol η 0 → η τ , has the spectral decompositionŝ where |m is the m th simultaneous eigenstate ofĤ(η τ ) andN, with eigenvalues E ′ m (η τ ) and N m .
The definition of work in this scenario requires two projective measurements: the first projects onto the eigenbasis of the initial HamiltonianĤ(η 0 ), with the system in thermal equilibrium. The system then evolves under the unitary dynamics U(τ, 0), generated by the protocol η 0 → η τ , before the second measurement projects onto the eigenbasis of the final HamiltonianĤ(η τ ). The probability of obtaining E n (η 0 ) for the first measurement outcome followed by E ′ m (η τ ) for the second is then The work distribution is defined as [25,26] and the average work W done on the system is given by the first moment of P η 0 →ητ (W ). It is useful to introduce the characteristic function of work [26] as with · · · = tr[· · ·ρ(η 0 )] denoting the thermal equilibrium average over the initial state. In terms of χ(t, τ ) the average work is expressed as W = −i dχ(t, τ )/dt| t=0 . Microscopically, the second law of thermodynamics is revised to the form with Ω(η) = − ln[Z(η)]/β being the thermodynamic grand potential, so as to encompass the explicit statistical nature of work. The deficit between average work and the variation in the grand potential can be accounted for by the introduction of the irreversible work contribution W IRR > 0, being such that W = ∆Ω + W IRR . When combined with the first law of thermodynamics, this relation can be rewritten as where ∆S is the change in entropy of the system and ∆S IRR = β W IRR (∆S REV = β ∆Q ) is the irreversible (reversible) entropy change. We note that for a closed quantum system, the heat transfer into the system is zero and the sole contribution to the entropy change is the irreversible entropy. To these concepts and their link to the irreversibility of a sudden transformation, we shall come back in Sec. 6, after we have evaluated the characteristic function of the work distribution and its first moments for the specific problem at hand.

Model Hamiltonian
Consider a gas of non-interacting cold fermions, confined by a one-dimensional trapping harmonic potential. The trap has a characteristic length x 0 and a resonance frequency ω. Then, each particle of the gas, in its unperturbed state, is described by the harmonic oscillator Hamiltonian The unperturbed system Hamiltonian is spanned by the fermion field in real spaceΨ ξ (x), with ξ denoting the spin degrees of freedom, i.e., ξ ranging over (2s + 1) values.
Let us now consider doing work on the Fermi gas by the turning on of an external potential. Let us further assume that the switch is done in a sudden way. In this paper we will take this external perturbation to be a localized potential at the centre of the trap which we can model by a Dirac δ function with strength V 0 such that V (x, t) = θ(t)πV 0 x 0 δ(x), where θ(t) is the Heaviside step function. Then, the total Hamiltonian for the gas after the switching on of the potential isĤ =Ĥ 0 +V , wherê 3.2. The characteristic function of work and the vacuum persistence amplitude Consider the work parameter η t = θ(t)πV 0 x 0 and the protocol η − → η + , from t → 0 − (with η − = 0) to t ≥ 0 + (with η + = πV 0 x 0 ). Then, the initial state of the systemρ(η − ) is given by Eq. (1) with the initial HamiltonianĤ(η − ) =Ĥ 0 , introduced in Eq. (6), and the particle number operatorN = ξΨ † ξ (0)Ψ ξ (0). The partition function in the initial state will be denoted Z. The final HamiltonianĤ(η + ) =Ĥ 0 +V includes the perturbation operator introduced in Eq. (7).
In this particular case a sudden quench occurs, and the characteristic function of work (4) reduces to [26] In the dynamical response theory of many particle system, a key quantity to calculate is the so-called vacuum persistence amplitude [28] i.e., the probability amplitude that the gas will retrieve its equilibrium state at time t, after the switching on of the perturbation. It is quite obvious that a simple relationship exists between the characteristic function (8) and the vacuum persistence amplitude (9): Exploiting this relationship we have the full access to the statistics of work done via calculation of the vacuum persistence amplitude. Initially the fermions lie in their equilibrium configuration, set byĤ 0 , until a sudden perturbationV (t) =V θ(t) is felt by the gas. The initial equilibrium depends on the harmonic oscillator Hamiltonian given in Eq. (5) whose unperturbed eigenfunctions are with H n (x/x 0 ) being the Hermite polynomials of order n. By the parity of these wave functions, the external potential induces excitations which connect only unperturbed one-fermion states labeled by even numbers n = 2r, with r = 0, 1, · · · , ∞. Then, the matrix elements of the external potential in the unperturbed oscillator basis read where we have introduced the Euler gamma function ratio γ r = Γ(r + 1/2)/Γ(r + 1). We express the fermion field in terms of the annihilation operatorĉ nξ for the (unperturbed) n-th single particle state of energy ε n = ω(n + 1/2) and spin ξ. Then, usinĝ with ψ ξ being the spinor wave function, the unperturbed Hamiltonian results in the particle number operator readsN = n,ξĉ † nξĉ nξ , while the effect of the perturbation on the gas is represented bŷ

Linked Cluster expansion
A typical approach in many-body physics to find a manageable expression for the vacuum persistence amplitude (9) is to turn to the interaction picture, where the impurity potential reads In this representation, we may use the identity (10) to express the characteristic function of work as Then, we may perform a linked cluster expansion and reduce χ β (t) to an exponential sum of connected Feynmann diagrams: ... . (13) These loops contain separable products of lines connected to vertices (V rr ′ ), with level occupation numbers given by the Fermi-Dirac distributions As is standard in discrete-level systems, and intrinsic semiconductors, we let the chemical potential µ lie in the middle between the Fermi energy ε F and the lowest unoccupied one-fermion state at the absolute zero (β → ∞). Without loss of generality, we may assume the Fermi level number n F to be even, label it by 2r F (r F being a positive integer) so that ε F = ω(2r F + 1/2). Then, we can parametrize µ as µ = ω(2r µ + 1/2), i.e., f ± r = [1+e ±2β ω(r−rµ) ] −1 , and compute r µ for finite β by constraining the average particle number to be: From the plots of Fig. 1, we observe that µ takes its maximum value at β ω → ∞. Then, it decreases with decreasing β ω, and it becomes largely negative for very low β ω where the classical limit applies. We now apply Wick's theorem and focus on the lowest order loops in (13), i.e., Λ β 1 (t) and Λ β 2 (t). In terms of the auxiliary functions we express these contributions as It is useful for the following to parameterize V 0 in terms of the energy scale √ 2 ωε F , which results from the geometric mean between the spacing of even energy levels and the Fermi energy. To this end, we introduce as a sensible parameter, and to keep the interaction potential small, we investigate the range α = 0 − 1. It turns out that, in the free gas limit (ω → 0), this interaction strength coefficient reduces to the so called "critical parameter" of the Mahan Nozières De Dominicis (MND) theory of the edge singularity [11,12], which, in the X-ray absorption problem, also gives a measure of the asymmetry of the absorption spectrum. Analogously, we will find below, for our trapped gas case, that the parameter α determines the skewness of the work distribution.

Connected diagrams and work distribution
The full derivation of Eqs. (19) and (20), and the numerical calculations involved, can be found in a recent paper by the authors [14], and the reader interested in specific details is directed there. In the present context, we quote the main results and apply them to the determination of the work distribution P η − →η + (W ) and its characteristic function (13). The one-vertex loop is just the adiabatic response of the gas, being of the form represents the first-order correction to the equilibrium energy, given by Eq. (22) brings a phase factor to χ β (t), which corresponds to shifting the work distribution P η − →η + (W ) by E 1 β . The two-vertex loop can be split into three parts with well defined trends and physical meaning, namely, The first one is Λ β provides the second-order correction to E β 0 and brings a further shift to which produces a Gaussian damping in χ β (t) and a Gaussian broadening in P η − →η + (W ). The third one accounts for the shake-up of the gas due to the sudden switching of the impurity. It is a periodic function of time with frequency 2ω (as a direct consequence of the harmonic form of trapping potential), which provides the non trivial part of P η − →η + (W ). The basic quantities of the problem, given in Eqs. (22), (24), (25), and (26), contain summations running over all one-fermion eigenstates of the trap, weighted by the Fermi factors f ± r , and enter the characteristic function of work (10): By Eq. (3), the work distribution is given by the convolution product Numerical computations of P η − →η + (W ) are shown in Fig. 2, where we recognize an asymmetric, broadened profile, signature of the singular behavior of the Fermi gas. The monotonic structure turns into a satellite structure of sub-peaks, separated by 2 ω and related to even-level transitions in the gas, as β ω gets above ∼ 0.5.  These secondary peaks are a direct manifestation of the single particle transitions occurring within the gas, with the fermions jumping between even harmonic oscillator states (separated in energy by even multiples of 2 ω), in response to the local perturbation. All of these transition constitute the so called "shake-up" process, which is, thus, explicitly witnessed by the work distribution.
At high temperatures, this feature disappears, being hidden by the increased gaussian broadening of the main peak, which (in our two-vertex approximation) is found at W = E β 1 + E β 2 , and which describes the transition between the equilibrium states induced by the switching of the scattering center.

Low thermal energy approximation
Consider now a Fermi gas with a large number of particles at sufficiently low temperatures, being such that we may approximate the chemical potential by Eq. (17). In this regime, we can expand f ± r in power series of e ±β ω(ε 2r −µ∞) and select the lowest order of this series. As detailed in the appendix, we can work in the r F ≫ 1-limit and find some manageable expressions for the auxiliary functions λ β ± (t) introduced in Eq. (18). By Eq. (19), the initial value λ β + (0) determines the one-vertex loop and hence the first order energy shifts (22). On the other hand, the product λ β + (t)λ β − (t) enters a double time-ordered integral which gives (20), together with the second-order energy shifts (24), the Gaussian standard deviation (25), and the shake up sub-diagram (26). This integral can be carried out analytically, in the r F ≫ 1-limit, by adding a small shift to the time domain on the imaginary axis. We, therefore, obtain the following analytical approximation for the characteristic function (27) where Here, τ 0 is a regularization parameter, i.e., an imaginary time shift, that we may interpret as the typical time interval needed by the system to respond to the abrupt switching on of the impurity potential at zero temperature. We see that the function (29) is undefined in the τ 0 → 0-limit, which is a consequence of the large r F expansion used in the thermal series for λ β ± (t). The exact expressions provided by Eqs. (19) and (20) are well set and can be used to compute (27) and then (28) numerically, without any divergence problem, as we did in Fig. 2.
We remark that the imaginary time regularization leading to Eq. (29) is required to provide an analytical support to the theory. Furthermore, it has been observed that [14] correctly tends to the Nozières and De Dominicis core hole propagator [12] when the harmonic trap frequency is lowered to zero, keeping the number of particles in the gas finite. In this limit, the regularization parameter τ 0 , becomes exactly the one needed in the MND theory [29]. The mathematical details of the derivation of Eqs. (30), (31), and (32) are discussed in Appendix A. The leading behavior of χ β in Eq.(29) vs temperature is provided by the gaussian damping factor, whereas both the energy shifts and the periodic part of the characteristic function are well approximated by their absolute zero expressions in a wide range of temperatures, corresponding to β ω 0.2 for r F 10. As shown in Fig. 3A, the firstorder shift E ∞ 1 , given by Eq. (30), is the correct large β ω-limit for E β 1 , numerically calculated from Eq. (22). To compute the low temperature form of δ β , we truncate the series in Eq.   30) and (31), respectively. The critical exponent α is fixed to α = 0.2, while several particle numbers are considered in the range N = 22 − 1002, i.e., r F = 5 − 500.
We see that already the M = 1 curve accurately reproduces the numerical data for β ω 5. The approximation with M = 100 components works particularly well in the extended range β ω 0.5, where the numerical δ β curves are independent on r F . Now, consider the regularised shift E ∞ 2 , as introduced in Eq. (30), and the regularised shake-up diagram Λ ∞ 2P (t), reported in Eq. (32). We need to adjust τ 0 for each r F in such a way that these two quantities are the correct absolute zero limits for the corresponding numerical quantities, i.e., It turns out that τ 0 decreases with increasing r F , which makes sense because the more particles the system has the more allowed transitions are offered to respond to the sudden perturbation (Fig. 4A). However, the set of τ 0 values that accurately fit the condition (33) are slightly different from those that realize the condition (34). This is because the asymptotic forms of E ∞ 2 and Λ ∞ 2P (t) contain terms proportional to r . ., which we have neglected working in the large r F limit (Appendix A). As a reasonable compromise, we take the average between these two optimized sets (Fig. 4A), which is an agreement with both conditions (Fig. 4A, C, D) within an error less than 5%. Interestingly, E β 2 and E β 1 have similar trends and absolute values for α = 0.2, while E β 2 is more sensitive to temperature than E β 1 for β ω < 0.05. On the other hand, the modulus of the sub-diagram Λ β 2P presents zeroes at ωt = kπ and maxima at ωt = kπ/2, with k = 0, ±1, ±2, · · ·. The intensities of such maxima (Fig. 4C) increase with increasing the Fermi number (2r F ), the critical exponent (α), and the thermal energy (β −1 ). The phase of Λ β 2P is discontinuous at the extremes of |Λ β 2P | and less dependent on these parameters.

Cumulant expansion
We now have all the ingredients to start our discussion of the statistical properties of the work distribution of a non-interacting Fermi gas held in a harmonic trap, following a sudden local quench of a point like scattering potential [14]. In the present section we will provide a physical interpretation for the features of the work distribution. In particular, we will work out the cumulant expansion of the work distribution and look for the link between different cumulants and thermodynamical quantities. The characteristic function obtained in Eq. (8) has been introduced in Eq. (4) as the mean value of the random variable e iW t/ in the work distribution P η − →η + (W ) or, equivalently, as the Fourier transform of the work distribution. The work distribution, then, is nothing but the absorption spectrum of the system due to the suddenly switchedon impurity potential [14]: it has been set in Eq. (28) and shown in Fig. 2. We notice that because of the assumed normalization of P η − →η + (W ), we have χ β (0) = 1. We call moment of order n, or equivalently n−th moment of the distribution, the mean value W n . Once the characteristic function is known, we can use the differentiation theorem of Fourier transforms to evaluate each of the above defined moments as This relation holds provided that χ β (t) is continuous and differentiable n-times, with all of the derivatives vanishing at t → ±∞ (that we will see not to be always the case). It will be more convenient for us to work with the cumulant expansion of ln χ β (t) instead of computing the moments of χ β (t). The cumulants are defined analogously to the moments in Eq. 35 with ln χ β (t) replacing χ β (t). This will have two advantages. First it will be easier to characterize the distribution since, as we will see, ln χ β (t) makes it possible to calculate important quantities such as the mean value, variance and skewness straightforwardly. Second, in our case χ β (t) = Π n e Λ β * n (t) so that we will be able to link the properties of the distribution (as given by its mean value, variance and skewness) to the dynamical functions Λ β n (t). Using the form given by Eq. (8) for the characteristic function, we express the cumulant generating function as The first three cumulants are, then, given by where we defined the variance σ 2 and the skewness κ. In our case, and by considering only two sets of diagrams, see Eq. (36), the above quantities are simply calculated as We see that the two connected diagrams found previously allow us to calculate the first three cumulants in a straightforward manner. In particular they are determined by the analytical properties of the periodic sub-diagram (26) at t → 0. Now, while Λ β 2P (t) is a continuous and differentiable function, its higher order time derivatives, i.e., d n Λ β 2P (t)/dt n for n > 1, are ill-defined in the t → 0-limit. As a corollary, we have that the higher-order time derivatives of the characteristic function are not defined at t = 0, then Eq. (35) lacks formal justification for n ≥ 2. This is a direct consequence of the point-like modeling of the impurity potential. A possible work around of the problem will be proposed in the following paragraphs by introducing suitable cut-off frequencies on the perturbation matrix elements (11). In particular, we shall use the fact that the regularized characteristic function (29) admits cumulants of any order to renormalize the second and third ones given in Eqs. (41) and (42). Indeed, the periodic subdiagram quoted in Eq. (29) has been proved to accurately reproduce the numerical expression (26) in the absolute zero-limit (Fig. 4). Then, it makes sense to match the cumulants, as given by Eq. (29), with the corresponding quantities obtained from the numerical form of ln χ β (t), Eq. (27). Such a condition does not alter the physics of the work distribution, while it will allow us to carry out the analysis of the physical behavior of the skewness, which will turn out to have an expression independent of this regularization procedure.

The mean value
First let us consider the mean value (40). Using the expressions for E β 2 and Λ β 2 (t) given in Eqs. (24) and (26), respectively, we may write Then, considering the identity Interestingly Eq. (43) states that the mean work is given by the first-order energy shift of the quenched Fermi gas. Such a relation continues to hold at absolute zero with the regularized approximations given by Eq. (29), (30), and (31). Indeed, the second order time derivative of the regularized periodic sub-diagram reads so that, using Eq. (40) again, we get: κ 1 (∞) = E ∞ 1 . Although the calculations leading Eq. (43) have been derived from an approximated expression for the vacuum persistence amplitude, we may prove Eq. (43) to be formally exact. To see this let us consider the expression (12) for the characteristic function together with the definition of the moments in Eq. (35): This result tells us that for a sudden quench the average work done is the mean value of the perturbation after the switch-on at instant t = 0. However, this is also the first-order energy reported in Eq. (22): The behavior of E β 1 , and hence of the average work W vs the inverse thermal energy β, as been thoroughly discussed in Sec. 4.2 and shown in Fig. 2A

The variance
Using Eq. (41), we may now evaluate the variance of the work distribution: Here we replace δ β with its explicit form given in Eq. (25), and compute the second order time derivative which leads the compact expression To provide an interpretation to this relation, we replace the definition of α and the expression (11) for the impurity potential matrix elements: By analogy with the Fermi's Golden rule, we can look at the quantity with ρ 0 = 2 −1/2 −1/2 ω −1/2 denoting the density of even fermion states in the continuous limit. We recognize T r→r ′ to be the rate of transition from the occupied one-particle state |2r to the empty state |2r ′ . Following this interpretation, we rewrite the variance as It is thus clear that σ 2 gives information on the broadening due to Fermions scattering in empty states and thus on the spectrum of the system. What is also evident from Eq. (44) is that the sum over unoccupied particle states, i.e., the r ′ -series, does not converge, because the perturbation V rr ′ is a weakly decreasing function of r and r ′ . Indeed, σ 2 is proportional to the initial product of the auxiliary functions λ β ± (t), introduced in Eq. (18): Now the asymptotic behaviors f − r ≈ 1 and γ r ≈ r −1/2 for r ≫ 1 lead to λ β − (0) → ∞ (see Appendix A). On the other hand, we can substitute regularized characteristic function (29) in Eq. (41) to get the asymptotic trend This result, combined with the fact that g β vanishes at the absolute zero, gives We see that the divergent behavior of σ 2 in Eq. (47) is absorbed by the regularization parameter τ 0 . Indeed, considering gases with very large particle numbers, i.e., working in the τ 0 ω ≪ 1-limit, we find κ 2 (∞) ≈ α 2 τ −2 0 . To recover a consistent definitions of κ 2 (β), as given by Eq. (41) and hence Eq. (47), we proceed similarly to the Nozières and De Dominicis [12] work around of the Fermi-edge singularity; we introduce an exponential frequency cut-off on the impurity potential and change the γ r factors to γ r e −2rω/ω 0 . This is equivalent to adding a lifetime width 1/ω 0 to the unperturbed propagator (14). By doing so, Eq. (44) becomes well defined and factorisable as and the auxiliary functions are evaluated on the complex time-domain. In Appendix A we show the cuf-off frequencies ω 0 is related to the regularization time τ 0 by the condition that κ 2 (∞), as calculated from Eq. (50) with β → ∞, matches with the regularized expression (49).

The skewness
The skewness κ = κ 3 (β) /κ 2 (β) 3/2 is related to both the second an the third cumulant of the distribution. To begin, we address our attention to the third cumulant (42), i.e., , keeping in mind that κ 2 (β) is a non negative quantity. We then compute the third order time derivative Next, we use the definitions for the impurity potential matrix elements and the transition rates, as in Eq. (46), and express Thus, for the skewness we get which encompasses a direct thermodynamical meaning. The first term is nothing but the energy taken from the system by emptying its states whereas the second is the energy given to the system by filling its empty states following the thermodynamic transformation, which in our case is a sudden quench. So this quantity tells us whether the transformation effect is to increase or decrease the internal energy of the gas and is clearly related to the asymmetry of the work distribution To provide a consistent definition of κ 3 (β), we need to tackle the divergent behavior of the r ′ -series in Eq. (52), which is more evident by expressing the third cumulant as To avoid this divergence, we adopt the previous procedure once again and redefine Eq. (52) by adding an imaginary time-shift to the auxiliary functions λ β ± (t), i.e., The new cuf-off frequency ω ′ 0 is fixed by the condition that κ 3 (∞), as calculated from Eq. (54) with β → ∞, tends to the regularized expression computed with the absolute zero distribution (29). In Appendix A we provide details on how the renormalization of the second and third cumulants of the work distribution, given by Eqs. (50) and (54), are carried out. In Fig, 5A we show the cut-off times 1/ω 0 and 1/ω ′ 0 that let us match the numerical behavior of κ 2 (β) and κ 3 (β) with the expressions given by Eqs. (49) and (55), respectively. Using such values we can obtain the behavior of the variance and the skewness vs β ω for different values of r F (Fig, 5B  and 5D). The dependence of the second and third moment on α is trivial, as both κ 2 (β) and κ 3 (β) are directly proportional to the critical index. Then, the skewness goes like which does not suffer from the divergent behavior of κ 2 (∞) and κ 3 (∞). Furthermore, in the τ 0 ω ≪ 1-limit, the third cumulant behaves as κ 3 (∞) ≈ 2α 3 τ −3 0 and the skewness turns out to depend only on the critical exponent κ ≈ 2/ √ α (Fig. 5D). Being completely independent of the regularization procedure, this is a physically meaningful result, showing how the critical parameter α does indeed determine the asymmetry of the work distribution, ultimately due to the Fermi edge behavior.

Irreversible work
In the previous section we have seen that the work distribution contains information about both the unperturbed system and the system following the thermodynamic transformation. In this section we shall first recover a Jarzynski-like equality, which will help us in identifying the hypothetical final equilibrium state of the system to compare with. We shall then find an approximate and analytic expression, and discuss a numerical method for the computation of the irreversible work, which is a figure of merit of the irreversibility of the transformation. We are interested in highlighting the microscopic origin of irreversibility for the sudden quench studied above. First let us recover the Jarzynski-like equality. We start from the evaluation of e −βW by Eq. (3), i.e., In going from the first to the second row we have used particle number conservation, p m|n ∝ δ Nn,N ′ m , and the completeness relation n p m|n = 1, [25]. We recognize that the n-sum in Eq. (57) is nothing but the grand canonical partition function for the gas in equilibrium with the impurity potential, with same temperature and chemical potential as the initial one. Hence, we have: where and denote the grand potentials for the unperturbed and perturbed equilibrium states, respectively. The fact that the Jarzynski relation makes explicit connection with a (hypothetical) final equilibrium state is meaningful as it gives us a reference for the (hypothetical) reversible version of the transformation encompassed by the switching on of the external perturbation. The reversible version of this protocol, thus, would involve a change in the number of particles, which is necessary in order to maintain the initial chemical potential even if the single particle energies are modified during the protocol [27]. By means of the Jensen inequality e −βW ≥ e −β W we can derive the statement of the second law of thermodynamics : As expected, the average work done on the system is greater than the change in the grand potential of the initial state and the hypothetical final equilibrium state. This relation suggests us to define a new variable that we shall name W IRR = W − ∆Ω, with ∆Ω = Ω ′ − Ω. The new variable will have the very same distribution as the original one but for the mean value which will be shifted by an amount ∆Ω: In the following, we propose an analytic and a numerical approaches to compute this excess work (62) and discuss it in order to highlight the nature of irreversibility in our system. We have already set all the elements to compute Ω, i.e., the unperturbed level energies ε 2r , distribution functions f ± r , and chemical potential µ, and we have already discussed the behavior of the average work W (see Fig. 2A). What is left to calculate is the perturbed grand potential Ω ′ , for which we will provide both an approximate analytic and a numerical computation.

Perturbation approach
The energies entering the hypothetical final state belong to the spectrum of the Hamiltonian which describes a particle of the gas in equilibrium with both the harmonic trap and the impurity potential. As discussed in Sec. 3, the odd single particle energy levels are left unperturbed ε ′ 2r+1 = ε 2r+1 . As for the even energies ε 2r , which cannot be calculated analytically, we resort to perturbation theory assuming the height πV 0 x 0 of the δ-potential to be small with respect to the unperturbed energies. Then, we express: 2r is computed from the non-degenerate Rayleigh-Schrödinger perturbation theory based on the assumption Substituting (64) in the perturbed grand potential (60), and carrying out a power series expansion for small δε 2r , to the second order we get: Now, considering the assumption (66), we find By Eq. (65), the first term at the right hand side of this relation is simply the average work W calculated above (see Sec. 5.1). The second term is the second order contribution ∆U (2) to the change in the total energy of the system. The third one may be expressed in terms of the Gaussian broadening (25) as βδ 2 β 2 ω 2 /2. We are now able to calculate the excess work (62) as Notice that by Eq. (65) the second-order corrections ε (2) 2r are negative, which makes −∆U (2) a positive quantity, and let the Jensen inequality (61) be always verified. Eq. (68) can be brought into a more interesting and general form, by noticing that We thus conclude that the irreversible work thus takes the suggestive form: This relation is independent on the model used to characterize the impurity potential. Interestingly, we may look at the occupation numbers n 2r ξ = 0, 1 as independent random variables, distributed according to the probability distribution e −β(ε 2r −µ)n 2r ξ 1+e −β(ε 2r −µ) , the function of the configurations {n r ξ } is a random variable too. As a rsult, the more peaked this random function is the smaller its contribution to the irreversible work. This means that for an adiabatic change there is no spread of the work distribution, since in that limit the average work done would be a Dirac delta function of ∆Ω.

Numerical Approach
We now turn the attention to the numerical calculation of the perturbed grand potential (60), which requires the knowledge of the eigenvalues of the singleparticle Hamiltonian (63). The odd harmonic oscillator eigenfunctions ψ 2r+1 (x) and eigenenergies ε 2r+1 are left unaffected by the δ-potential, due to the fact that ψ 2r+1 (0) = 0 (see Sec. 3). On the other hand the perturbed even eigenfunctions of (63), with the physically correct asymptotic behavior, are the parabolic cylinder functions with associated level energies ε 2r = ω(2r + 1/2). The latter are the just perturbed energies denoted ε ′ 2r above. As shown for example in Ref. [30], the stationary Shrödinger equation for the δ-potential implies The last two relations yield an implicit condition between the strength of the δ-function and the quantum numbersr: Since the Γ-function has poles for negative integer values, Eq. (72) leads tor → r for V 0 → 0, andr → r + 1/2 for V 0 → ∞. Then, the energy eigenvalues ε 2r converge to the unperturbed energies ε 2r when the potential barrier is set to zero, while they tend to ε 2r+1 for an infinite barrier. In the latter case, the perturbed energies become identical to the values of the odd eigenfunctions, leading to a double degeneracy of all eigenvalues.
For arbitrary values of V 0 , which means for α = 0 − 1 at any r F > 0, see Eq. (21), we need to solve Eq. (72) numerically and get a sequence of quantum numbersr =r(α, r F ), with r≤r ≤ r + 1/2. Here, it is important then to notice that for each fixed value of α and r F there exists a one to one correspondence between each unperturbed quantum number r and ar. Once ther are known, we can compute the perturbed energies ε 2r and the normalization constants ηr = η(α, r F ), which let us determine the perturbed wavefunctions ψ 2r (x). In Fig. 6A and 6B, we show how the single particle ground state wave functions and the energy levels change with increasing both α and r F . The critical parameter is allowed to take the values α = 0.1, 0.4, 0.6, in a spin-1/2 gas low and large particle numbers, i.e., r F = 5, 100. In panels A and C the Jarzynski inequality (61) is proved to hold for each sampled values of α and r F . In panels B and D the excess of work done by the system is shown to reach a value independent on the temperature of the system for β ω 0.01 − 0.1, depending on r F . Such a value increases with both r F and α.
Within this framework, the perturbed grand-canonical potential is given by The irreversible work can be thus computed from (62) using β∆Ω = (2s + 1) In Fig. 7A and 7C we see that the second law (61) is obeyed for any chosen valued of α and r F . The amount of irreversible work, shown in Fig. 7B and 7D, increases with decreasing temperature reaching a nearly constant saturation value for low enough temperatures, β ω 0.01 − 0.1 depending on the number of particles in the gas. A similar trend is followed by other crucial quantities discussed here, such as the chemical potential, the energy shifts, and the cumulants of the work distribution function. On the other hand, as suggested by the approximation (68), the irreversible work is highly sensitive to the total particle number and the critical exponent following the proportionality relation W IRR ∝ αr F , for r F ≫ 1.

Conclusions
In this work, we explored the physics and thermodynamics of an inhomogeneous Fermi gas perturbed by the sudden switch on of a local scattering potential. Exploiting the direct relationship between the characteristic function of work and the vacuum persistence amplitude, and by means of the linked cluster expansion technique, we obtained the full statistics of work done by performing such a local quench on the gas, showing that the first, second and third moments of the work distribution encapsulate the salient thermodynamic features of the model. Indeed, features of the textbook Fermi-edge singularity problem were found to be present in the higher moments. Furthermore, we obtained analytic and perturbative expressions for the excess or irreversible work done on the system as well as a Jarzynski like equality. From a general perspective, our work demonstrates the potential of combining ideas coming from recent developed techniques in non-equilibrium statistical mechanics with traditional approaches from many body physics for the analysis of out of equilibrium problems in quantum systems. It is also worth pointing out that the approach developed here is by no means restricted to a specific system and can be extended to other types of quench problems, global and local, in non-interacting and interacting many-particle systems. In addition to the theoretical framework developed, one can imagine that in the future the work distribution of the Fermi gas maybe extracted in an experimental setting by means of coupling to an auxiliary ancilla system which would then function as a probe, giving access to the relevant thermodynamic quantities by monitoring, e.g. its decoherence dynamics [31]. This was first suggested in the context of Fermi gases by Goold et al in [13]. In fact very recent proposals to verify the quantum fluctuation relations by means of interferometry of an ancillary probe qubit [32,33] have been realised in a recent experiment in a NMR setting [34]. The challenge for the future is to really push these experiments to the many-body domain. which by Eq. (23) contains all other basic quantities of the problem, see Eqs. (24), (25), and (26). To compute the Gaussian term, however, it is more straightforward to work on expression given in Eq. (25) and use the expansion (A.2). By doing so, we get in which As a first approximation, we focus on the temperature range where the chemical potential is approximated by Eq. (17), i.e., β ω 0.2 for r F > 5. In this regime, we apply the sum rules withF 2 1 being the regularised Hypergeometric functioñ Substituting into Eqs. (A.4) and (A.5), we find 1 − e 2iωt − e 2iωt(r F +1)F 2 1 r F + 1, e 2iωt , (A.10) ∓ e ±2iωt(r F +1) e 3τmω/2F 2 1 r F + 1, e 2ω(±it+τm) + e −3τmω/2F 2 1 r F + 1, e 2ω(±it−τm) As a second approximation, we take systems with large numbers of particles, and employ the large-r F expansions: so that we straightforwardly get , (A.14) .
This expression is largely dominated by the β-independent value, reported in Eq. (30) and shown in Fig. 3A. Turning to the Gaussian coefficient, a simple change of summation indices, with µ = µ ∞ , leads to The transformed summations in this last line are dominated by low r terms. In a many fermion environment, by the asymptotic relation (A.13), we approximate and neglect terms going like e −2ωτmr F , to obtain This relation, being identical to Eq. (31), allows us to express express Then, the standard deviation at low temperatures is independent of the number of particles in the gas (Fig. 3B).
As for the second-order shift Λ β 2S (t) and the Fermi-edge component Λ β 2P (t), we first consider the β-independent coefficients (A.15) of the series (A.3), write down the product and plug it into the time-ordered integrals (A.7). Then, we need to deal with the short-time singularity of Eq. (A.18) by adding an imaginary time regularization to the t ′′ -integral, i.e., we need to shift the t ′′ integration domain by iτ 0 . The resulting integral is dominated by the Fermi-edge terms reported in Eq. (32): Since at the absolute zero the Gaussian part is absent, we may readily interpret As shown in Figs. 4B-D, these results are in excellent agreement with the numerical calculations at β ω → ∞ by suitable adjustments of τ 0 . Such a regularization parameter depends on r F and takes the physical interpretation of the average time needed by the system to respond to the abrupt impurity perturbation at zero temperature (Fig. 4A). It is also interesting to note that if we let the harmonic frequency in Eq. (A.20) go to zero, by fixing α and keeping the number of particles in the gas (2r F ≈ ε F / ω) finite, we retrieve the Nozieres-De Dominicis result [11,12,14] Λ MND (t) = −α ln(it/τ 0 + 1), (A. 21) which leads to the propagator e Λ MND (t) = (it/τ 0 + 1) −α , originally calculated for a suddenly switched on core-hole in a free electron gas. Such a limiting procedure is illustrated in Fig. A1, where we clearly see that as ω → 0, the periodicity of the propagator e Λ ∞ 2P (t) vanishes and its modulus reduces to a single peak with tails going like 1/t. Inclusion of finite temperature corrections [14] may be done by considering all possible products of λ ∞ 0± (t) and λ ∞ m>0± (t), as given by Eq. (A.12), then using the expansion (A.9), and finally performing the t ′ and t ′′ -integrals excluding terms proportional to t 2 . However, the zero temperature result (A.20) is largely dominant in characteristic function χ β (t), and in the system response ν β (t), within the ranges of temperatures β ω = 0.4 − ∞ and particle numbers N 20 considered here.  Figure A1. Real and imaginary parts of the excited impurity propagator e Λ ∞ 2P (t) in the trapped (ω = 1) and free (ω → 0) fermion gas. The critical exponent is fixed to α = 0.1, while several regularisation times (τ 0 = 0.05, 0.1, 0.15) are tested. once the regularisation parameter τ 0 is adjusted in order to obey the conditions (33) and (34). The adjusted values of τ 0 have been also reported in Fig. 4A vs r F .
In Sec. 5 we have provided a method to determine the first three cumulants of the work distributions, see Eqs. (37)-(39). In particular, we have seen that these quantities depend on the auxiliary functions λ β ± (t) and their first-order time-derivatives at the time the impurity potential is activated. The first cumulant is well defined, being proportional to λ β + (0), and coincides with the first-order energy shift (A.6) discussed above, as reported in Eq. (43). On the other hand the second and third cumulants involve λ β − (0) and dλ β − (t)/dt| t=0 , see Eqs. (47) and (53). These two quantities contain weighted sums of the γ r and rγ r -factors by the hole occupation numbers f − r . For large r, the two series diverge like r −1/2 and r 1/2 , respectively. It is not surprising that a normalizable distribution function, with a well defined mean value, has divergent higher order moments, and this is due to both the sudden dynamics and the spatial modelling of the impurity potential. However, the regularization procedure introduced in Sec. 4.2, and applied above to derive Eq. (A.23), leads to a well defined characteristic function admitting finite cumulants of any order. To have a consistent theory, we also want the cumulants of the work distribution, as computed from Eq. (A.22), to correctly tend to the same values obtained from Eq. (A.23). With this in mind, we have proposed in Sec. 5 a renormalization procedure in which the auxiliary functions are extended on the complex-time domain t → t + iτ and shifted by a small imaginary time, i.e., τ = ±1/ω 0 in Eq. (50) and τ = ±1/ω ′ 0 in Eq. (54). Then, using Eqs. (A.4), (A.10) and (A.11) above, we can use the extended auxiliary functions λ ∞ ± (t + iτ ) to obtain the absolute zero expressions κ 2 (∞) = 2αε F ω λ ∞ + (0)λ ∞ − (0) = 2αε F √ πe −2ω(r F +1)/ω 0 √ 1 − e −2ω/ω 0F 2 1 (r F + 1, e −2ω/ω 0 ) − 2αε F e −4ω(r F +1)/ω 0F 2 1 (r F + 1, e −2ω/ω 0 ) 2 (A. 24) and To link the these expressions with the regularized ones that we have derived in Sec. 5, we have fixed the cut-off frequencies ω 0 and ω ′ 0 by constraining: The adjusted values of the corresponding cut-off times 1/ω 0 and 1/ω ′ 0 have been reported in Fig. 5A.