Quantum Work of an Optical Lattice

TL;DR

This paper investigates the global non-equilibrium behavior of a Lieb-Liniger gas released from an optical lattice by calculating the work distribution, revealing universal properties, interaction effects, and the role of bound states.

Contribution

It provides a Bethe Ansatz-based analysis of the work distribution and Loschmidt amplitude for the Lieb-Liniger gas post-quench, including attractive interactions and bound state effects.

Findings

Universal edge exponents in work distribution

Interaction-dependent decay of Loschmidt echo

Bound states enable work extraction with low probability

Abstract

A classic example of a quantum quench concerns the release of a interacting Bose gas from an optical lattice. The local properties of quenches such as this have been extensively studied however the global properties of these non-equilibrium quantum systems have received far less attention. Here we study several aspects of global non-equilibrium behavior by calculating the amount of work done by the quench as measured through the work distribution function. Using Bethe Ansatz techniques we determine the Loschmidt amplitude and work distribution function of the Lieb-Liniger gas after it is released from an optical lattice. We find the average work and its universal edge exponents from which we determine the long time decay of the Loshcmidt echo and highlight striking differences caused by the the interactions as well as changes in the geometry of the system. We extend our calculation to the attractive regime of the model and show that the system exhibits properties similar to the super Tonks-Girardaeu gas. Finally we examine the prominent role played by bound states in the work distribution and show that, with low probability, they allow for work to be extracted from the quench.

Equations

$$\mathcal{P}(W)=\sum _n\delta \left(W-(E_n-{ϵ}_{\text{i}})\right)|⟨n|{\Psi }_i⟩{|}^2.$$

$$\displaystyle H=\int\mathrm{d}x\,\Bigg{\{}b^{\dagger}(x)\left[-\frac{\partial_{x}^{2}}{2m}\right]b(x)+c\,b^{\dagger}(x)b(x)b^{\dagger}(x)b(x)\Bigg{\}}$$

$$\int {\mathrm{d}}^{N}x\prod _{i<j}\frac{k_i-k_j-ic\text{sgn}(x_i-x_j)}{k_i-k_j-ic}\prod _l^N{e}^{i{k}_l{x}_l}{b}^{†}(x_l)|0⟩.$$

$$k_j=\frac{2\pi }{L}{n}_j-\frac{1}{L}\sum _j^N\phi (k_j-k_k).$$

$$|{\Psi }_i⟩=\int {\mathrm{d}}^{N}x\prod _{j=1}^N{\left[\frac{m\omega }{\pi }\right]}^{\frac{1}{4}}{e}^{-\frac{m\omega }{2}(x_j-{\bar{x}}_j{)}^2}{b}^{†}(x_j)|0⟩$$

$$1_N=\sum _{n_1<\cdots <n_N}\frac{|\{n\}⟩⟨\{n\}|}{\mathcal{N}(\{n\})}$$

$$\mathcal{N}(\{n\})=\det \left[{\delta }_{jk}\left(L+\sum _{l=1}^N{\phi }^{\prime }(k_j-k_l)\right)-{\phi }^{\prime }(k_j-k_k)\right].$$

$$1_N=\sum _{n_1,\dots ,n_N}\frac{|\{n\}⟩\left(\{n\}\right|}{\mathcal{N}(\{n\}}.$$

$$\left|\{n\}\right)=\int {\mathrm{d}}^{N}x\theta (\vec{x})\prod _l^N{e}^{i{k}_l{x}_l}|0⟩$$

$$|{\Psi }_i⟩={\left[\frac{4\pi }{m\omega }\right]}^{\frac{N}{4}}\sum _{n_1,\dots ,n_N}\frac{e^{-{\sum }_{j=1}^N\left[\frac{k_j^2}{2m\omega }+i{k}_j{\bar{x}}_j\right]}}{\mathcal{N}(\{n\})}|\{n\}⟩.$$

$$\mathcal{P}(W)=\sum _{n_1<\cdots <n_N}\delta \left(W-(E(\{n\})-{ϵ}_{\text{i}})\right)\frac{|⟨\{n\}|{\Psi }_i⟩{|}^2}{\mathcal{N}(\{n\})}$$

$$\mathcal{P}(W)={\int }_{-\infty }^{\infty }\frac{\mathrm{d}t}{2\pi }{e}^{iWt}\mathcal{G}(t).$$

$$\mathcal{G}(t)$$

$$G(\{n\})=\sum _{P\in S_N}(-1{)}^P{e}^{-i{\sum }_j^N{k}_j({\bar{x}}_j-{\bar{x}}_{Pj})-i{\sum }_{(j,k)\in P}\phi (k_j-k_k)}$$

$$k_j={\left[1+\frac{2\rho }{c}\right]}^{-1}\frac{2\pi }{L}{n}_j,$$

$$\mathcal{N}(\{n\})=L^N\left(1+(N-1)\frac{2\rho }{c}\right)$$

$$\mathcal{G}(t)$$

$$\mathcal{P}(W)$$

$$\mathcal{P}(W)\sim \frac{e^{-\frac{2W}{\omega }}}{W\Gamma (N/2)}{\left[\frac{2W}{\omega }\right]}^{\frac{N}{2}}.$$

$$⟨e^{-\beta W}⟩={\left(\frac{1}{1+\frac{\omega \beta }{2}}\right)}^{\frac{N}{2}}\left[1+\sum _{P\ne 1}(-1{)}^P\frac{e^{-\frac{\omega {\alpha }_P^2}{2+\omega \beta }}}{\Gamma (N/2)}\right]$$

$$\mathcal{P}(W)\sim W^{\frac{N^2}{2}-1}.$$

$$1_N={\int }_{\Gamma }\frac{{\mathrm{d}}^{N}k}{(2\pi {)}^N}|\{k\}⟩\left(\{k\}\right|$$

$${\mathcal{P}}_{c<0}(W)$$

$${\mathcal{P}}_{\text{bound}}(W)\approx N\sqrt{\frac{2\pi \omega }{m}}\frac{e^{-|c|\delta -\frac{2W}{\omega }}}{\Gamma \left(\frac{N}{2}-1\right)}{\left[\frac{2(W+\frac{|c{|}^2}{4m})}{\omega }\right]}^{\frac{N}{2}-2}.$$

$$e^{-i(N-1)\delta k_1-i{\sum }_j\phi (k_1-k_j)}$$

$$\mathcal{P}(W)\sim W^{\frac{N}{2}}.$$

Full text

I Introduction

The quantum quench is one of the simplest protocols of non-equilibrium quantum physics. An adiabatically closed system is initially prepared in some state $\left|\Psi_{i}\right>$, typically an eigenstate of a Hamiltonian $H_{i}$. At a given moment, $t=0$, the parameters of the system are suddenly changed and the system time evolves under a new Hamiltonian $H$. The sudden quench of the system parameters from the initial to the final Hamiltonian excites states throughout the spectrum and in doing so a truly non equilibrium situation is createdCalabrese and Cardy (2005); Polkovnikov et al. (2011); Mitra (2018). A classic quantum quench experiment concerns the release of gas or a Bose Einstein condensate from an optical lattice which is suddenly removed allowing the gas to expand Bloch et al. (2008); Greiner et al. (2002); Will et al. (2010). Experiments of this type address many questions that concern the dynamics of these systems, their entanglement, entropy production or thermalization, to name a few.

Low dimensional systems are of particular interest in this regard. The enhanced quantum fluctuations in such systems give access to the study of various strongly correlated phases that are hard to reach in higher dimensions. Many of these low dimensional systems are described by integrable Hamiltonians which facilitates their study by powerful analytic methods such as the Bethe Ansatz and conformal field theory Andrei (2016); Vidmar and Rigol (2016); Calabrese and Cardy (2016); Essler and Fagotti (2016); Caux (2016). A major focus has been the study of the local properties of the quenched system and in particular the behavior of local observables and correlation functions. In parallel, it was understood that a quench constitutes a thermodynamic process and within this context one can examine the concepts of work, entropy and heat of a far from equilibrium quantum system Talkner et al. (2007); Silva (2008); Goold et al. (2018); Kim et al. (2011); Chen et al. (2018). These global properties of post quench systems will be our main concern here. In particular, we will consider a gas of neutral bosonic atoms in a one dimensional trap described by the Lieb-Liniger Hamiltonian. The atoms are initially in the ground state of a 1D optical lattice and then suddenly released. We will calculate the work distribution of the quench, examining both the repulsive and attractive regimes as well as open and periodic boundary conditions.

As always, the work is given by the difference between two measurements of the energy, one pre- and the other post- quench, $W=E_{f}-\epsilon_{i}$. 111$dW=dE$, with $dQ=0$ as the system is isolated. But while the initial energy $\epsilon_{i}$ is given, the final energy may be any of the eigenvalues of the post quench Hamiltonian, $E_{n}$ which can be measured with probability, $P_{n}=|\,\langle n|\Psi_{i}\rangle\,|^{2}$. This renders quantum work a random variable with a probability distribution defined as, Talkner et al. (2007); Silva (2008)

[TABLE]

Here $\left|n\right>$ are the eigenstates of $H$ with energy $E_{n}$ and $\epsilon_{\text{i}}$ is the initial energy.

Much is known about the form of this distribution function in a number of models and quenchesSilva (2008); Gambassi and Silva (2011, 2012); Sotiriadis et al. (2013); Smacchia and Silva (2013); Chenu et al. (2018) including certain limits and approximations of the optical lattice quench discussed aboveSotiriadis et al. (2014); Pálmai and Sotiriadis (2014); Palmai (2015); Rylands and Andrei (2019). In the study of ref 26, the optical lattice is lowered but not removed entirely unlike the case we will consider here. The retention of the lattice is a crucial component as then the final Hamiltonian remains gapped and there exists a basis of long lived quasi-particles. These two features allow for a very intuitive picture of the work distribution to emerge. $\mathcal{P}(W)$ is defined for $W\geq\delta E$ with $\delta E=E_{0}-\epsilon_{\text{i}}$ being the energy difference between the ground state of $H$ and the initial state. It possesses a delta function peak at $W=\delta E$ weighted by the fidelity, $\mathcal{F}=|\left<\Psi_{i}\ \vphantom{0}\right|\left.0\vphantom{\Psi_{i}}\right>|^{2}$ signifying a transition from the initial state to the ground state. Separated from this there exists a continuum of excited states into which $\left|\Psi_{i}\right>$ can transition during the quench. The lower threshold for the continuum is at $W=2m+\delta E$ with $m$ being the mass of the lightest quasi-particle. This signifies the emission of two quasiparticles with opposite momentum from the initial state. At the threshold, $\mathcal{P}(W)$ exhibits an edge singularity similar to the Anderson and Mahan effects in the X-ray edge problemMahan (2000). When quenching from the ground state of the optical lattice in a non interacting model the distribution diverges at the threshold with exponent $\alpha=-1/2$Sotiriadis et al. (2013). On the other hand when interactions are present (even weak ones) this changes drastically to a square root singularity, $\alpha=1/2$, exemplifying the strongly correlated nature of interacting one dimensional systems Palmai (2015); Rylands and Andrei (2019). Above this, similar edge singularities will occur when new excitation channels open up e.g. the emission of 2n particles with zero momentum causes an edge singularity at $W=2nm+\delta E$ with an exponent $n/2$ when interactions are present. Below the threshold, in the region $\delta E<W<2m+\delta E$ additional delta function peaks appear if the theory supports bound states. These occur at $W=m_{b}+\delta E$, $m_{b}$ being the masses of the bound states which have the same parity as the initial state. For the ground state quench only bound states comprised of an even number of particles will appearRylands and Andrei (2019). In the thermodynamic limit the distribution becomes peaked about the average work with fluctuations vanishing as $1/\sqrt{N}$ with the most interesting features found in the region of the threshold singularity Sotiriadis et al. (2013).

Here, in contrast to study described above , we will examine the work done when the lattice is removed entirely with the post quench evolution governed by the Lieb-Liniger model, see figure 1. Since this model is gapless, the perspective of the work distribution function just elucidated is no longer correct. In a gapless theory the quench creates a macroscopic number of excitations making studies of the work statistics much more difficultVasseur et al. (2013); Perfetto et al. (2019) prompting some natural questions. What is the fate of the edge singularities when the gap is reduced to zero? Do the exponents change and and is the dependence on interactions still so dramatic? How do bound states present themselves if the model is gapless, do they just melt into the continuum? In what follows we will examine these questions as well as discuss the prominent role played by boundary conditions in this quench.

The remainder of the paper is organized as follows. In section II we introduce the model and our quench protocol. A useful identity is presented and used to calculate the time evolution of the initial state. In section III the Loschmidt amplitude is calculated for arbitrary coupling strength and particle number. From this we determine the work distribution function and examine it for strong repulsive interactions, finding the average work as well as the universal edge exponents. This analysis is extended to the attractive regime where we examine how bound states change the distribution and allow for negative values of the work to be measured. In the penultimate section we contrast the behavior in systems where periodic boundary and open boundary conditions are imposed. Finally we summarize our work, discuss generalizations of the results as well as relevance to experiment.

II An alternate identity

An excellent description of a cold atomic gas in the absence of an external potential is furnished by the Lieb-Liniger (LL) modelBloch et al. (2008). The Hamiltonian is

[TABLE]

where $b^{\dagger}(x),b(x)$ are creation and annihilation fields of bosons which have a point-like density-density interaction of strength $c$ and we set $\hbar=1$. In this article we consider both the repulsive, $c>0$, and attractive, $c<0$, regimes. The model is well known to be integrable for all couplings $c$ Lieb and Liniger (1963); Lieb (1963) and its (unnormalised) $N$-particle eigenstates can be expressed in the form

[TABLE]

The energy and momentum of such a state is $E=\sum_{j}^{N}k_{j}^{2}/(2m),~{}P=\sum_{j=1}^{N}k_{j}$ and if periodic boundary conditions are imposed the single particle momenta are quantized according to the Bethe equations,

[TABLE]

Here $\varphi(x)=2\arctan{(x/c)}$ is the two particle phase shift, $L$ is the system size and $n_{j}$ are distinct integers or half integers which serve as the quantum numbers labeling the eigenstates of the periodic system, and unless otherwise stated it should be understood that $k_{j}$ is the single particle momentum corresponding to $n_{j}$ according to (3).

Our quench protocol consists of releasing the system from an initially deep optical trap. This deep trapping potential is modelled by taking the initial state to be

[TABLE]

which is the ground state of an optical lattice of frequency $\omega$ lattice spacing $\delta$. The bosons are initially taken to be located at positions $\bar{x}_{j}$ with $\bar{x}_{j}-\bar{x}_{j+l}=\delta$. We restrict to the situation where there is at most one boson per site with the sites filled consecutively thus allowing for any value of $\rho=N/L$ provided $\rho\leq\delta$. It is further assumed that the trap is deep enough that any overlap between neighbouring sites is negligible. The initial state is then evolved according to the LL Hamiltonian.

Standard practice is to study the evolution of the initial state, $e^{-iH_{\rm LL}t}|\Psi_{i}\rangle$, by inserting a resolution of the identity in the basis of the many-body system eigenstates,

[TABLE]

with $\mathcal{N}(\{n\})$ being the norm of the Bethe states. For repulsive interactions this is given by the Gaudin formulaGaudin and Caux (2014); Korepin et al. (1993)

[TABLE]

After calculating the overlaps $C_{n}=\langle\{n\}|\Psi_{i}\rangle/\mathcal{N}$ the time evolution can be trivially performed according to $\left|\{n\},t\right>=e^{-iE(\{n\})t}\left|\{n\}\right>$. The bottleneck in this procedure occurs in the calculation of the overlaps which proves to be rather difficult. Outside of the Tonks-Girardeau (TG) limit of $c\to\infty$, where calculations are simplified Pezer and Buljan (2007); van den Berg et al. (2016); Cartarius et al. (2015); Girardeau and Wright (2000); Collura et al. (2013); Yukalov and Girardeau (2005), exact overlaps are scarce Brockmann (2014); De Nardis et al. (2014). To simplify the calculation we shall use an alternate resolution of the identity that is particularly convenient when working with initial states like (4) which are ordered in real spaceGoldstein and Andrei (2013) ,

[TABLE]

Here we have introduced the notation $\left|\{n\}\right)$ to describe an eigenstate of (II) restricted to a certain ordering in real space,

[TABLE]

where $\theta(\vec{x})$ is a Heaviside function which is non zero only for $x_{1}>x_{2}>\dots x_{N}$, and the momenta $k_{j}$ are determined by the quantum numbers $\{n\}$. It is important to note also that the ordering in the sum over quantum numbers of the system present in (5) has been removed in (6). This alternate resolution of the identity is implicitly ordered in real space as opposed to momentum space and therefore is the natural choice for calculating overlaps such as those with $\left|\Psi_{i}\right>$. Using the properties of the Bethe statesKorepin et al. (1993) it can be confirmed that this expression satisfies all the properties of a resolution of the identity.

The overlaps between (7) and (4) are now straightforwardly calculated, allowing us to express the initial state as

[TABLE]

Before proceeding further we make some comments on the above expression. The apparent ease with which we have arrived at (8) was facilitated entirely by the correct choice of identity and was further simplified by the fact that there was at most a single boson per site.

The central goal of this paper is to study the amount of work done, $W$, when the optical lattice is lowered. More precisely we will calculate the work probability distribution Silva (2008); Talkner et al. (2007). In the notation of (5) this is

[TABLE]

In the present circumstances $\epsilon_{\text{i}}=N\omega/2$ and from here on we measure the work done staring from this value, $W\to W-\epsilon_{\text{i}}$. To proceed, we introduce the Loschmidt amplitude (LA) $\mathcal{G}(t)=\left<\Psi_{i}\vphantom{\Psi_{i}}\right|e^{-iHt}\left|\Psi_{i}\vphantom{\Psi_{i}}\right>$ which is the Fourier transfrom of the work distribution

[TABLE]

The LA is a quantity of significance in its own right and central to a number of fields. The zeros of the LA define dynamical quantum phase transitions Heyl and Kehrein (2012); Heyl (2018) whilst the square of the LA, $|\mathcal{G}(t)|^{2}$, alternately known as the Loschmidt echo or return rate is prominent in studies of quantum chaotic systemsChenu et al. (2018); Gorin et al. (2006). Predominantly, we shall employ it as a calculational tool to determine $\mathcal{P}(W)$.

Using the expression (8) we find that

[TABLE]

where $G(\{n\})=\det{\left[e^{-ik_{j}(\bar{x}_{j}-\bar{x}_{k})-i\theta(j-k)\varphi(k_{j}-k_{k})}\right]}$ and $\theta(j-k)$ is a Heaviside function. Written out explicitly this is

[TABLE]

where the sum, $\sum_{P\in S_{N}}$, is over elements of the symmetric group and $(j,k)\in P$ is shorthand for pairs whose relative position is exchanged by the permutation, $j<k,~{}P(j)>P(k)$. This sum over permutations can be given the interpretation of particles exchanging positions after expanding from their original lattice positions with every exchange of particles being accompanied by the two particle phase shift, $\varphi$.

The formula (11) gives the exact Loschmidt amplitude for arbitrary $c,N$ and $L$ however its generality makes it somewhat cumbersome. One can simplify it by expanding in $c\gg m\omega$. Under this assumption $\varphi(x)\approx 2x/c+\mathcal{O}(1/c^{3})$ and we find that

[TABLE]

and in (12) we have used the fact that only states with zero momentum are present in the sum (11). For a finite size system the Loschmidt amplitude displays recurrences with a period $\tau=\left(1+2\rho/c\right)^{2}L^{2}/\pi\omega$. These recurrences disappear in the infinite volume limit.

III Work in infinite volume

We turn now to the evaluation of the work done by Fourier transforming the Loschmidt amplitude of the open system. This is the case where the occupied part of the lattice is much smaller than its overall size $\rho\ll\delta$. We shall also consider in the next section the case of a fully filled lattice, $\rho=\delta$ and periodic boundary conditions.

III.1 Repulsive interactions

In the thermodynamic limit, $N,L\to\infty$ the sum over quantum numbers becomes a product of integrals which can be evaluated giving,

[TABLE]

Where we have introduced $\alpha_{P}^{2}=m\delta_{\text{eff}}^{2}\|P\|^{2}/2$ with $\delta_{\text{eff}}=\left[1+\frac{2}{c\delta}\right]\delta$, an effective distance between lattices sites and $\|P\|^{2}=\sum_{j}^{N}(j-P(j))^{2}$. Again we can interpret the sum over permutations as the a sum over particles exchanging positions, with $\|P\|^{2}=2$ for example corresponding to a neighbouring pair exchanging positions while $\|P\|^{2}=8$ could be 4 nearest neighbour exchanges or 1 next nearest neighbour exchange.

The large repulsive interaction will inhibit the spreading of the particles which is felt through an increase in the effective distance between sites. Including higher order terms in this large $c$ expansion will cause further dressing of this distance. A similar dressing of the distance occurs in the scattering of soliton-like objects in integrable modelsVlijm et al. (2015); Doyon et al. (2018); Doyon and Spohn (2017).

Performing the Fourier transform of $\mathcal{G}(t)$ we find that the work distribution function is

[TABLE]

where $J_{n}(x)$ is a Bessel function of the first kind. Here we see that the sum over the $P$ is analogous to the sum over the number of excited particles for the gapped case which was discussed in the introduction. A notable distinction from the gapped case is the absence of a delta function peak as well as any threshold singularity at finite $W$.

We plot (15) for different values of $N$ in Fig. 2 and also the non-interacting result for the same values in Fig. 3. Comparing the two figures we see some common features as well as some striking distinctions. We note that the average of the distributions $\left<W\right>=\int\mathrm{d}W\,W\mathcal{P}(W)$ appears to be independent of the presence of interactions. Since the quench is extensive in nature we can expect that $\left<W\right>\sim N$ which is seen in the figures through the rightward shift of the distributions. Using this along with the properties of the Bessel function we find that the dominant contribution of the distribution in this region comes from the identity permutation

[TABLE]

This is the moment generating function of the Gamma distribution whose average is $\left<W\right>=N\omega/4$. As anticipated it is independent of the interaction strength. The lack of dependence on $c$ can be understood from the latter formula along with the fact the initially the bosons have negligible overlap. This in agreement with the result calculated using $\left<W\right>=\left<\Psi_{i}\vphantom{\Psi_{i}}\right|H\left|\Psi_{i}\vphantom{\Psi_{i}}\right>$. Along similar lines one can show that the $m^{\text{th}}$ moments $\mathcal{P}(W)$, with $m\ll N$ are also independent of the interaction for example the variance and skewness are $N\omega^{2}/8$ and $\sqrt{2/N}$ respectively. Going beyond this we can calculate the average exponentiated workJarzynski (1997); Talkner et al. (2007)

[TABLE]

from which we can derive all the moments of the work distribution through repeated differentiation of $\beta$.

The large $W\gg\left<W\right>$ regime appears to be also unaffected by the interactions which can be interpreted by translating back to the language of the the LA using (10). At short times the bosons expand from their initial trap positions without encountering one another and so are independent of the interactions.

At small values of $W$ the distribution is strongly affected by the interactions. We can see large resonant peaks which diminish as the particle number is increased. Moreover the distribution decays much quicker as $W\to 0$ when interactions are present. By expanding the Bessel functions and using the identity $\sum_{P}(-1)^{P}(\|P\|^{2})^{n}=0,~{}\forall n<\binom{N}{2}$Speer (2019), we have that near the threshold when the system is interacting

[TABLE]

This exponent can be interpreted as coming from the coalescence of all the edge singularities in the gapped case. As with the gapped case this dramatically different when there are no interactions present. Repeating the same calculation in the non interacting case gives instead that $\mathcal{P}(W)\sim W^{\frac{N}{2}-1}$. The value of the exponent is not dependent on the value of the interaction strength other than if it is non zero. This highlights how even weakly interacting theories exhibit strong correlations in one dimension.

The region $W\sim 0$ of the distribution gives insight to the long time behavior of the Loschmidt echo. In the presence of interactions we have that as $t\to\infty$, $|\mathcal{G}(t)|^{2}\to 1/t^{N^{2}}$ whilst in the non interacting case we have $|\mathcal{G}(t)|^{2}\to 1/t^{N}$. The much faster decay in the interacting case results from the process of dynamical fermionisationIyer et al. (2013); Iyer and Andrei (2012) wherein the interacting bosons acquire fermionic correlations as the system evolves. Consequently the bosons spread out through the trap and result in a rapidly vanishing overlap with $\left|\Psi_{i}\right>$.

III.2 Attractive interactions

We now turn to the attractive regime which is of significant interest. The properties of the attractive model both in and out of equilibrium are much less studied than its repulsive counterpart. This dearth of theoretical results stems from the increased complexity of the Bethe Ansatz solution in the attractive model. When $c<0$ the model supports bound states and the ground state consists of a single bound state of all $N$ particles McGuire (1964). While the eigenstates given by (2) and Bethe equations (3) remain valid, complex values of $k$ which correspond to bound states are allowed. The resolutions of the identity appearing in (5) and (6) also remain formally valid provided these complex valued solutions are accounted for Piroli et al. (2016); Zill et al. (2018). A large stumbling block however is that the normalisation of the Bethe states in the attractive regime is not known in closed form.

In the low density limit however it has been shown that for both repulsive and attractive interactions the spatially ordered identity (6) becomes Rupasov and Iudson (1984); Iyer and Andrei (2012); Iyer et al. (2013); Andrei (2016)

[TABLE]

where we label the eigenstates by $\{k\}$ rather than $\{n\}$. The contours of integration, $\Gamma$ lie on the real line for repulsive interactions and are spread out in the imaginary direction for the attractive case with Im$(k_{j+1}-k_{j})>|c|$.

Making use of this here in conjunction with the same $|c|\gg m\omega$ expansion we find that the work done in the attractive regime separates into two contributions,

[TABLE]

The first term $\mathcal{P}_{\text{free}}(W)$ is the contribution from particles which do not form bound states, it is identical to the expression in repulsive case given in (15) only now $c<0$. The major difference imposed by this is that the effective distance between the particles is smaller $\delta_{\text{eff}}<\delta$, the attractive interactions promoting the clustering of particles.

The simple analytic continuation to negative coupling of the first term is reminiscent of the the super Tonks-Girardeau gasAstrakharchik et al. (2005); Batchelor et al. (2005); Haller et al. (2009). This highly correlated state of the LL model is created by preparing a repulsive LL gas in the Tonks-Girardeau limit, $c\to\infty$ Tonks (1936); Girardeau (1960) and then abruptly changing the interaction strength from the being large and positive to large and negative. The result is a metastable nonequilibrium state which exhibits enhanced correlations. Many of the properties of this state can be shown to emerge from a simple analytic continuation of the coupling to large negative values Kormos et al. (2011); Astrakharchik et al. (2004); Muth and Fleischhauer (2010). In effect the negligible overlap of each particle of our initial state mimics the density profile of the TG gas and so super-TG like behaviour is not unexpected. We should stress that the expression (19) is valid at arbitrary negative values $c$ and so not limited to super-TG regime.

The second term $\mathcal{P}_{\text{bound}}(W)$ is entirely different. It is due to the bound states and is calculated by deforming the contours in (19) to the real line and picking up contributions due to the poles at $k_{i}-k_{j}=ic$ present in in (2). An $n$-particle bound state can be shown to contribute $\mathcal{P}_{n-\text{bound}}(W)\propto|c|^{n-1}e^{-n|c|\delta}$ with factors from multiple bound states being multiplicative.

This exponential factor means that the probability that the initial state transitions to one containing bound states is highly suppressed and in the true super-TG limit vanish entirely. Despite this, for finite $|c|$ the bound states have a strong signature in work distribution function. Since forming a bound state will lower the energy of the systemMcGuire (1964) the work distribution becomes non vanishing at negative values of $W$. There is a non zero probability that work can be extracted from the system. Importantly this does not violate the 2nd law of thermodynamics as the average work remains positive $\left<W\right>$ Jarzynski (1997, 2011). In fact, it has been shown observed recently that the probability of extracting work from a single electron transistor can be as high as 65% whilst still satisfying the 2nd law Maillet et al. (2018).

To see this we examine the leading term of $\mathcal{P}_{\text{bound}}(W)$ which arises due to the formation of a single two particle bound state

[TABLE]

Which is non vanishing for $-|c|^{2}/4m<W$. Determining the full bound state contribution is a straightforward yet involved calculation which we we will not deal with here.

IV Finite Density and the role of boundary conditions

It is interesting to study also the case where $\rho=\delta$ so that the initial state consists of a fully occupied lattice. In this scenario the boundary conditions play an important role and will change the behavior of the system in the region $W\ll\left<W\right>$. At finite density the Bethe equations can be used to reduce many of the terms in the sum over permutations in $G(\{n\})$ (11). For example the permutation correspond to the ordering $P=(23\dots N1)$ gives

[TABLE]

where we have used $L=N\delta$ and (3) to get the right hand side. Taking the same limits as before one arrives at (14) for the Loschmidt amplitude however the boundary conditions should be taken into account when calculating $\alpha_{P}^{2}$. For instance when $\rho\ll\delta$ the permutation $P=(23\dots N1)$ gives $\alpha_{P}^{2}=m\delta^{2}_{\text{eff}}N(N-1)/2$ however for periodic boundary conditions as a result of (22) we get instead $\alpha_{P}^{2}=m\delta^{2}N/2$. The terms corresponding to no or few particles exchanging positions are unaffected by the boundary conditions however they serve to reduce those which involve widely separated particles exchanging positions.

The work distribution function is again given by (15) in this case however as with the Loschmidt amplitude one must account for the boundary conditions when calculating $\alpha_{P}^{2}$. In the region of $W\sim\left<W\right>$ this has little effect as this area is dominated by terms corresponding to no or few exchanges of particles and the average work is as before. At small values of work however all permutations contribute and there is a difference. By expanding the Bessel functions to first order we find that instead

[TABLE]

Consequently the long time decay of the echo is given by $|\mathcal{G}(t)|^{2}\to 1/t^{N+2}$. In this instance the strongly interacting particles have no large trap to expand into as was the case before resulting in slower decay of the echo.

V Conclusion

In this paper we have studied the work done on an interacting cold atomic gas during a quantum quench. Specifically we studied the work done on the system when it is fully released form an optical lattice. Using Bethe Ansatz methods we derived an exact expression for the Loschmidt amplitude valid at all values of the interaction strength and using this studied the work distribution function at large $|c|\gg m\omega$. The edge exponents and average work were calculated and we showed that when the gas is released in a much larger trap the edge exponents differ dramatically when the system is interacting. As a consequence the Loschmidt echo decays much more rapidly in the interacting system. This calculation was then extended to the attractive regime where it was seen that the formation of bound state in the post quench system result in a small, non vanishing probability of extracting work form the optical lattice quench. Furthermore the distribution displayed properties which are indicative of the super Tonks-Girardeau gas, meaning that the non bound state contribution is obtained from analytically continuing the repulsive result to negative coupling. Finally we examined the case where the initial lattice has unit filling and periodic boundary conditions. Here it was seen that although the edge exponent did change when interactions were present the effect was not as dramatic as the case of open boundary conditions.

The calculations presented here were carried out for the simplest case of at most one particle per lattice site in the initial state. The results, especially in the attractive regime suggest that starting from higher filling number or using a model with more bound state channels such as the Gaudin-Yang gas would produce interesting resultsYang (1967); Gaudin and Caux (2014); Guan and Andrei (2018).

The experimental measurement of the work statistics of a closed quantum system after a quench has been proposed and carried out in a number of different settings Türeci et al. (2011); Latta et al. (2011); Dorner et al. (2013); Maillet et al. (2018) including in cold atom gasesCerisola et al. (2017). Although this has so far not been carried out in the setting we have proposed our results provide firm predictions for what would be observed.

Acknowledgements.

This research was supported by DOE-BES (DESC0001911) (CR) an by NSF Grant DMR 1410583 (NA).