Work Statistics and Entanglement Across the Fermionic Superfluid‐Insulator Transition

Entanglement in many‐body systems may display quantum phase transition signatures, and analogous insights are emerging in the study of work fluctuations. Here, the fermionic superfluid‐to‐insulator transition (SIT) is considered and related to its entanglement properties and its work distribution statistics. Using the attractive fermionic Hubbard model with randomly distributed impurities, the work distribution is analyzed under two quench protocols triggering the SIT. In the first, the concentration of impurities is increased; in the second, the impurities' disorder strength is varied. The results indicate that, at criticality, the entanglement is minimized while the average work is maximized. This study demonstrates that, for this state, density fluctuations vanish at all orders, resulting in all central moments of the work probability distribution being precisely zero. For systems undergoing a precursor to the transition (short chains with finite impurity potential) numerical results confirm these predictions, with higher moments further from the ideal results. For both protocols, at criticality, the system absorbs the most energy with almost no penalty in terms of fluctuations: ultimately this feature can be used to implement a quantum critical battery. The impact of temperature on this critical behaviour is also investigated and shown to favor work extraction for high enough temperatures.


I. INTRODUCTION
The behavior of the work that can be extracted or absorbed by a many-body quantum system is an important question in quantum thermodynamics.In interacting many-body quantum systems, correlations that are inherently quantum may manifest in the form of entanglement, and may show interesting properties at criticality, including universal scaling [1].This has inspired a series of investigations on the work statistics across a quantum phase transition (QPT) [2][3][4][5][6].Many previous studies have focused on characterizing the scaling properties of the work statistics following a sudden quench across the QPT [4,[6][7][8][9].Recently, the effects of the finite-time dynamics have started to be addressed, with special attention to the statistics beyond the second moment [10][11][12][13].These studies opened the way for novel applications, such as heat engines implemented using many-body systems as working medium, [14][15][16][17][18][19], which can be regarded as an example of quantum advantage when exploiting quantum correlations [20].It has been demonstrated that interactions between particles may allow for a boost in their efficiency, which is much larger than that of a noninteracting system of the same size [20][21][22].In the particular case of a working medium undergoing a QPT, criticality may lead to supralinear scaling of power [17,23], however, even at the level of simple models, the practical implementation of these engines requires special control techniques to operate thermal cycles at finite-time, in addition to minimizing fluctuations of work and output power [24][25][26].
Entanglement has emerged as a powerful tool to detect and characterize QPTs [1,27,28].Pioneer works exploring bi-partite entanglement in critical models showed that QPTs may be associated to an entanglement ex-tremum or to the entanglement non-analytical behaviour [1,[29][30][31].Following these ideas, various measures of entanglement have been employed as witnesses of changes in phases of matter.In the context of localization in the presence of disorder, metals [32][33][34][35], bosonic systems [36][37][38][39], spinless [40,41] and spinfull [42][43][44] fermions were explored.The fermionic superfluid to insulator transition (SIT) [45][46][47] is one of the QPTs displaying interesting entanglement properties that has been least explored in applications to quantum thermodynamics.The SIT in disordered systems is a paradigmatic problem in quantum physics, as it spurs from the competition between free and coherent mobility, and localization of particles.Advances in experiments with optical lattices and ultracold atoms have made it possible to quantum-simulate the SIT with ultracold bosons [48] and fermions [49].For strongly correlated bosons in a Tonks-Girardeau gas, modeled by the Bose-Hubbard model, it has been shown that the superfluid-insulating transition is accompanied by a pinning QPT [50], which can be achieved at infinitesimally weak lattice potentials [51].In fermionic systems, superfluidity can emerge from attractive Coulomb interactions [52,53], while an insulating behavior is due to repulsion [54].In the presence of impurities, the interplay between attractive interactions and disorder leads to a richer phase diagram, in which an insulating, superfluid, and localized behavior can be present.In fermionic systems the conditions under which the SIT occurs are still under investigation.Refs.[55,56] investigated the average single-site entanglement in disordered chains with attractive Coulomb electron-electron interactions.These results indicated the absence of a critical potential V for the emergence of the SIT triggered by disorder, whilst the existence of a critical concentration C C and a critical particle density for which the entanglement is minimum.
These recent results could help fine tuning the parameters controlling the disorder landscape and attractive couplings to induce the fermionic SIT in, e.g., optical lattices.
Motivated by the prospects of exploiting the SIT in thermodynamic cycles [17] and by its interesting entanglement structure in fermionic systems [56], we compare the average work statistics of two types of sudden quenches of interest to the SIT.In the first protocol, the concentration of impurities is changed by adding an extra impurity to the initial state.In the second protocol, we consider the situation in which the impurity potential strength is instantaneously increased.For a system initially prepared at temperature T = 0, we show that if the initial ground-state has minimum entanglement, the average work will be maximum at the critical concentration C C , whilst its variance will be minimized at the same rate in which fluctuations in density-density correlations decay.We also inspect the skewness of the work distribution, which provides a measure of non-Gaussianity and has been shown to be very sensitive to phase transitions for systems driven in finite time [11,57].We observe a similar sensitivity for the SIT, with the skewness providing more information than the variance about the subtle interplay between impurity strength and Coulomb interaction in the second protocol.Finally, the effects of temperature are also discussed.As the temperature increases, features of the quantum phase transition are suppressed, and we observe that thermal fluctuations favor work extraction (instead of absorption) whenever thermal excitations become larger than the typical energy scales of the system.

II. MODELING THE SUPERFLUID-INSULATOR TRANSITION IN FERMIONIC SYSTEMS
The superfluid-insulator transition in fermionic systems can be described by the Hubbard model with attractive electron-electron interactions in the presence of disorder [56].In a 1D chain of L sites, the Hamiltonian reads ) where J > 0 is the hopping parameter, U < 0 is the attractive on-site Coulomb interaction and V i is the strength of the impurity potential at site i, the local disorder.The number operator is ni,σ = ĉ † i,σ ĉi,σ where ĉ † i,σ is the creation operator for a fermion with z-spin component σ =↑, ↓ in the i-th site, such that n = n ↑ + n ↓ and m = n ↑ − n ↓ are particle density and magnetization, respectively, with n σ = i Tr[n i,σ ρ] with ρ the state of the system.
For L → ∞ and in the absence of disorder (V i = 0), increasing attractive interaction drives the system from a superfluid state of weakly-coupled pairs (BCS-like) to a phase of tightly coupled dimers (Bose-Einstein limit).In contrast, by triggering the local disorder V i at a fixed interaction U , the pairs tend to localize and the system undergoes a superfluid to insulator transition.Localization can be reached in a superfluid by either increasing the disorder intensity or lowering the particle density [56].
It has been recently shown [56] that, by considering N i point-like impurities of same strength V i = V , randomly distributed along the L-size chain, the level of localization in the SIT depends on the concentration C = (N i /L) × 100%.For sufficiently strong disorder strength |V | |U |, it is possible to achieve full localization at the critical concentration C C = 100n/2.For |V | → ∞, this full localization is characterized by realspace pairs localization, marked by zero entanglement [56].
Here we consider the same disorder landscape − a certain concentration C of pointlike impurities randomly placed along the chain with disorder strength V −.Both long (L = 100) and short (L = 8) superfluid chains with m = 0 are considered.To avoid features related to specific configurations, for all quantities we perform an average over 100 samples for the long chain case, and over all possible impurity configurations for L = 8 chains.
For L = 100, we consider the ground-state average single-site entanglement, which has been successfully applied to explore quantum phase transitions in the Hubbard model [31, 42-44, 56, 58].We follow the same approach of Refs.[42][43][44]56] obtaining the ground-state single-site entanglement via the linear entropy.This is averaged first over the L sites, , where ρ i is the reduced density matrix of site i, and then over 100 samples of different impurity configurations compatible with each of the protocols, L = (1/M ) m L m .

A. Statistics of work in a sudden quench
One can perform (or extract) work on an isolate quantum system by changing the parameters g in the timedependent Hamiltonian Ĥ( g t ) taking the system outof-equilibrium.In this process, all possible transitions between the eigenestates of the initial Ĥ( g 0 ) and final Ĥ( g f ) Hamiltonians may be involved, determining the change in energy as well as its fluctuations.Experimentally, the quantum work is accessible by means of spectroscopic methods [59][60][61], and the work distribution has been reconstructed in a liquid-state nuclear magnetic resonance platform using small molecules as working fluid [61,62].
In a sudden quench, the transitions are instantaneous and the system does not have time to adapt.If the system is prepared in a superposition of the eigenstates |n 0 of Ĥ( g 0 ) with weights p n , the work probability distribu-tion P (W ) is calculated as follows [63] where { 0 n } and { 0 + m } are the eigenvalues of Ĥ( g 0 ) and Ĥ( g f ), respectively, and p m|n = | m 0 + |n 0 | 2 is the probability to find the system in the m-th eigenstate |m 0 + of Ĥ( g f ) at t > 0 + given the dynamics started in state |n 0 at t = 0.
Associated with the work distribution, we have the kth order central moments (W − W ) k which, for k > 1 can be calculated recursively, starting from the average work Average work performed on the system correspond to W > 0; while extracted work is signalled by W < 0.
The second moment, k = 2, is the variance, associated with the energy fluctuations.The third moment quantifies the skewness (k = 3) of the work distribution, which is related to the deviation from Gaussianity.

B. Sudden quench protocols across the SIT
We will study two sudden-quench protocols across the SIT and the associated average work.In the first protocol we vary the impurity concentration C for fixed values of V .For each initial concentration, we consider in turn each of the possible corresponding impurity configurations.For each initial impurity configuration, we consider in turn all possible configurations (L = 8) or 100 samples of different impurity configurations (L = 100) that can be achieved by adding one extra impurity, see Fig. 1(a).For each initial and final configuration, we calculate the average work and the moments of its distribution.For each initial concentration, we then average these quantities over all possible couples of initial and final configurations.
In the second protocol initial and final Hamiltonians will have the same number of impurities but their position may vary and the final impurity potential V f will be in modulus larger than the initial potential V i , see Fig. 1(b).The results will be averaged over all possible initial and final configuration with same impurity concentration C.
For both protocols, in the sudden quench limit, all central moments can be expressed in terms of correlation functions of the density, as the terms of the Hamiltonian associated with the kinetic energy and the Coulomb repulsion remains the same.Here we discuss in details the first three central moments.The average work is a functional of the local densities of the initial state: where ∆v j = (V f j − V 0 j ) The fluctuations of work, corresponding to the second central moment, depend on density-density fluctuations, as where we have used Finally, the third central moment, the skewness, can be expressed as the sum of densities' correlations where we have used In particular, for the ground-state |Ψ 0 equations 5, 9, and 13 are easily computed from the density correlations Ψ 0 |n j |Ψ 0 , Ψ 0 |n j n |Ψ 0 , and Ψ 0 |n j n ni |Ψ 0 .
The value of W , σ 2 W and µ 3 will vary depending on the initial (and final) impurity configurations of H 0 (of H f ).We will then consider their average over all possible configurations N c compatible with the relevant initial and final impurity concentrations (L = 8).For L = 100, N c will comprise 100 configurations.These averages are defined as While W is the average of the average work, for the sake of simplicity, in the rest of the paper, we will refer to it simply as the 'average work'.

A. The critical state
Vanishing entanglement at criticality reveals that, for each configuration, the ground-state is localized at the impurity sites, each of which is exactly doubly occupied in the limit |V | |U | and T = 0.In the occupation basis, this state is factorized and can be written as with j =↑↓ if site i contains an impurity and j = 0 otherwise.
For this type of localized state, density fluctuations vanish at all orders that is where i extends to any subset of sites.By inspection of eqs.(7), and (11) it is then clear that variance and skewness vanish exactly at criticality.As all higher central moments can be similarly written in term of density fluctuations, it follows from eq. ( 18) that all central moments of the work probability distribution exactly vanish at the critical state.

B. Quench in C: adding an extra impurity
We start our analysis with the average work -eqs.4 and 14-at T = 0 and the first protocol.
Previous results [56] for the average ground-state entanglement revealed that the emergence of the full localization in the SIT occurs at the critical concentration C C = 100 n/2 (for attractive disorder).As the disorder strength becomes comparable to the Coulomb attraction, the system starts to exhibit localization and, for |V | → ∞, it becomes fully localized as the entanglement vanishes.Figure 2(a) shows the average single-site entanglement L of the initial state [64] as a function of concentration for a chain of L = 100 sites, average density n = 0.8, U = −5J and various values of disorder strength.Figure 2(b) depicts W , for the same system undergoing a sudden quench that changes the concentration from C i to C i+1 .Our results show that W for the initial critical concentrations is maximum, corresponding to minimum entanglement.
This implies that the ground-state is localized at the impurity sites, and of a form very close to eq. (17).For this state, the only non-zero elements contributing to the sum in eq. 5, correspond to impurities which appear in both initial and final impurity configurations, hence providing up to N i /2 contributions of value |2V |.This implies that, in finite chains, the average work is amplified by increasing the disorder strength V , both because the value of each contribution is increased and because a stronger impurity increases localization for finite-chains (and finite temperature).With this in mind, by rescaling the average work by the disorder strength |V 0 |, we can check when full localization is achieved for this protocol and finite chains.This is explored in Figure 2 We now focus our attention on the statistics of the distribution and its changes with the temperature.Here, we will consider short chains (L=8) at half-filling n = 1 initially prepared in a thermal state.For T = 0J/k B and T = 2J/k B and all impurity potentials, work can be extracted from the system when starting from zero impurity as initial concentrations.This is a consequence of eq. 5 as it implies that V 0 i = 0 always, while V f i < 0 at the impurity site, that will result in negative work.However this type of contribution competes with positive contributions when starting from nonzero initial impurities concentrations, and work is to be performed on the system in this case -see Fig. 4(a)  and (d).
We note that eq. 5 is a consequence of considering a sudden quench dynamics, while finite-time dynamics would allow for redistribution of particle density in response to the new impurity distribution, as observedalbeit for repulsive impurities and repulsive particle interaction -in refs [65,66].This implies that it cannot be excluded that work could be extracted from the system at hand when considering finite-time dynamics.
At T = 0J/k B , both variance and skewness show a clear signature of criticality at concentration of 50 % see Fig. 4 (b) and (c).The variance almost vanishes, demonstrating, even in a short chain, (almost) no fluctuations for the system at criticality.This is consistent with eq. ( 9) for a system in the critical state eq.( 17), for which density fluctuations vanish at all orders.Indeed, our numerical results show that also the skewness vanishes at C C = 50%.In addition, it changes its sign from positive, at low concentration, to negative, after C C .A similar result was observed in Hubbard chains driven by a ramp field tuned at finite time [11].Results in [11] show this signature in the skewness to be stronger away from sudden quench and towards adiabaticity.In the present work, which used a sudden-quench dynamics, the change-in-sign signature in the skewness appears only for relatively high impurity potentials (T = 0J/k B ), and it disappears for T = 2J/k B (compare panels (c) and (f) of Fig. 4).At this intermediate temperature, the skewness preserves only a kink close to the critical concentration value and only at high impurity potentials, while the signature of criticality in the first and second moments is still well marked (Fig. 4 (d) and (e)), albeit only at relatively high impurity potentials for the variance.It would be interesting to check if a finite-time dynamics would restore a stronger signature in the third moment, as observed in [11].Panels (g), (h) and (i) of Fig. 4 confirms that any signature of critical behaviour is lost at high temperatures.In particular: work can now be extracted from the system at all concentration values and potentials, fluctuations become maximal around C C = 50%, and the skeweness has a smooth behaviour at all concentrations and impurity potentials.Here we consider quenching a different parameter, the impurities disorder strength.For each (fixed) impurity concentration and attractive Coulomb interaction of U = −5J, the impurities disorder strength will be quenched between an initial value V 0 and the final value V f = 2U .We consider initial disorder strengths from Let us first consider the results for T = 0 when plotted with respect to the varying parameter V 0 (Fig. 5, upper row).Comparing mean, variance and skewness, it is difficult to identify consistent signatures of a critical impurity strength.However, the mean (panel (a)) suggests a special role for V 0 = −5J, at which all curves with impurity concentrations C ≤ 50% cross.Indeed, the value −5J is the strength chosen for the Coulomb interaction, so it represents the watershed between either impurity strength or Coulomb attraction being the dominant interaction.Interestingly, a sudden quench in V allows for work extraction ( W < 0) for all concentrations when A trace of this critical role is present in the skewness (panel (b)), where curves close to the critical concentration cross at V 0 ≈ −5.However, the main feature of the skewness is a clear kink at V 0 = −1 for intermediate concentrations.In addition, for all concentration and small enough impurity potential the skewness changes sign.The variance presents some crossovers at small V 0 between curves corresponding to different concentrations, but no features at V 0 = U .
All moments confirm a special role for the data sets collected at the critical concentration C C = 50% (first row, green lines with diamond markers): the impurity concentration C is identified as a critical parameter even through this protocol.Hence, in the second row of Fig. 5, we plot the same data, but as a function of the impurity concentration C of each calculation.Now all three moments show a clear signature at the critical concentration, all developing a kink at C = 50% as V 0 increases.For large enough impurity strength, the variance almost vanishes confirming the quasi-absence of density-density The variance and skewness are symmetric with respect to the critical concentration, and present additional features at C = 25% and its symmetric C = 75%.The skewness kinks of panel (c) reflects in panel (f) with an initial potential V 0 < 1J being too small for the survival of any sign of criticality (panel (f), darkest blue curve), at least for this chain size.This is confirmed by the variance data (panel (e)).Sensitivity of the skewness to the special value V 0 = −5J = U translates into the extreme at critical concentration C = 50% changing from a minimum For increasing |V 0 |, this maximum tends towards zero skewness, as predicted in Sec.III A (not shown).In the mean, Increasing the temperature to T = 2J/k B has little effect in the average work extracted, compare Fig. 5(a) and (b) to Fig. 6(a) and (b).However, the second and third moments are more sensitive to the thermal fluctuations, with the signature corresponding to the critical concentration (kinks at C C in panels (e) and (f)) being washed out for V 0 < 3J.Thermal fluctuations also affect the behavior characterizing the watershed value V 0 = U in Fig. 5 (panels (a), (c), and (f)): in Fig. 6 these signatures are shifted at higher values of V 0 .Thermal fluctuation increase the probability of higher energy transitions, and indeed the work distribution becomes wider (second moment), but also more asymmetric (third moment): compare ranges of the y-axis in, e.g., panels (e) and panels (f) of Figs. 5 and 6.
Finally, high-temperature results (T = 30J/k B , Fig. 7) washes away any sign of critical behaviour, as expected.Interestingly, for this lesser quantum-correlated system, work can now be extracted for any parameter combination (i.e.ˆ W < 0), suggesting that work has to be done to the system in order to create quantum correlations.perature.We stress that here we do not averaged over all possible configurations, but just plot results for a single calculation.Even so, it can be observed how, with increasing temperature, the distributions at any C i become wider and then, at T = 30, more regular (single relevant maximum ) though still displaying a substantial skewness.

IV. CONCLUSION
In this paper, we investigated the statistics of work in a fermionic quantum system undergoing a sudden quench across the superfluid-insulator transition.We demonstrate that at the critical concentration all central moments of the work probability distribution exactly vanish.In fact, in a sudden quench, the work is reduced to a sum of local operators, while the SIT critical state has zero average-site entanglement.This implies that all density fluctuations exactly vanish at all orders and hence all central moments of the work probability distribution follow suit.We note that this property would extend to quenching protocols in other phases for which local entanglement is negligible.
We numerically compared two paths for the SIT in fermionic systems with randomly distributed impurities and explored the first three moments of the work distri-bution from very low to high temperatures.Our results indicate that the SIT increases work absorption both when triggered by a change in the impurity concentration or by a quench in the disorder strength.Indeed, we demonstrate that the average work is maximized when, at criticality, the single-site entanglement is minimized and explain this analytically.This property could be employed for the implementation of efficient energy storage (quantum batteries).
Work fluctuations become more pronounced away from the critical concentration.Interestingly, deviations from Gaussianity, as revealed by the skewness of the distribution, behave qualitatively differently in the two protocols.When triggering the SIT by varying the impurity concentration the skewness changes sign at the critical point, resembling the behaviour seen in [11]; however the skewness remains negative on both sides of the critical concentration when the transition is triggered by varying the impurities' potential strength.In this second protocol, both average work and skewness are also highly sensitive to the interplay between impurity and Coulomb potential, with the regime in which the impurity potential dominates been clearly marked in both moments.We verified that all signatures of the SIT are washed by increasing temperatures.
Future directions include the development of finitetime protocols in which the trade-off between the work and its fluctuations can be controlled dynamically, such as proposed in [67], and its application to improve the performance of fermionic thermal machines operating at finite temperature.

FIG. 1 .
FIG. 1. Illustration of sudden quench protocols across the superfluid-insulator transition in a chain of L = 8 sites.Panel (a) shows a quench in the concentration of impurities and a few possible configurations starting from one of the 28 possible configurations with Ni = 2 impurities.For 3 impurities (final state), there are a total of 56 configurations and we show only 8 of them.Panel (b) shows a possible initial state for a quench in the disorder strength at fixed concentration with all possible final configurations with Ni = 1 impurity.
(c), which demonstrates that the rescaled work does not change (full localization for this protocol) for |V 0 | > |U |.

Fig. 3
shows examples of P (W ) at different concentrations and for the three chosen temperatures, T = 0J/k B (panel (a)), T = 2J/k B (panel (b)), and T = 30J/k B (panel (c)).The intermediate temperature has been chosen to be low enough so that the most populated state is the ground-state, while a few low-lying states contribute to the dynamics.Results for the moments of the work distribution are shown in Fig. 4, for T = 0J/k B , T = 2J/k B and T = 30J/k B , as indicated.Here the moments are averaged over all configurations.

FIG. 2 .
FIG. 2. Panel (a): Single-site entanglement averaged over 100 of the possible impurity configurations versus concentration.Panel (b): W /J, averaged over 100 of the possible impurity configurations, versus initial concentration Panel (c): Average work from panel (b) scaled by the disorder strength |V0|.

FIG. 3 .
FIG. 3. Example of work distributions resulting from quenches in C for chains with L = 8 sites, with fixed initial and final impurity potential V0 = V f = −5J at T = 0 (a), T = 2J/kB (b) and T = 30J/kB (c).The initial and final configurations are randomly picked among all possibilities for Ci and Ci+1.However, the chosen initial and final configurations' pairs remain the same at different T 's.

|V 0 |
< |U | and a decreasing range of concentrations as V 0 becomes dominant over the Coulomb interaction.

FIG. 4 .
FIG. 4. Average quantum work (first column), variance (second column) and skewness (third column) versus impurities' concentration of the initial state for systems with L = 8 sites, particle interaction U = −5J and for different disorder intensities.The impurities' concentration of the final state is always the next concentration.Temperature varies as follows: first row: T = 0J/kB; second row: T = 2J/kB; third row: T = 30J/kB;

FIG. 5 .
FIG. 5. First row: Average quantum work (left), variance (middle) and skewness (right) versus V0 (the impurities' strength of the initial state) for different disorder concentrations, as labeled.The impurities' strength of the final state is always V f = −10J.Other parameters are: L = 8, U = −5J, and T = 0J/kB.Second row: Same data as for first row but plotted versus the impurities' concentration of the systems considered.

Figure 8
Figure 8 displays one example of P (W ) for each tem-

FIG. 8 .
FIG. 8. Work distribution for chains with L = 8 sites resulting from a quench in V from V0 = −J to V f = −10J with fixed concentration Ci at T = 0J/kB, T = 2J/kB and T = 30J/kB.The initial and final configurations are randomly picked among all possibilities for Ci.However, initial and final configurations for the same Ci at different T are the same.