Quench action and large deviations: work statistics in the one-dimensional Bose gas

TL;DR

This paper analyzes the large deviations of work in a one-dimensional Bose gas after an interaction quench, revealing non-Gaussian fluctuations near the mean and algebraic decay for larger deviations, using the quench action and Bethe ansatz.

Contribution

It provides an exact computation of the rate function for work fluctuations in an interacting Bose gas post-quench, extending large deviation analysis to quantum integrable systems.

Findings

Fluctuations near the mean are non-Gaussian.

Large deviations below the mean are exponentially suppressed.

Deviations above the mean decay algebraically with an unknown exponent.

Abstract

We study the statistics of large deviations of the intensive work done in an interaction quench of a one-dimensional Bose gas with a large number N of particles, system size L and fixed density. We consider the case in which the system is initially prepared in the non-interacting ground state and a repulsive interaction is suddenly turned on. For large deviations of the work below its mean value, we show that the large deviation principle holds by means of the quench action approach. Using the latter, we compute exactly the so-called rate function, and study its properties analytically. In particular, we find that fluctuations close to the mean value of the work exhibit a marked non- Gaussian behavior, even though their probability is always exponentially suppressed below it as L increases. Deviations larger than the mean value, instead, exhibit an algebraic decay, whose exponent can not be determined directly by large-deviation theory. Exploiting the exact Bethe ansatz representation of the eigenstates of the Hamiltonian, we calculate this exponent for vanishing particle density. Our approach can be straightforwardly generalized to quantum quenches in other interacting integrable systems.

Equations

$$H(c)=-\frac{{\hslash }^2}{2m}\sum _{j=1}^N\frac{{\partial }^2}{\partial x_j^2}+2c\sum _{j<k}\delta (x_j-x_k).$$

$${\psi }_N\left(x_1,\dots ,x_N\right)$$

$$e^{-i{\lambda }_{j}L}=\prod _{k\ne j}^N\frac{{\lambda }_k-{\lambda }_j+ic}{{\lambda }_k-{\lambda }_j-ic},\text{mbox}withj=1,\dots ,N.$$

$${\lambda }_j=\frac{2\pi I_j}{L}-\frac{2}{L}\sum _{k=1}^N\arctan \left(\frac{{\lambda }_j-{\lambda }_k}{c}\right),$$

$$E\left[\{{\lambda }_j{\}}_{j=1}^N\right]=\sum _{j=1}^N{\lambda }_j^2.$$

$${\rho }^t(\lambda )=\frac{1}{2\pi }+\frac{1}{2\pi }{\int }_{-\infty }^{\infty }\mathrm{d}\mu K(\lambda -\mu )\rho (\mu ),$$

$${\rho }^t(\lambda )=\rho (\lambda )+{\rho }^h(\lambda ),$$

$$K(\lambda )=\frac{2c}{{\lambda }^2+c^2}.$$

$$\eta (\lambda )=\frac{{\rho }^h(\lambda )}{\rho (\lambda )}.$$

$$D[\rho ]=\underset{N,L\to \infty }{\lim }\frac{N}{L}={\int }_{-\infty }^{+\infty }\mathrm{d}\lambda \rho (\lambda ),$$

$$e[\rho ]=\underset{N,L\to \infty }{\lim }\frac{E}{L}={\int }_{-\infty }^{+\infty }\mathrm{d}\lambda \rho (\lambda ){\lambda }^2.$$

$$\mathcal{Z}(\beta )=\mathrm{tr}\left[e^{-\beta H}\right]=\sum _{\{{\lambda }_j\}}{e}^{-\beta E\left[\{{\lambda }_j\}\right]}.$$

$$\mathcal{Z}(\beta )=\int \mathcal{D}\rho e^{-L{S}_{\mathrm{th}}[\beta ,\rho ]},$$

$$S_{\mathrm{th}}[\beta ,\rho ]$$

$$s_{YY}[\rho ](\lambda )$$

$$f=-\frac{T\ln \mathcal{Z}}{L}=Dh-\frac{T}{2\pi }{\int }_{-\infty }^{\infty }\mathrm{d}\lambda \ln \left(1+{\mathrm{e}}^{-\epsilon (\lambda )/T}\right).$$

$$\epsilon (\lambda )$$

$${\psi }_N^{(0)}(x_1,x_2...x_N)=\frac{1}{L^{N/2}}.$$

$$\{{\lambda }_j{\}}_{j=1}^N=\{{\lambda }_j^{+}{\}}_{j=1}^{N/2}\cup \{-{\lambda }_j^{+}{\}}_{j=1}^{N/2},$$

$$\{{\lambda }_j{\}}_{j=1}^N=\{{\lambda }_j^{+}{\}}_{j=1}^{(N-1)/2}\cup \{-{\lambda }_j^{+}{\}}_{j=1}^{(N-1)/2}\cup \{0\},$$

$$⟨\{{\lambda }_j\}|\mathrm{BEC}⟩=\sqrt{\frac{(cL{)}^{-N}N!}{{\det }_{j,k=1}^N{G}_{jk}}}\frac{{\det }_{j,k=1}^{N/2}{G}_{jk}^Q}{\prod _{j=1}^{N/2}\frac{{\lambda }_j}{c}\sqrt{\frac{{\lambda }_j^2}{c^2}+\frac{1}{4}}}.$$

$$G_{jk}$$

$$G_{jk}^Q$$

$$⟨{\psi }_0|\mathcal{O}(t)|{\psi }_0⟩=\sum _{n,m}⟨n|\mathcal{O}|m⟩⟨{\psi }_0|n⟩⟨m|{\psi }_0⟩e^{i(E_n-E_m)t},$$

$$\underset{t\to \infty }{\lim }⟨{\psi }_0|\mathcal{O}(t)|{\psi }_0⟩=\sum _n|⟨n|{\psi }_0⟩{|}^2⟨n|\mathcal{O}|n⟩$$

$$=\sum _{\{{\lambda }_j\}}|⟨\{{\lambda }_j\}|{\psi }_0⟩{|}^2⟨\{{\lambda }_j\}|\mathcal{O}|\{{\lambda }_j\}⟩,$$

$$\sum _{\{{\lambda }_j\}}|⟨\{{\lambda }_j\}|{\psi }_0⟩{|}^2$$

$$\underset{t\to \infty }{\lim }⟨{\psi }_0|\mathcal{O}(t)|{\psi }_0⟩=\int \mathcal{D}\rho ⟨\rho |\mathcal{O}|\rho ⟩e^{-L{S}_{\mathrm{QA}}[\rho ]},$$

$$S_{\mathrm{QA}}[\rho ]$$

$$|⟨\{{\lambda }_j\}|{\psi }_0⟩|\simeq e^{-L{S}_O[\rho ]}.$$

$$S_O[\rho ]$$

$$\gamma =\frac{c}{D}.$$

Full text

I Introduction

A defining feature of the quantum theory is the presence of statistical fluctuations in the outcome of the measurement of any physical observable: given a collection of systems prepared in the same quantum state, the result of identical measurements will be generally different, leading to a non-trivial probability distribution function for the outcomes. For many-body systems, it is often very difficult to determine the latter and to give predictions beyond the corresponding mean value. This is especially true out of equilibrium, where exceptional computational challenges arise, even in prototypical solvable models CaEM16 . On the other hand, old atomic experiments have recently proven that the full probability distribution of certain observables can be probed in mesoscopic systems HLSI08 ; AJKB10 ; JABK11 ; KPIS10 ; KISD11 , motivating new efforts in the theoretical study of the fluctuations of quantum measurements in many-body physics, both in ChDe07 ; LaFe08 ; ia-13 ; sk-13 ; e-13 ; k-14 ; mcsc-15 ; sp-17 ; CoEG17 ; nr-17 ; hb-17 and out EiRW03 ; gadp-06 ; er-13 ; lddz-15 ; GrEC18 ; bpc-18 ; BaPi18 ; CoEs19 of equilibrium.

In this respect, a natural quantity which can be investigated is the work done upon changing some of the system’s parameters Jarz97 ; TalLutz07 ; Silv08 ; GooldPlastina18 , as, e.g., in a quantum quench cc-06 ; PolkSeng11 . Crucially, its statistics has now proven to encode important information on the internal dynamics of the system Silv08 ; PaSi09 ; SmSi13 ; MaSi14 ; PaSo14 ; Palm15 ; Rotondo18 and to display interesting features such as, most prominently, an emergent universal behavior BDKP11 , in quenches near a critical point GaSi12 ; GaSi11 ; GaSi11 ; GaSi12 ; SoGS13 . Furthermore, the statistics of the work is a valuable tool for studying dynamical phase transitions HePK13 ; Heyl14 ; Heyl18 and for detecting them HePK13 ; Heyl14 ; AbKe16 .

Despite its importance, the explicit calculation of the statistics of the work remain a difficult task even in simplified protocols such as quenches cc-06 ; PolkSeng11 and analytic results in the presence of interactions exist only in a few special cases PaSo14 ; Palm15 ; RyAn18 . Yet, some of its features can be understood based on general arguments and on the analysis of non-interacting models Silv08 . For instance, several studies have shown the existence of a universal edge singularity at the lowest threshold of the probability distribution $P(W)$ of the extensive work $W$ Silv08 ; SoGS13 , and its robustness against different non-equilibrium protocols SmSi13 . These features have been also verified in integrable quantum field theories PaSo14 ; Palm15 ; RyAn18 , by exploiting the fact that $P(W)$ can be obtained as the Fourier transform of the so-called post-quench Loschmidt amplitude Silv08 ; HePK13 .

In addition to $W$, one could also study the probability distribution of the intensive work $w=W/L$, where $L$ is the system size. As $W$ is an extensive variable, one generically expects that both its mean value $\langle W\rangle$ and its fluctuations $\langle(\Delta W)^{2}\rangle=\langle W^{2}\rangle-\langle W\rangle^{2}$ grow proportionally to $L$ as $L$ increases, the latter statement being true whenever $W$ can be seen as resulting from the sum of a number proportional to $L$ of almost independent contributions. This implies that the typical fluctuations in the value of $W$ are of order $\sqrt{L}$, i.e., those in $w$, of order $1/\sqrt{L}$, vanish as $L\rightarrow\infty$. Correspondingly, the distribution function $p(w)$ of $w$ approaches a delta function $\delta(w-\bar{w})$ which selects the average value $\bar{w}=\langle W\rangle/L$. On the other hand, fluctuations of the value of $w$ away from $\bar{w}$ corresponds to fluctuations of order $L$ in $W$, i.e., to large and atypical fluctuations, which are increasingly rare as $L$ grows. For free Bosonic and Fermionic models it was found GaSi12 that $p(w)\sim{\rm exp}[-LI(w)]$ where $I(w)$ is the so-called rate function of large deviation theory Touc09 , a non-negative function which vanishes for $w=\bar{w}$ and which controls the rate of exponential suppression of large deviations. Importantly, it was shown in Ref. GaSi12 that $I(w)$ provides insight into the universal properties of the system for $w\ll\bar{w}$. Furthermore, it was argued that its qualitative behavior can be inferred based on the knowledge of a few parameters of the quench. The analysis of Ref. GaSi12 also revealed that, for free bosonic models starting from a critical initial state, a further universal behavior appears in the regime $w>\bar{w}$, where $p(w)$ displays a transition from the aforementioned exponential decay to an algebraic decay. This transition is analogous to well-known phenomenon of Bose-Einstein condensation in quantum statistical mechanics Huang .

In the presence of interactions, the picture presented in Ref. GaSi12 remains qualitatively correct but the rate function $I(w)$ is very hard to compute in practice. However, accounting for interactions is obviously important, for instance for quantitative comparisons between theoretical predictions and experiments, in which the interactions are usually non-negligible. In this paper, we show that the statistics of the large deviations of the work done during a global quench can be determined quite generically for any interacting integrable system by means of the recently introduced quench action method CaEs13 ; DWBC14 ; Caux16 . In order to exemplify this approach, we study a prototypical example, namely an interaction quench in the Lieb-Liniger model LiLi63 ; Lieb63 describing a one-dimensional gas of $N$ Bosons with point-wise repulsive interaction. Specifically, we consider the protocol DWBC14 where the system, initially prepared in its non-interacting ground state, is brought out of equilibrium by suddenly turning on a repulsive interaction between the particles.

The aim of this paper is threefold. First, by means of a case study, we show that the quench action method predicts a statistics $p(w)$ of the large deviations of the intensive variable $w$ which takes the exponential form $p(w)\sim{\rm exp}[-LI(w)]$ for $L\rightarrow\infty$ at fixed density $D=N/L$, i.e., $p(w)$ naturally satisfies the so-called large deviation principle Touc09 . In addition, the quench action method allows the calculation of $I(w)$. Second, from this result we carry out a quantitative analysis of the interaction quench described above, pointing out its most interesting features (strongly depending on the presence of the interactions) which cannot be captured by the qualitative picture presented in Ref. GaSi12 . Third, going beyond large deviation theory, we analyze the region $w>\bar{w}$ where $I(w)$ vanishes identically and $p(w)$ has an algebraic decay upon increasing $w$. Although, in this case, the quench action method is not sufficient to quantitatively describe $p(w)$, we are able to determine the exponent of this algebraic tail by performing a finite-size calculation in the limit of vanishing densities $D$ of the Bosons. To the best of our knowledge, this provides the first quantitative description of the “condensed regime” characterizing the work statistics of interacting Bosonic systems for quenches starting from a critical initial state. Our results are expected to be relevant for experimental realizations of one-dimensional Bose gases in cold-atomic settings KiWW04 ; KiWW06 ; AEWK08 ; HGMD09 ; EWAR10 ; KHML10 ; FCFF11 ; DBAD12 ; FPCF15 and for measurements of the work statistics in non-equilibrium protocols HSKDL08 ; DCHFGV13 ; MDCP13 ; CBKZH13 ; Batalhaoetal ; Anetal .

The rest of this paper is organized as follows. In Sec. II we introduce the Lieb-Liniger Hamiltonian and its solution via the Bethe ansatz (Subsecs. II.1 and II.2), together with details on the quench protocol and the quench action approach (Subsecs. II.3 and II.4). In Sec. III, we recall known results on the statistics of the work done in a quantum quench. In Sec. IV we show that the quench action method naturally leads to a probability distribution $p(w)$ of the intensive work which takes the typical form of the large deviation principle, i.e., $p(w)\sim{\rm exp}[-LI(w)]$, where $I(w)$ can be computed within this formalism. Section V contains all the results regarding the rate function obtained via the quench action approach, while in Sec. VI we determine the power-law decay characterizing the fluctuations exceeding the mean value $\bar{w}$ in the limit of vanishing densities of the Bosons. Finally, we report our conclusions in Sec. VII, while the most technical aspects of our work are presented in several appendices.

II The model and the quench protocol

II.1 The Hamiltonian

We consider the Lieb-Liniger model Lieb63 ; LiLi63 describing a gas of $N$ Bosons at positions $\{x_{1},...,x_{N}\}$ with mass $m$ and point-wise repulsive interactions, with Hamiltonian

[TABLE]

The interaction strength $c$ is related to the scattering length $a_{\rm 1D}$ in one dimension through $c=-\hbar^{2}/ma_{\rm 1D}$ olshanii-98 and it can be varied via Feshbach resonances iasm-98 . In the following we set $\hbar=2m=1$ and assume that the $N$ Bosons are confined within a one-dimensional ring of length $L$, realizing periodic boundary conditions.

The Hamiltonian in Eq. (1) is integrable and, as a consequence, it can be diagonalized exactly by means of the Bethe ansatz Lieb63 . In particular, the $N$-particle eigenfunctions can be written as

[TABLE]

where the sum runs over the $N!$ permutations $P$ of $N$ elements. The parameters $\{\lambda_{j}\}_{j=1}^{N}$ are the so-called rapidities and, in analogy with the quasi-momenta which are relevant in the case of free quantum gases, they parametrize the different eigenstates of the Hamiltonian. When the system has a finite extension $L$, the rapidities have to satisfy a set of quantization conditions which are known as Bethe equations LiLi63

[TABLE]

For repulsive interactions $c>0$, it can be shown that all $\lambda_{j}$’s are real; accordingly, it is convenient to consider the logarithm of Eq. (3), i.e.,

[TABLE]

where we introduced the quantum numbers $I_{j}$. These numbers parametrize the sets of rapidities $\{\lambda_{j}\}$ and are integers (half-integers) for odd (even) $N$; note that they have to be chosen in such a way that $I_{j}\neq I_{k}$ for $j\neq k$ korepin_book . The knowledge of the rapidities $\{\lambda_{j}\}_{j=1}^{N}$ completely specifies the eigenstates $\psi_{N}$ of $H(c)$ and their properties. For example, the corresponding energy eigenvalue can be written as

[TABLE]

In the following, we will denote the normalized eigenstate of $H(c)$ corresponding to a set of rapidities $\{\lambda_{j}\}$ as $|\{\lambda_{j}\}\rangle$.

II.2 The thermodynamic description

When the number $N$ of particles is very large, the explicit form of the wave function in Eq. (2) becomes difficult to deal with, and the Bethe equations (4) harder to solve numerically. For these reasons, in order to study the thermodynamic limit of the model, it is necessary to employ an appropriate “thermodynamic Bethe ansatz” formalism YaYa69 , which we briefly review in this section. Here we only report the aspects that are directly relevant to the present study, while the interested reader is referred to Ref. takahashi_book for a thorough treatment.

In the thermodynamic limit $L,N\to\infty$ with fixed density $D=N/L$, it can be seen that the rapidities $\lambda_{j}$ of a given eigenstate, whose number also grows to infinity, arrange themselves on the real line according to a smooth distribution function $\rho(\lambda)$, with $\lambda\in(-\infty,+\infty)$. Complementary to the latter, one can also introduce a distribution of “holes” $\rho^{h}(\lambda)$, namely of unoccupied states. The functions $\rho(\lambda)$ and $\rho^{h}(\lambda)$ are analogous to the distributions of momenta and vacancies for free Fermi gases. However, contrary to the non-interacting case, $\rho$ and $\rho^{h}$ are related in a non-trivial way. In particular, they satisfy the following integral equation

[TABLE]

where we defined the total distribution function

[TABLE]

and the kernel

[TABLE]

Equation (6) can actually be derived by taking the thermodynamic limit of the Bethe equations (4) takahashi_book . For future use, we also introduce the following standard definition

[TABLE]

It is widely believed that the knowledge of the rapidity distribution function $\rho(\lambda)$ is sufficient to compute all of the thermodynamic properties of the corresponding eigenstate. For example, in the thermodynamic limit, the densities $D[\rho]$ and $e[\rho]$ of particles and energy per unit length can be obtained, respectively, as

[TABLE]

One of the advantages of the thermodynamic description introduced above is the possibility to replace discrete sums over eigenstates with functional integrals over rapidity distribution functions. This is best illustrated by the computation of the thermal partition function at temperature $T=\beta^{-1}$ with $k_{B}=1$, i.e.,

[TABLE]

Note that, while each term on the right-hand side is known, the sum runs over all the possible sets of rapidities $\{\lambda_{j}\}$ and hence it is very difficult to evaluate in practice for large $N$. In the thermodynamic limit, however, one can rewrite Eq. (12) as a functional integral takahashi_book

[TABLE]

where the functional

[TABLE]

plays the role of a thermal free energy. Here the first term on the right-hand side is derived using Eq. (11) for the thermodynamic limit of the energy in Eq. (5) and corresponds to the exponential in Eq. (12). The second term, instead, is the so-called Yang-Yang entropy YaYa69 $S_{\rm YY}=\int{\rm d}\lambda\,s_{YY}(\lambda)$ with density

[TABLE]

which accounts for the fact that each rapidity distribution function $\rho(\lambda)$ emerges from several “microscopic realizations”, i.e., that there are many sets of rapidities $\{\lambda_{j}\}$ associated with the same function $\rho(\lambda)$ takahashi_book . For $L\rightarrow\infty$ the functional integral in Eq. (13) can be computed by a saddle-point evaluation, yielding the following expression for the free energy density $f$ associated with the thermal partition function in Eq. (12) takahashi_book ; korepin_book

[TABLE]

Here $\varepsilon(\lambda)$ is the solution to the integral equation

[TABLE]

and $h$ is a Lagrange multiplier introduced in order to enforce the assigned density of particles, according to Eq. (10). As we will see in Sec. II.4, the formalism discussed above will also be essential for introducing the quench action approach.

II.3 The quench protocol

As anticipated, in this work we consider an interaction quench in which the system is initially prepared in the ground state of the non-interacting Hamiltonian $H(c_{0}=0)$, usually denoted by $|\rm{BEC}\rangle$. The corresponding wave-function $\psi_{N}^{(0)}(x_{1},x_{2}\ldots,x_{N})=\langle x_{1},x_{2},\ldots\,,x_{N}|\rm{BEC}\rangle$ reads

[TABLE]

At time $t=0$, a finite inter-particle repulsive interaction $c>0$ is turned on and the gas is subsequently left to evolve unitarily. Our motivation to investigate the above quench is twofold: first, the simplicity of the initial state allows one to derive analytic predictions which would be difficult to obtain in general; second, we will see that this kind of quench leads to interesting features in the work statistics.

The non-equilibrium dynamics arising from the interaction quench $c_{0}\rightarrow c$ described above has been extensively investigated in the literature GrRD10 ; KSCC13 ; KoCC14 ; DeCa14 ; DePC15 ; PiCE16 ; ZWKG15 . From the analytical point of view, an important result has been the discovery in Ref. DWBC14 of an exact formula (later proven in Ref. Broc14 ) for the overlaps between the initial state in Eq. (18) and the eigenstates of the Hamiltonian in Eq. (2). This formula, which is an essential ingredient for the application of the quench action approach described in Sec. II.4, will be used several times in this work, and is hence reviewed in what follows.

It was first shown in Refs. DWBC14 ; Broc14 that the initial state in Eq. (18) has a non-vanishing overlap only with eigenstates corresponding to sets $\{\lambda_{j}\}$ of rapidities which are parity invariant, i.e., such that $\{\lambda_{j}\}=\{-\lambda_{j}\}$. This implies that the set $\{\lambda_{j}\}$ can be decomposed as

[TABLE]

if $N$ is even and

[TABLE]

if $N$ is odd, where $\lambda^{+}_{j}>0$. For these states, the overlap formula is extremely simple. Explicitly, for even $N$, it reads

[TABLE]

Here we introduced the matrices $G_{jk}$ and $G_{jk}^{Q}$, with elements

[TABLE]

where $K^{Q}(\lambda,\mu)=K(\lambda-\mu)+K(\lambda+\mu)$, and $K(\lambda)$ is defined in Eq. (8). An analogous result holds for the case of odd $N$ BDWC14 . As we will see further below, our analytic study ultimately hinges on the existence of the exact formula (21).

II.4 The quench action method

In this section we discuss the quench action approach, which is the last piece of technical background needed in order to carry out our analysis of the statistics of the work. In the following, we only review some relevant aspects, referring the reader to the literature for a more comprehensive treatment Caux16 .

This integrability-based method has been introduced in Ref. CaEs13 to tackle the difficult problem of computing the thermodynamic limit of time averages

[TABLE]

after a global quench, where $|\psi_{0}\rangle$ is the initial state and $\mathcal{O}$ a generic local observable. Here we denoted by $|n\rangle$ and $|m\rangle$ the eigenstates of the Hamiltonian with energies $E_{n}$ and $E_{m}$, respectively. While, in principle, the quench action approach can be used to compute the full evolution of the expectation value in Eq. (24) DeCa14 ; DePC15 ; BeSE14 ; PiCa17 , it is particularly effective if one is interested only in its infinite-time limit Caux16

[TABLE]

where we used the fact that for the Lieb-Liniger model the eigenstates are parametrized by sets of rapidities $\{\lambda_{j}\}$. Indeed, the quench action approach provides a simple prescription to evaluate the thermodynamic limit of Eq. (25), which is based on replacing the spectral sum

[TABLE]

in Eq. (25) with a functional integration, in analogy with what we did in Eq. (13). By doing so, one arrives at the formal expression

[TABLE]

where we denoted by $|\rho\rangle$ an eigenstate whose rapidities $\{\lambda_{j}\}$ approach the distribution $\rho(\lambda)$ in the thermodynamic limit. Here, the functional $S_{\rm QA}[\rho]$, usually called the quench action (QA), plays a role analogous to the thermal free energy in Eq. (14). Explicitly, it reads Caux16

[TABLE]

where $S_{O}[\rho]$ is the functional associated with the overlap term in the spectral sum in Eq. (26), i.e.,

[TABLE]

Note that, since the Hamiltonian in Eq. (1) conserves the particle number which is well-defined in the initial state, a Lagrange multiplier $h$ has been introduced in Eq. (28) (where the prefactor $1/2$ is for later convenience): this allows us to extend the functional integration over the whole space of rapidity distribution functions. Note also that the Yang-Yang entropy appearing in Eq. (28) bears an additional prefactor $1/2$, which is due to the fact that only parity-invariant eigenstates contribute to the spectral sum in Eq. (26) DWBC14 . The integral in Eq. (27) can now be computed via the saddle-point method, yielding an exact result for the infinite-time average of the expectation value of local observables in the thermodynamic limit.

A crucial point in the procedure outlined above is the availability of an analytic expression for the functional $S_{O}[\rho]$ in Eq. (29). If an analytic expression of the form (21) is known for the overlaps, that for $S_{O}[\rho]$ can be easily derived as shown explicitly in Ref. DWBC14 ; in particular, for the case of the initial state in Eq. (18) one has

[TABLE]

where we introduced the normalized interaction strength

[TABLE]

Unfortunately, for arbitrary initial states it remains an open problem whether formulas analogous to Eq. (21) can be derived, so that, in general, the explicit expression for $S_{O}[\rho]$ is unknown KoPo12 ; Pozs14 ; PiCa14 ; BNWC14 ; LeKZ15 ; LeKM16 ; Pozs18 ; LeKL18 . Note, however, that it was recently shown that for any integrable model it is always possible to find a class of “integrable initial states” for which this can be done PiPV17 ; PoPV18 . As we will comment on later, the results derived in this work can thus be generalized straightforwardly to other integrable systems, at least for quenches from the latter class of initial states.

III The statistics of the work

Before presenting our results, we review some generic features of the work statistics obtained in free models TalLutz07 ; Silv08 ; GaSi11 ; GaSi12 ; SoGS13 ; SmSi13 . The following discussion will be useful for a comparison with the interacting case analyzed in this paper.

Consider a quantum system with $N$ degrees of freedom, initially in the ground state $|\psi_{0}\rangle$ of its Hamiltonian $H(c_{0})$. In the following, we denote its ground-state energy eigenvalue with $E_{0}^{c_{0}}$. The probability distribution $P(W)$ of the extensive work $W$ done in quenching a global parameter $c_{0}\rightarrow c$ is defined as TalLutz07

[TABLE]

where $|\psi_{n}^{c}\rangle$ are the eigenstates of the post-quench Hamiltonian $H(c)$ with corresponding energies $E_{n}^{c}$. One immediately notices from Eq. (32) that the work has a minimum threshold value $W_{\rm{rev}}=E_{0}^{c}-E_{0}^{c_{0}}$. This has the meaning of reversible work, i.e., the work performed at zero temperature when the transformation $c_{0}\rightarrow c$ is done in a reversible way. As a consequence, we will refer the work to this threshold, focussing on the irreversible contribution $W_{\rm{irr}}=W-W_{\rm{rev}}\geq 0$, which is related to the irreversible entropy production PolkSeng11 ; DorGoold12 . For convenience, we will henceforth indicate $W_{\rm{irr}}$ by $W$, dropping the subscript.

In the following, we will be interested in the moment generating function $G(s)$ of the extensive work $W>0$

[TABLE]

and in the corresponding scaled cumulant generating function (SCGF) $f(s)$, defined by Touc09

[TABLE]

Note that, from Eq. (33), $\langle W\rangle=\langle\psi_{0}|H(c)|\psi_{0}\rangle-E_{0}^{c}\,,$ which motivates the expectation that $\langle W\rangle\propto N$. As anticipated in the introduction, in order to investigate the large deviations of the random variable $W$, it is convenient to focus on the intensive irreversible work $w=W/N$, with probability density $p(w)$. Upon increasing $N$ one generically expects that $p(w)$ satisfies the so-called large deviation principle Touc09 , i.e., that $p(w)\sim{\rm exp}[-NI(w)]$, where $I(w)$ is referred to as the rate function. This function is non-negative, convex and, in general, displays a unique zero at the average and most probable value $\bar{w}$: upon increasing $N$ the function $p(w)$ becomes peaked around $\bar{w}$, being exponentially suppressed for $w\neq\bar{w}$. Importantly, $I(w)$ can be computed by means of the Gärtner-Ellis theorem, which states that $I(w)$ is given by the Legendre-Fenchel transform of $f(s)$ in Eq. (34), namely

[TABLE]

where the infimum has to be taken within the domain in which $f(s)$ is defined. Note that, once the large deviation principle $p(w)\sim{\rm exp}[-NI(w)]$ is satisfied, the Gärtner-Ellis theorem can be heuristically derived by a saddle-point approximation of the inverse Laplace transform of $G(s)$ GaSi12 ; Touc09 . However, it might be difficult to prove this principle a priori in specific cases, and thus this is usually done a posteriori.

In the case of free (Fermionic and Bosonic) models, the rate function is quadratic in a neighborhood of $\bar{w}$, so that small deviations from the average intensive work have a Gaussian distribution GaSi12 : this is what one would expect from a naive application of the central limit theorem. On the contrary, for large deviations from $\bar{w}$, $I(w)$ might differ significantly from its quadratic approximation, displaying interesting features. Most prominently, as shown in Refs. GaSi12 ; SoGS13 , the behavior of $I(w)$ for $w\ll\bar{w}$ becomes universal if the post-quench Hamiltonian is close to criticality, a fact which can be rationalized via a quantum-to-classical correspondence.

In the case of systems composed by free Bosonic excitations, it has been shown in Ref. GaSi12 that $p(w)$ may feature a different kind of universal behavior also for $w>\bar{w}$. In particular, as the pre-quench initial state is varied from being non-critical to critical, a transition in the form of $p(w)$ takes place such that $I(w)$ vanishes identically for $w\geq\bar{w}$ when the pre-quench initial state becomes critical. This has been identified as a “condensation” transition, in analogy to the Bose-Einstein condensation in equilibrium statistical mechanics Huang and it implies that $p(w)$ displays an algebraic decay upon increasing $w$. Although no general expression has been reported so far for such a power-law tail of $p(w)$ in this condensation regime, the latter has been shown to appear also in different non-equilibrium protocols SmSi13 . In this work we will provide a quantitative prediction for the corresponding exponent in the interaction quench introduced in Sec. II.3.

IV From the quench action approach to large deviation theory: the rate function

We now present our analysis and predictions for the statistics of the work done by the quench introduced in Sec. II.3. We begin by showing that the quench action approach allows us to demonstrate, directly and rather generally, that the large deviation principle $p(w)\sim e^{-LI(w)}$ holds and then compute the rate function $I(w)$.

We start from the expression of the moment generating function $G(s)$ in Eq. (33) in which we insert the resolution of the identity operator $\mathbb{I}$ in terms of the post-quench Bethe eigenstates $|\{\lambda_{j}\}\rangle$

[TABLE]

obtaining

[TABLE]

One then notices that Eq. (37) has a structure analogous to that of Eq. (25). Accordingly, it can be expressed as the r.h.s of Eq. (27) which involves the quench action $S_{\rm QA}[\rho]$, namely

[TABLE]

where we introduced the ground-state energy density $e_{0}(c)=E_{0}^{c}/L$, and the modified quench action

[TABLE]

with $S_{\rm QA}[\rho]$ given in Eq. (28), and $e[\rho]$ in Eq. (11). In the thermodynamic limit, the functional integral in Eq. (38) can be evaluated via the saddle-point method, leading to

[TABLE]

Here the function $\rho^{\ast}_{s}$ is determined by the saddle-point condition

[TABLE]

Note that by straightforward manipulations, Eq. (41) can be cast into the explicit form [see also Eqs. (28), (30), and (15)]

[TABLE]

involving, instead of $\rho_{s}^{\ast}$,

[TABLE]

where $\eta(\lambda)$ is defined in Eq. (9). Equation (42) has to be interpreted as follows. For each value of $s$, one finds a unique solution for the function $\varepsilon_{s}^{\ast}(\lambda)$, and hence for $\eta_{s}^{\ast}(\lambda)$. Then, by recalling that $\rho^{t}(\lambda)=\rho(\lambda)(1+\eta(\lambda))$, one plugs the latter function into Eq. (6), in order to obtain a final prediction for $\rho^{\ast}_{s}(\lambda)$. Note that the Lagrange multiplier $h(s)$ in Eq. (42) has to be chosen such that the prescribed density $D$ is obtained after using Eq. (10).

Within the saddle-point approximation in Eq. (40) one finds that, from Eq. (34),

[TABLE]

and therefore, in order to calculate $I(w)$ according to Eq. (35) one has to find the infimum, as a function of $s$, of $sw-f(s)$. When this is attained in a differentiable point $s_{w}$, it is determined by the condition

[TABLE]

Due to the concavity of $f(s)$ Touc09 , the stationary point $s=s_{w}$ can only correspond to a minimum. Using now $({\rm d}/{\rm d}s)=({\rm d\rho}/{\rm d}s)(\delta/\delta\rho)$, and exploiting Eq. (41), one can easily show that this condition is in fact equivalent to requiring

[TABLE]

As a consequence, if $f(s)$ is in addition strictly concave, the expression in Eq. (35) simplifies as

[TABLE]

Importantly, in this derivation, we never had to evaluate the quench action $S_{\rm QA}[s,\rho]$ at complex values of $s$, where it has been shown that it might display singular points PiPV17_II . Note also that the specific form of the overlap term $S_{O}[\rho]$, entering only in Eq. (42), does not play any role in this derivation. As a consequence, the latter can be generalized straightforwardly to any integrable model where the quench action approach can be applied. We note that the function $f(s)$ in Eq. (44) has been defined with a rescaling by the system size $L$ and not by the number of Bosons $N$ as in Eq. (34). The two definitions are clearly equivalent since $N=DL$ and the density $D$ is assumed to be fixed.

We can now proceed towards the explicit evaluation of the rate function $I(w)$, using Eq. (35). First, note that exploiting Eqs. (6) and (42), the action in Eq. (39) can be rewritten in the compact form

[TABLE]

Next, one needs to solve Eqs. (42) and (6) with the constraint in Eq. (10). This can be easily done numerically by standard iterative procedures. The resulting solution for $\varepsilon_{s}^{\ast}(\lambda)$ can then be plugged into Eq. (48) and integrated numerically. Finally, in order to obtain $f(s)$ in Eq. (44), one also needs to compute the ground-state energy $e_{0}(c)$. In fact, this can be written in terms of the solution of an integral equation (see, e.g., Ref. takahashi_book ). In particular, we have

[TABLE]

where $\rho_{\rm GS}(\lambda)$ satisfies the Lieb equation

[TABLE]

and where the real number $Q$ is determined self-consistently by requiring

[TABLE]

We have now all the necessary ingredients to evaluate the rate function $I(w)$, which is obtained by numerically performing the Legendre-Fenchel transform in Eq. (35). The latter expression is indeed better suited for a numerical evaluation of $I(w)$ than Eq. (47) since $f(s)$ is in general known only numerically from Eq. (44). The Legendre transform in Eq. (47) will be instead used in order to determine analytically the asymptotic behavior of $I(w)$ both close to $\bar{w}$ and for low values of $w$. We have implemented the numerical procedure outlined above, which presents no difficulty, and we have worked out analytically the asymptotic of $I(w)$; our results are summarized and discussed in the next section.

V Results

In this section we present our results for the scaled cumulant generating function in Eq. (44) and for the rate function $I(w)$. We begin by reporting in Sec. V.1 their numerical evaluation based on the exact formulas presented in the previous section, and then work out analytically their asymptotic behavior in Sec. V.2. Finally, we devote Sec. V.3 to a detailed analysis of the so-called Tonks-Girardeau (TG) limit $c\rightarrow\infty$ , which lends itself to a fully analytical treatment.

V.1 The exact rate function: numerical results

We start by presenting our numerical predictions for the scaled cumulant generating function $f(s)$, which are reported in Fig. 1($a$). The data are obtained using Eq. (44), after numerical solution of Eqs. (42), (6), and (50), which is done by standard iterative procedures Klauser2011 .

We see from Fig. 1($a$), that $f(s)$ displays many of the generic features predicted in Ref. GaSi12 . In particular, it is a concave function defined in a semi-infinite interval $[\bar{s},\infty)$ of the real line. In this case $\bar{s}=0$, since the function $f(s)$ diverges for $s<0$. This is also consistent with the fact that $f(s)$ has a singular point in its second derivative at $s=0$. To see this, one can compute the second derivative of $f(s)$, yielding

[TABLE]

as it can be verified by using the Wick theorem and by noticing that divergent terms arise on the r.h.s. of Eq. (52). On the contrary, the first derivative of $f(s)$ is finite in $s=0$ and it gives the average intensive work $\bar{w}$ performed in the quench. In particular, from Eq. (44), we have

[TABLE]

As $s$ approaches zero, $f(s)$ displays a non-analytic behavior, which will be characterized in Sec. V.2 and which is responsible for the divergence of the higher-order cumulants. We also note that the qualitative features of $f(s)$ do not change upon varying $c$. However, the average work $\bar{w}$, and hence the first derivative in $s=0$, increases upon increasing the interaction $c$. This is expected because as the repulsion among the Bosons increases, a larger energy is present in the initial BEC state where all the Bosons have zero momentum according to Eq. (18). In particular, in the Tonks-Girardeau limit $c\rightarrow\infty$ the mean intensive work $\bar{w}$ diverges.

From Fig. 1($a$) one also sees $f(s)$ slowly approaches its asymptotic value for $s\to\infty$, which increases upon increasing $c$. This asymptotic behavior, however, is difficult to analyze numerically and we postpone its discussion to Sec. V.2, where it will be determined analytically.

In Fig. 1($b$) we report our predictions for the rate function $I(w)$ corresponding to $f(s)$ in panel ($a$), which can be obtained after numerical Legendre-Fenchel transform of $f(s)$, as explained in the previous section. As we discussed above, the rate function $I(w)$ vanishes at $w=\bar{w}$, while it is identically zero for $w>\bar{w}$. This means that fluctuations for $w>\bar{w}$ must have a sub-exponential dependence $C(L,c)$ on the system size $L$. In fact, assuming a power-law decay

[TABLE]

one can constrain the value of $\beta$ by taking into account the divergence of cumulants beyond the first one, see Eqs. (52) and (53). In particular, it must be $2<\beta<3$. While the pre-factor $C(L,c)$ and the exponent $\beta$ can not be determined from large deviation theory, we will show that they can be calculated from the Bethe ansatz in Sec. VI, at least in the limit of vanishing densities of Bosons.

Finally, from Fig. 1($b$) we see that $I(w)$ rapidly approaches zero as $w\to\bar{w}^{-}$. In fact, due to the limitations in the accuracy of the numerical solutions, it is difficult to characterize this decay numerically, as $I(w)$ becomes very small when $w\simeq\bar{w}$. However, as we show in the next section, this regime can be successfully tackled analytically, so that the behavior of $I(w)$ near $\bar{w}$ can be determined exactly.

V.2 Asymptotic behavior of the rate function: analytic results

As we mentioned in Sec. III, whenever the central limit theorem applies Touc09 because $w$ can be seen as the sum of a large number of microscopic works done separately on the single particles, the rate function $I(w)$ has a quadratic expansion around $w=\bar{w}$. This is, for instance, the case for free fermionic models Silv08 ; GaSi11 ; GaSi12 ; SoGS13 ; SmSi13 , where a Gaussian distribution describes $p(w)$ for small deviations from its mean value $\bar{w}$. In the case under study, however, we show in the following that the behavior of $I(w)$ near $\bar{w}$ is not Gaussian, meaning that the central limit theorem does not apply.

In order to study the behavior of $I(w)$ for $w\rightarrow\bar{w}^{-}$, we exploit Eq. (47) which applies to our case since $f(s)$ is strictly concave and therefore $f^{\prime}(s)$ is invertible. Due to the concavity of $f(s)$, it is easy to show that the behavior of $I(w)$ near $\bar{w}$ is determined by the expansion of $f(s)$ in a neighborhood of $s=0$. In other words, we are left with the problem of determining the form of $f(s)$ for small $s$. To this end, we start from Eq. (44) and define

[TABLE]

Differentiating Eq. (44) with respect to $s$, we obtain

[TABLE]

Next, differentiating Eq. (42) with respect to $s$, multiplying each side of the resulting equation by $\rho_{s}^{\ast}(\lambda)$ and finally integrating in $\lambda$ we obtain

[TABLE]

where the Bethe equations (6) have been used. Putting everything together, we obtain the important relation

[TABLE]

Accordingly, the small-$s$ behavior of $f(s)$ is determined by that of $\int_{-\infty}^{+\infty}{\rm d\lambda}\,\rho_{s}^{\ast}(\lambda)\lambda^{2}$. Note that for $s=0$ we obtain that the derivative of $f(s)$ is the energy of the initial state, as it should. In Appendix A we show that

[TABLE]

so that we finally obtain

[TABLE]

We can now plug this expression into Eq. (47) and compute the first term in the expansion of $I(w)$ for $w\simeq\bar{w}^{-}$. By doing so, we obtain that

[TABLE]

i.e., the first term of the expansion around $\bar{w}$ is cubic instead of quadratic. As anticipated, we therefore find that small fluctuations have not a Gaussian distribution, in stark contrast with the free case.

Next, we proceed to studying the limit of $I(w)$ for small values of $w$. From Eq. (47), we see that the latter is determined by the behavior of $f(s)$ at $s\rightarrow\infty$, which we now work out analytically. This can be done by following the derivation of Refs. Takahashi1973 ; BePi2018 , where analogous calculations were done in the context of thermal equilibrium. We start by rewriting Eq. (48) as

[TABLE]

where

[TABLE]

For large $s$, the function $\varepsilon_{s}^{\ast}(\lambda)$ has two symmetric zeros which we call $Q^{\prime}$ ($-Q^{\prime}$), while we name $Q$ ($-Q$) the zeros of $\varepsilon_{\infty}(\lambda)$, defined as the solution of the limit $s\rightarrow\infty$ of Eq. (42), namely by

[TABLE]

where we made the assumption that

[TABLE]

Assuming the validity of the latter equation we write the following expansion for $h(s)$ at large $s$

[TABLE]

We computed $h^{\prime}$, $h_{0}$ and $h_{-1}$ by performing a fit against the numerical data for $h(s)$. The expansion in Eq. (66) has been numerically checked. Next, we write Eq. (63) as

[TABLE]

The first term in Eq. (67) can be studied by expanding the integrand around the points $Q^{\prime}(-Q^{\prime})$

[TABLE]

Note that the second and third term in Eq. (67) vanish as $O\left((Q-Q^{\prime})^{2}\right)$, i.e.,

[TABLE]

We now make use of the following identities, which are proven in Appendix B:

[TABLE]

where $U_{1}(\lambda)$ and $U_{2}(\lambda)$ are obtained as the solution to the following integral equations:

[TABLE]

Plugging the identities in Eqs. (70)–(73) into Eqs. (68) and (69) and then into Eq. (67), straightforward manipulations finally yield, for $s\rightarrow\infty$,

[TABLE]

with

[TABLE]

Here $v_{s}$ is the sound velocity of the system defined by

[TABLE]

The expression of $f_{0}$ in Eq. (75) coincides with the prediction of Ref. GaSi12 , since $f_{0}$ can be actually rewritten as

[TABLE]

where $|\psi_{0}^{c}\rangle$ is the ground state of the post-quench Hamiltonian $H(c)$ in Eq. (1).

The expression in Eq. (74) finally provides access to the behavior of $I(w)$ for small values of $w$. Indeed, by plugging Eq. (74) into Eq. (47) we obtain

[TABLE]

We compare this expansion with the exact rate function $I(w)$ obtained by the numerical evaluation of the formulas derived in the previous section for several values of the interaction $c$. An example is displayed in Fig. 2, where a good agreement between the two curves is manifest. Note that the leading behavior in the expansion in Eq. (78) is in agreement with the predictions of Ref. GaSi12 . Indeed, based on a quantum-to-classical correspondence, it was argued in Ref. GaSi12 that when the post-quench Hamiltonian is critical, the rate function displays the generic behavior

[TABLE]

where $d$ is the dimensionality of the system. In our case, the criticality condition is verified since the Lieb-Liniger spectrum is gapless. Accordingly, large deviations for small values of the work encode signatures of universality as predicted by the quantum-to-classical correspondence GaSi12 . A precise determination of the classical counterpart of the quantum quench analyzed here, however, goes beyond the scope of the present paper, and will not be discussed further.

V.3 The Tonks-Girardeau limit

In this section we focus on the Tonks-Girardeau limit Gira60 , corresponding to the quench where the final interactions are taken to be infinitely large. In fact, on the one hand, in this regime the formulas derived in the previous sections simplify, so that one can push the analytical control even further. On the other hand, in this limit, qualitative differences emerge in the statistics of the work, which are worth exploring per se, especially given the great relevance of this regime for cold-atomic experiments KiWW04 ; HGMD09 .

From the computational point of view, in the limit $c\rightarrow\infty$, the kernel $K(\lambda)$ in Eq. (8), entering the integral equations which characterize the quench action formalism, vanishes identically, largely simplifying the analysis. In particular, from Eq. (42), the solution for $\eta_{s}^{\ast}(\lambda)$ can be explicitly written as

[TABLE]

Here we have introduced the following parametrization of the Lagrange multiplier $h(s)$

[TABLE]

which is particularly convenient because, as shown in Ref. DWBC14 , for the quench action equations corresponding to $s=0$ in the present discussion,

[TABLE]

and therefore $h^{\ast}(0,D)=0$.

Next, from Eq. (6) one finds, in the TG limit, $\rho^{t}(\lambda)=1/(2\pi)$ and thus from $\rho^{t}(\lambda)=\rho(\lambda)(1+\eta(\lambda))$

[TABLE]

Accordingly, the density constraint in Eq. (10), which determines the parameter $h^{\ast}(s,D)$, can be written as

[TABLE]

The function $h^{\ast}(s,D)$ determined by this condition does not depend on $c$, since the latter does not appear in Eq. (84), and it actually depends on $s$ and $D$ via the combination $\tilde{s}=sD^{2}$. Equation (84) can be easily solved numerically: we report the corresponding result for $h^{\ast}(\tilde{s})$ in Fig. 3. Interestingly, the function $h^{\ast}(\tilde{s})$ appears to be almost linear in $\tilde{s}$. In fact, it is not difficult to compute the asymptotic behavior of $h^{\ast}(\tilde{s})$ for $\tilde{s}\to\infty$. Inserting Eq. (66) written up to order $O(s^{-1})$, $h_{0}^{\prime}=h_{0}-h(0)$ with $h_{0}$ and $h(0)$ defined in Eqs. (66) and Eq. (82), respectively,

[TABLE]

into Eq. (83), in the large-$s$ limit $\rho^{\ast}_{s}(\lambda)$ becomes a step function; in particular, imposing the correct density, we find

[TABLE]

and therefore the simple relation $h^{\prime}=2\pi^{2}$ in the TG limit. The behavior of $h^{\ast}(\tilde{s})$ can be analytically studied also for $\tilde{s}\rightarrow 0$, although the computations are more involved. For this reason, we present them in Appendix D, and we report here the final result, which reads

[TABLE]

The expressions in Eqs. (85) and (87) allow us to obtain directly the corresponding expansions for the scaled cumulant generating function $f(s)$. First, note that in the TG limit the latter can be expressed explicitly by plugging Eq. (80) into Eqs. (44) and (48), finding

[TABLE]

This expression can be easily evaluated numerically, as it amounts to a simple integral, once the function $h^{\ast}(sD^{2})$ is known. Notice that, rescaling $\lambda=Dy$ in the integral of Eq. (88), it turns out that $f(s)/D$ is actually a function of $\tilde{s}=sD^{2}$ only, as it has been already noticed also for $h^{\ast}(s,D^{2})$. We report the resulting data for $f(s)/D$ in Fig. 4. The same rescaling in terms of $D$ does not apply to the case at finite $c$ of Secs. V.1 and V.2 as one realizes by looking, e.g., at the small-$s$ expansion in Eq. (60). Note also that, differently from the case of finite interactions $c$, the average value of the work is infinite in the Tonks-Girardeau limit. Indeed, the expansion of $f(s)/D$ near $s=0$ differs from Eq. (60), as one realizes by plugging Eq. (87) into Eq. (88). The details of this calculation are reported in Appendix C, where we find

[TABLE]

with

[TABLE]

Analogously, the asymptotic behavior of $f(s)$ for $s\rightarrow\infty$ can be derived by plugging in Eq. (88) the expression in Eq. (85) for $h^{\ast}(\tilde{s})$. The intermediate steps are reported in Appendix D, while the final result reads

[TABLE]

with

[TABLE]

Note that the expression for $f_{1}$ in Eq. (93) is equal to the limit $c\rightarrow\infty$ of Eq. (75), as it should.

Finally, we discuss the rate function $I(w)$ in the TG limit which can be computed by numerically performing the Legendre-Fenchel transform of Eq. (88), displayed in Fig. 1(b). Given the scaling form of $f(s)/D$ as a function of $sD^{2}$, one readily obtains from Eq. (35) a scaling form for $I(w)/D$ as a function of the variable $\tilde{w}=w/D^{3}$ only. Contrary to the case of finite interactions, $I(w)/D$ never vanishes, as the average work $\bar{w}$ grows to infinity as $c\to\infty$. It is thus meaningful to study the asymptotic behavior of $I(w)/D$ for large values of $w/D^{3}$. This can be easily done by plugging into Eq. (47) the expansion in Eq. (89), finding

[TABLE]

which is plotted in Fig. 5 together with the exact numerical values of $I(w)/D$.

Before concluding this section, we note that an analogous analysis can be done for the limit $w\rightarrow 0$ of the rate function $I(w)$, by plugging Eq. (91) into Eq. (47). In this case, we find that $I(w)$ behaves as in Eq. (78), showing that fluctuations for small values of $w$ are not qualitatively affected by considering the TG regime $c\rightarrow\infty$.

VI Algebraic behavior of $p(w)$ for large $w$

In the previous sections we have quantitatively analyzed the rate function $I(w)$, characterizing the exponential decay of the distribution function $p(w)$ for $w<\bar{w}$ as $L$ grows. On the other hand, for the quench considered here, we saw that $I(w)$ vanishes identically for $w>\bar{w}$, so that, in this regime, the decay of $p(w)$ as a function of $L$ is sub-exponential. As we have anticipated in Sec. IV, assuming (for large $w$) a power-law decay $p(w)\sim C(L,c)\,w^{-\beta}$ one can bound the value of $\beta$ based on the divergence of cumulants of $w$ higher than the first. However, the leading behavior in $L$ of $C(L,c)$, characterizing $p(w)$ for $w>\bar{w}$, and the value of the exponent $\beta$, can not be obtained by large deviation theory and a more sophisticated analysis has to be carried out.

In order to determine such a leading behavior, one could start from an exact expression for $p(w)$ at finite sizes, and then perform the correct asymptotic analysis. This strategy, however, appears to be unpractical, as the exact computation of $p(w)$ at finite sizes is a significant challenge. Nevertheless, in this section we show that this problem can be solved assuming a vanishingly small density of particles. In particular, we consider the limit $L\to\infty$, with the number $N$ of particles kept fixed. We will show that, in this regime, one can extract the exact leading dependence on $L$ of $p(w)$ for large values of $w$.

From the technical point of view, the reason why the problem becomes tractable in this regime lies in the simplified structure of the solution to the Bethe equations (3), parametrizing the eigenstates of the Hamiltonian. Indeed, fixing the quantum numbers $I_{j}$ in Eq. (4), one immediately obtains the following expansion as $L$ grows

[TABLE]

with $j=1\,,\ldots N$. Namely, at the leading order in $L$, the rapidities coincide with the quasi-momenta of a free quantum gas confined within a ring of length $L$ (with the condition $I_{j}\neq I_{k}$ for $j\neq k$). The corresponding eigenvalues of the Hamiltonian become

[TABLE]

For simplicity, we consider below the case in which the number $N$ of particles is odd, so that the quantum numbers $I_{j}$ are integer, i.e., we choose

[TABLE]

A completely analogous derivation applies to the case of even $N$.

As a first ingredient for the computation of $p(w)$, we consider the zero-density limit of the overlap in Eq. (21), which has been already studied in the literature BDWC14 ; DoCa12 ; CaLe14 . In particular one finds

[TABLE]

As a consistency check, one should verify that, using Eq. (98) and keeping only the leading term of the rapidities $\lambda_{j}\simeq 2\pi I_{j}/L$, one obtains the correct values for the normalization and the energy of the initial state, i.e.,

[TABLE]

and

[TABLE]

In fact, Eq. (99) can be established analytically on the basis of Eq. (98) by using the identity (see, e.g., Ref. Zagier2012 )

[TABLE]

while we checked that also Eq. (100) is fulfilled by numerically performing the sum for small particle numbers.

We have now all the ingredients to determine the leading behavior of $p(w)$. Our strategy consists in a direct computation based on the definition in Eq. (32) which, for the Lieb-Liniger model, reads

[TABLE]

Let us fix the value $w$, so that the extensive work is $W=wL$, and let $\epsilon\ll W$ be a small energy shell (more precisely, we choose $\epsilon=\tilde{\epsilon}L$ with $\tilde{\epsilon}\ll w$). Then, the definition in Eq. (102) directly yields

[TABLE]

Note that since $w=O(L^{0})$, we are in the regime $w\gg\bar{w}$, since in the zero-density limit $\bar{w}=O(L^{-2})$. We can then proceed to evaluate the sum over the energy shell in Eq. (103) and obtain the behavior of $P(W)$ and hence of $p(w)$. In order to simplify the discussion, we start by illustrating the main idea of this derivation in the simplest case where $M=1$, i.e., $N=3$. The generalization to an arbitrary number $N$ of particles, which does not bear conceptual complications, is presented in Appendix E. At the end of this section we will report the final result of this analysis.

For $M=1$, Eq. (103) can be rewritten as

[TABLE]

where $I_{\rm min}$ and $I_{\rm max}$ are determined, via Eq. (96), by the boundaries of the energy shell in Eq. (103), i.e.,

[TABLE]

Here we dropped the ground state energy $E_{0}^{c}\sim 1/L^{2}$ since it is sub-leading with respect to $W$. Consequently we have

[TABLE]

where

[TABLE]

The series in Eq. (107) can be easily bounded as

[TABLE]

where

[TABLE]

Plugging Eqs. (105)–(108) into Eqs. (110) and (111), we immediately get

[TABLE]

Accordingly, at the leading order in $\epsilon$, we find $P_{\rm min}=P_{\rm max}=P(W)$. We recall now that the probability distribution function $p(w)$ of the intensive work is related to $P(W)$ by $p(w)=L\,P(wL)$. Accordingly, from Eq. (113), we obtain

[TABLE]

which has the form anticipated in Eq. (54) with $C(L,c)=c^{2}/L^{5/2}$, $\beta=5/2$ and clearly displays that the dependence $C(L,c)$ of $p(w)$ on $L$ is not exponential. The value of $\beta$ furthermore satisfies the bound $2<\beta<3$, anticipated in Sec. V.1. Finally, we show in Appendix E that the same reasoning can be repeated for an arbitrary number $N$ of particles and that the behavior in Eq. (114) is valid for all $N$ (with an $N$-dependent numerical prefactor). To our knowledge, this constitutes the first quantitative calculation of the power-law tail of $p(w)$, first predicted in Ref. GaSi12 , for interactive bosonic systems starting from a critical initial state. Note that, even if the result in Eq. (114) holds for an arbitrary finite number of particles $N$, its validity is still limited to the regime of vanishing density $D$: indeed, by construction, the large-$L$ limit is taken while $N$ is kept finite. In fact, in the case of finite $D$, while the bound $2<\beta<3$ continues to hold, we can not make any statement on the exact value of the exponent without further assumptions.

VII Conclusions

In this work we studied the large-deviation statistics of the intensive work $w$ done by an interaction quench of the one-dimensional Lieb-Liniger model, focusing on the case in which the initial state is the ground state of the non-interacting gas. By means of the quench action approach, we have shown that, for $w<\bar{w}$, the large-deviation principle applies to the probability $p(w)$, as it depends exponentially on the system size $L$ with $p(w)\sim{\rm exp}[-LI(w)]$, and that the Gärtner-Ellis theorem employed in Ref. GaSi12 can be used in order to determine the corresponding rate function $I(w)$. We have provided a fully quantitative analysis of the latter, working out analytically its behavior for small values of the intensive work $w$, cf. Eq. (78), and close to the average work $\bar{w}$, cf. Eq. (61). Interestingly, we have shown that for $w\simeq\bar{w}^{-}$ fluctuations are not Gaussian, in contrast to what would be expected from a direct application of the central-limit theorem. Furthermore, we analyzed the probability distribution function of the intensive work $p(w)$ for $w\gg\bar{w}$ where the large-deviation principle is violated and $p(w)$ has a sub-exponential dependence on $L$. Using an exact Bethe ansatz representation of the eigenstates of the Hamiltonian, we have derived the power-law decay of $p(w)$ in the regime of vanishing particle density, see Eq. (114), providing the first quantitative calculation of the power-law tail of $p(w)$ for interactive bosonic systems starting from a critical initial state.

Contrary to other works SoGS13 ; PaSo14 ; RyAn18 , our approach to derive $p(w)$ is not based on the Fourier transform of the so-called Loschmidt echo evaluated at real times. This is an important point, as the latter has proven to be especially hard to compute PiPV17_II ; PiPV18 , due to the presence of points of “non-analyticity” arising in its real-time dynamics HePK13 . For this reason the quench action approach presented in this work can be straightforwardly generalized to a wide class of quantum quenches in other interacting integrable models. A particular interesting example would be the prototypical $XXZ$ Heisenberg chain, where, for instance, the quantum-classical correspondence mentioned in Sec V.2 could be investigated in detail.

Finally, it would be interesting to investigate the statistics of the work done by quenches to the attractive regime of the Lieb-Liniger model, where intriguing phenomena, such as the formation of multi-particle bound states PiCE16 , have been predicted. While the study of arbitrary attractive interactions might be challenging due to the emergence of singularities in the spectrum of the Hamiltonian takahashi_book , we believe the weakly attractive regime investigated in Refs. CaCR00 ; CaCR00_II ; KaSU03 ; FlFP16 ; PiCa16 to be within the reach of the techniques presented in this work.

VIII Acknowledgments

G.P. thanks D.Rossini for fruitful discussions. L.P. acknowledges support from the Alexander von Humboldt foundation and from the Deutsche Forschungsgemeinschaft under Germany’s Excellence Strategy – EXC-2111 – 390814868. A.G. acknowledges fruitful discussions with A.Silva and F.H.L.Essler.

Appendix A Small-$s$ asymptotics of the scaled cumulant generating function

In this appendix we derive Eqs. (59) and (60), characterizing the behavior of the scaled cumulant generating function $f(s)$ as $s\rightarrow 0$. The starting point is Eq. (58) in the main text and therefore the expansion in $s$ of $\rho_{s}^{\ast}(\lambda)$ is needed. The latter can be obtained from the expansions in $s$ of $a_{s}(\lambda)$ and $\rho_{s}^{t}(\lambda)$ that we now perform.

First, from the integral equation of $\eta_{s}^{\ast}(\lambda)$ we can write that for $a_{s}(\lambda)$, which reads

[TABLE]

As a first step, we write the formal expansions of $\log a_{s}(\lambda)$ and $h(s)$ as a function of $s$:

[TABLE]

Plugging these into Eq. (115), we obtain a system of integral equations, one for each successive order in the expansion in $s$. The order [math] gives the same quench action equations solved in Ref. DWBC14 . To clarify the procedure we write in addition the result at first order

[TABLE]

Since, for each fixed value of $s$, the driving term of Eq. (115) grows as $\lambda^{2}$ when $\lambda\rightarrow\infty$, it follows that $\log a_{s}(\lambda)$ increases at most as $\lambda^{2}$ in the same limit, and therefore the following expansion as a function of $\lambda$ can be written:

[TABLE]

where $\beta_{2j}(s)=O(s)$ and $j\leq 1$ is an integer number. In particular, one has

[TABLE]

Accordingly, from Eq. (119) one has

[TABLE]

Since $a_{s}(\lambda)$ equals $a_{0}(\lambda)$ in $s=0$ we have that

[TABLE]

with $\delta_{n,0}$ the Kronecker delta symbol. Next, since the driving term of Bethe equations (6) is $1/2\pi$, with a reasoning analogous to the one done to justify Eq. (119) one has the asymptotic expansion for $\rho_{s}^{t}$ for large $\lambda$

[TABLE]

with suitable coefficients $\gamma_{2n}(s)$ whose explicit expression we do not need for the present calculation. We can eventually use the results in Eqs. (121) and (123) into the integral in the r.h.s. of Eq. (58), that we conveniently write as

[TABLE]

Since the first integral on the r.h.s. has a finite support, one can expand $\rho_{s}^{\ast}(\lambda)$ in a power series in $s$. Each term can be integrated without divergences, so that after integration the result is expected to have the form of a power series in $s$. Hence, no term proportional to $\sqrt{s}$ will arise from this contribution. This is not the case for the second integral on the r.h.s. as we see in the following. First, we write it as

[TABLE]

where we used $\rho(\lambda)=\rho^{t}(\lambda)a(\lambda)/(1+a(\lambda))$ for $\rho^{\ast}_{s}(\lambda)$ and Eq. (123). It is not difficult to see that substituting $a_{s}(\lambda)$ with $a_{0}(\lambda)$ in the denominator of the integrand on the r.h.s. of Eq. (125) leads only to corrections of order $O(s)$. As we are interested in the emergence of terms of order $O(\sqrt{s})$, we are allowed to perform this substitution. Next we can use Eq. (121) to express $a_{s}(\lambda)$ in the numerator of the resulting integrand of Eq. (125) and then note that DWBC14

[TABLE]

such that Eq. (125) becomes

[TABLE]

where we used Eq. (122) for $\alpha_{0}$, which again leads to corrections of order $O(s)$, while the coefficients $z_{2n}(s)$ can in principle be derived from $\alpha_{2n}(s)$ and $\gamma_{2n}(s)$ introduced in Eqs. (121) and (123), respectively. The integration in Eq. (127) can be performed term by term and only the first one yields a contribution $O(\sqrt{s})$ which is easily computed. Putting the latter result together with Eq. (53), which fixes the term linear in $s$, one arrives at Eq. (59), which we also tested numerically.

Appendix B Large-$s$ asymptotics of the scaled cumulant generating function

In this appendix we provide a detailed derivation of Eqs. (70)–(73) reported in the main text.

We start from the equation (42) for $\varepsilon_{s}^{\ast}(\lambda)$ which in the large-$s$ limit reads as

[TABLE]

where we used the large-$s$ asymptotic of $h(s)$ in Eq. (66). Taking the difference between Eq. (128) and Eq. (64) we have

[TABLE]

where $\varepsilon_{\infty}^{-}(\lambda)$ is defined as

[TABLE]

The last integral in Eq. (129) can be decomposed in a way similar to Eq. (67), namely,

[TABLE]

with $Q^{\prime}$ and $Q$ having the same meaning as in Sec. V.2. The analysis of the integrals appearing on the right hand side is completely analogous to the one carried out in Sec. V.2 for Eqs. (68),(69) and (71). In conclusion from Eqs. (129),(131), $\delta\varepsilon_{s}^{\ast}(\lambda)=\varepsilon_{s}^{\ast}(\lambda)-\varepsilon_{\infty}(\lambda)$ satisfies the following integral equation

[TABLE]

With a reasoning analogous to the one done in Ref. BePi2018 one can show that the term containing $(Q^{\prime}-Q)^{2}$ is at least of order $O(s^{-2})$. Accordingly,

[TABLE]

where the function $U_{1}(\lambda)$ is the solution of the integral equation

[TABLE]

In particular, computing Eq. (133) for $\lambda=Q^{\prime}$ and expanding it in the difference $Q^{\prime}-Q$ we get

[TABLE]

Note that this result shows that $(Q^{\prime}-Q)^{2}$ is exactly of order $O(s^{-2})$. As a consequence, in order to determine $\delta\varepsilon_{s}^{\ast}(\lambda)$ up to the second order in $1/s$ in Eq. (132) we keep the terms containing $Q^{\prime}-Q$. Exploiting the first-order result reported in Eq. (133), we finally obtain

[TABLE]

completing the derivation of Eqs. (70)–(73) in the main text.

Appendix C Small-$s$ asymptotics of the scaled cumulant generating function in the Tonks-Girardeau limit

In this section we study the small-$s$ asymptotic behavior of the function $h^{\ast}(s,D^{2})=h^{\ast}(\tilde{s})$ defined in Eq. (81), we name for brevity $\tilde{s}=sD^{2}$ in this appendix.

We start by taking the derivative with respect to $\tilde{s}$ of both sides of Eq. (84), obtaining

[TABLE]

which is equivalent to

[TABLE]

The right-hand side of this equation can be rewritten as

[TABLE]

We analyze the two terms $H_{1,2}(\tilde{s})$ in the limit $\tilde{s}\rightarrow 0$. For $H_{2}$, we simply compute

[TABLE]

where we have used that $h^{\ast}(\tilde{s})\rightarrow 0$ as $\tilde{s}\rightarrow 0$, see Eq. (81). For $H_{1}$, instead, we have

[TABLE]

Similarly, the integral on the left-hand side of Eq. (139) can be straightforwardly evaluated for $\tilde{s}\rightarrow 0$. Taking into account Eqs. (141) and (142), Eq. (139) becomes, for $\tilde{s}\rightarrow 0$,

[TABLE]

and therefore, after integration in $\tilde{s}=sD^{2}$

[TABLE]

which is Eq. (87). The small-$\tilde{s}$ asymptotic behavior of the scaled cumulant generating function follows from this result. In particular, plugging Eq. (87) into Eq. (88) and performing the change of variable $\lambda=yD$ inside the integral, the latter is rewritten as

[TABLE]

Taking now the derivative with respect to $\tilde{s}$ of both sides, one has

[TABLE]

where $F^{\prime}(\tilde{s})$ stands for the derivative w.r.t. $\tilde{s}$ of the integral appearing in Eq. (145). The asymptotics of this integral can be worked out using the same steps as above for the function $h^{\ast}(\tilde{s})$. In particular, defining for convenience $C=16/\sqrt{2\pi}$, we write

[TABLE]

where

[TABLE]

We analyze $I_{1}(\tilde{s})$ and $I_{2}(\tilde{s})$ separately. For $I_{1}(\tilde{s})$ in the limit $\tilde{s}\rightarrow 0$, we have

[TABLE]

while for $I_{2}(\tilde{s})$

[TABLE]

For $J_{1,2}$ introduced above we have

[TABLE]

while

[TABLE]

Collecting all the terms in Eqs. (150) and (152), Eq. (146) yields

[TABLE]

After integration, we finally obtain

[TABLE]

i.e., Eq. (89).

Appendix D Large-$s$ asymptotics of the scaled cumulant generating function in the Tonks-Girardeau limit

In this appendix we provide details of the calculation leading to the large-$s$ expansion of the scaled cumulant generating function $f(s)$ reported in Eq. (91). As well as in Appendix. C we denote for brevity $\tilde{s}=sD^{2}$.

We start from Eq. (88) with $h^{\ast}(\tilde{s})$ given by Eq. (85) and $h^{\prime}=2\pi^{2}$. One rewrites it as

[TABLE]

where

[TABLE]

where we performed the change of variable $\lambda=Dy$ inside the integral in Eq. (88). The value $\lambda=Q=D\pi$ ($y=\pi$) is the TG limit of the analogous symbol introduced in Sec. V.2 and then in Appendix B. First, we rewrite $f_{<}(\tilde{s})/D$ as

[TABLE]

with $f_{0}$ given in Eq. (92). Expanding the logarithm in a power series, one has:

[TABLE]

where the last integral, after integration by parts, can be estimated as

[TABLE]

Plugging Eq. (161) into Eq. (160), we get

[TABLE]

The evaluation of $f_{>}(\tilde{s})/D$ proceeds along the same lines: in particular, after expanding the logarithm, we can write

[TABLE]

Summing now the results in Eqs. (162), (163), and using the identity (see, e.g., Ref. NIST )

[TABLE]

we finally arrive at Eq. (91).

Appendix E Algebraic behavior of $p(w)$ for large $w$: arbitrary particle number

In this appendix we show how to extend by induction the computation presented in Sec. VI for $N=3$ to an arbitrary number $N$ of particles. In order to simplify the discussion we will first present the explicit example $M=2$ ($N=5$), and then treat the general case. As a result of the analysis of this appendix, we conclude that Eq. (114) holds for arbitrary values of $N$.

In the case $M=2$, Eq. (103) reads

[TABLE]

where the domain $\mathcal{D}$ of the double sum is determined, via Eq. (96), by the boundaries of the energy shell in Eq. (165) and the fact the two quantum numbers $I_{1}$ and $I_{2}$ have to be different. In particular

[TABLE]

corresponding to the region highlighted in blue in the $I_{1}-I_{2}$ plane shown in Fig. 6. For convenience we will neglect from the start the contribution of the domain $\mathcal{D}$ where $I_{1}>I_{1}^{\rm max}=L\sqrt{(W-\epsilon)}/(4\pi)$, depicted in red in Fig. 6, which provides a contribution $O(\epsilon^{2})$. Then, from Eqs. (165),(95) and (98) one has

[TABLE]

where $I_{2}=I_{2}^{\rm min}(I_{1})$ and $I_{2}^{\rm max}(I_{1})$ are obtained from the on-shell condition in Eq. (166), in analogy with the case $M=1$, and read

[TABLE]

The sum $\mathcal{S}(I_{1})$ over $I_{2}$ in Eq. (167) can be bounded using the same argument as the one presented in Sec. VI for the case $M=1$ to get Eqs. (110) and (111); in particular at leading order in $\epsilon$

[TABLE]

In terms of $\mathcal{S}(I_{1})$ in Eq. (169), $\epsilon P(w)$ in Eq. (167) can be written as

[TABLE]

where $\mathcal{S}(I_{1})$ can be bounded as

[TABLE]

where the symbol $\propto$ henceforth indicates that we are neglecting numerical prefactors. From Eq. (171) it follows that

[TABLE]

and, equivalently, for the probability density function $p(w)=LP(wL)$ of the intensive work $w$

[TABLE]

In order to get the leading behaviour of $p(w)$ as a function of $L$, as a final step, we need to estimate the asymptotic of the sum over $I_{1}$ in Eq. (173) for $L$ large. We have

[TABLE]

with

[TABLE]

and

[TABLE]

Plugging this into Eq. (173) we arrive at the final result

[TABLE]

At this point it should be clear how to generalize the result in Eq. (177) by induction to the general case of $N=2M+1$ particles. Indeed, consider the expression for $P_{M}(W)$

[TABLE]

Here the sum in the first line is over the $M$-dimensional shell $\mathcal{D}$ defined by

[TABLE]

while $\mathcal{D^{\prime}}$, appearing in the second line, is the same as $\mathcal{D}$ with $W\rightarrow W^{\prime}=W-8\pi^{2}I_{1}^{2}/L^{2}$ and $I_{1}^{\rm max}=L\sqrt{(W-\epsilon)}/(\sqrt{8M}\pi)$. Then, we note that Eq. (178) can be written in terms of $P_{M-1}$ as

[TABLE]

Exploiting now the induction hypothesis, we have

[TABLE]

where again we neglect numerical prefactors depending on $N$. Then using Eqs. (174), (175) and (176) it follows that $p_{M}(w)$ displays the large-$L$ behavior in Eq. (177), concluding the derivation.