Quenches from gaussian bosons to Tonks-Girardeau gas: stationary states and effects of a confining potential

TL;DR

This paper investigates the non-equilibrium dynamics of a Bose gas quenched from a Gaussian initial state to the Tonks-Girardeau regime, analyzing stationary states and effects of confinement using correlation functions and phase-space mapping.

Contribution

It provides a comprehensive description of steady-state correlation functions after a quantum quench to the Tonks-Girardeau limit, including effects of a confining potential.

Findings

Correlation functions expressed as Fredholm determinants

Steady-state properties derived for quenched Bose gases

Mapping of dynamics to classical phase-space in confinement

Abstract

We consider the general problem of quenching an interacting Bose gas from the noninteracting regime to the strong repulsive case described by the Tonks-Girardeau limit, with the initial state being a gaussian ensemble for the bosons. A generic multi-point correlation function in the steady state can be fully described in terms of a Fredholm-like determinant suitable both for a numerical study and for an analytical study in certain limiting cases. Finally, we extend the study to the presence of a smooth confining potential showing that, in the thermodynamic limit, the time evolution of the two-point function can be mapped to a classical problem in a properly defined phase-space.

Equations

$${\rho }_{\text{GGE}}=\frac{1}{\mathcal{Z}}{e}^{-{\sum }_j{\lambda }_j{\mathcal{I}}_j},$$

$$H_{\text{LL}}=\int dx\left[\frac{{\partial }_x{\psi }^{†}(x){\partial }_x\psi (x)}{2m}+c{\psi }^{†}(x){\psi }^{†}(x)\psi (x)\psi (x)\right].$$

$$H_{\text{LL}}=\int dx\frac{{\partial }_x{\Psi }^{†}(x){\partial }_x\Psi (x)}{2m},c=+\infty .$$

$$\Psi (x)=P_x\psi (x)P_x.$$

$$\{\Psi (x),\Psi (x)\}=0,\{\Psi (x),{\Psi }^{†}(x)\}=1$$

$$⟨{\Phi }^{†}(y_1)\dots {\Phi }^{†}(y_l)\Phi (x_l)\dots \Phi (x_1)⟩=$$

$$⟨:{\psi }^{†}(y_1)\dots {\psi }^{†}(y_l)e^{-2{\int }_{\mathcal{D}}dzn(z)}\psi (x_l)\dots \psi (x_1):⟩,$$

$${\int }_{\mathcal{D}}dzn(z)={\int }_{-\infty }^{\infty }dzn(z)(1-e^{i\pi {\sum }_{j=1}^l{\chi }_j(z)})/2,$$

$$⟨{\Phi }^{†}(y)\Phi (x)⟩=⟨:{\psi }^{†}(y)e^{-2{\int }_{\mathcal{D}}dzn(z)}\psi (x):⟩=$$

$$⟨:{\psi }^{†}(y)e^{-2{\int }_{\mathcal{D}}dzn(z)}\psi (x):{⟩}_{\text{c}}⟨:e^{-2{\int }_{\mathcal{D}}dzn(z)}:⟩,$$

$$⟨:e^{-2{\int }_{\mathcal{D}}dzn(z)}:⟩=e^{{\sum }_{j=1}^{\infty }\frac{(-2{)}^j}{j!}{\int }_{\mathcal{D}}{d}^{j}z⟨:n(z_1)...n(z_j):{⟩}_c}.$$

$$⟨:{\psi }^{†}(y)e^{-2{\int }_{\mathcal{D}}dzn(z)}\psi (x):{⟩}_{\text{c}}=$$

$$\sum _{j=0}^{\infty }(-2{)}^j{\int }_{\mathcal{D}}{d}^{j}z\prod _{l=0}^j\mathcal{C}(z_{l+1},z_l)$$

$$⟨:e^{-2{\int }_{\mathcal{D}}dzn(z)}:⟩=e^{{\sum }_{j=1}^{\infty }\frac{(-2{)}^j}{j}{\int }_{\mathcal{D}}{d}^{j}z{\prod }_{l=1}^j\mathcal{C}(z_{l+1},z_l)},$$

$$⟨:{\psi }^{†}(y)e^{-2{\int }_{\mathcal{D}}dzn(z)}\psi (x):{⟩}_{\text{c}}=\left(\frac{{\mathcal{C}}_{\mathcal{D}}}{1+2{\mathcal{C}}_{\mathcal{D}}}\right)(y,x),$$

$$⟨:e^{-2{\int }_{\mathcal{D}}dzn(z)}:{⟩}_{\text{c}}=\underset{\mathcal{D}}{\det }(1+2{\mathcal{C}}_{\mathcal{D}}{)}^{-1}.$$

$$⟨{\Phi }^{†}(y_1)\dots {\Phi }^{†}(y_l)\Phi (x_l)\dots \Phi (x_1)⟩=$$

$$\frac{1}{{\det }_{\mathcal{D}}(1+2{\mathcal{C}}_{\mathcal{D}})}\sum _{\sigma }\prod _{j=1}^l\left(\frac{{\mathcal{C}}_{\mathcal{D}}}{1+2{\mathcal{C}}_{\mathcal{D}}}\right)(y_j,x_{\sigma (j)}),$$

$$⟨{\psi }^{†}(y_1)\dots {\psi }^{†}(y_l)\psi (x_l)\dots \psi (x_1)⟩=$$

$$\underset{\mathcal{D}}{\det }(1-2{\mathcal{F}}_{\mathcal{D}})\sum _{\sigma }\text{sign}(\sigma )\prod _{j=1}^l\left(\frac{{\mathcal{F}}_{\mathcal{D}}}{1-2{\mathcal{F}}_{\mathcal{D}}}\right)(y_j,x_{\sigma (j)}),$$

$$\mathcal{C}(x,y)=\int \frac{dk}{2\pi }\frac{e^{ik(x-y)}}{e^{\beta \left(\frac{k^2}{2m}-\mu \right)}-1},$$

$$⟨{\Phi }^{†}(y)\Phi (x)⟩=\mathcal{C}(y,x)e^{-2|x-y|\left(n+\mathcal{O}(n^2)\right)}+\mathcal{O}(n^2).$$

$$⟨{\Phi }^{†}(y)\Phi (x)⟩=\frac{n}{(1+2n|x-y|{)}^2}.$$

$$\underset{\mathcal{D}}{\det }(1+2{\mathcal{C}}_{\mathcal{D}}{)}^{-1}=e^{-{\int }_x^{y}dz\left[\log (1+2{\mathcal{C}}_{\mathcal{D}})\right](z,z)},$$

$$\underset{\mathcal{D}}{\det }(1+2{\mathcal{C}}_{\mathcal{D}}{)}^{-1}=e^{-|x-y|\int \frac{dk}{2\pi }\log \left[1+2\tilde{\mathcal{C}}(k)\right]+\mathcal{O}\left(|x-y{|}^0\right)},$$

$$⟨{\Phi }^{†}(y)\Phi (x)⟩$$

$$⟨{\Phi }^{†}(y)\Phi (x)⟩\propto e^{-\xi |x-y|}\left(\frac{e^{-\lambda |x-y|}}{\lambda }+\text{c.c.}\right),$$

$$\xi$$

$$\lambda$$

$$⟨n(x,t)n(y,t){⟩}_{\text{steady state}}=n^2-|⟨{\Phi }^{†}(x)\Phi (y)⟩{|}^2$$

$$H_{LL}\to H_{LL}+\int dxV(x){\psi }^{†}(x)\psi (x),$$

$$C(s,q)=\int \frac{dk}{2\pi }\frac{e^{iks}}{e^{\beta \left(\frac{k^2}{2m}+v(q)-\mu \right)}-1}.$$

$$i{\partial }_t{F}_t=\left\{\frac{{\partial }_s{\partial }_q}{Lm}-\left[v\left(q+\frac{s}{2L}\right)-v\left(q-\frac{s}{2L}\right)\right]\right\}{F}_t.$$

Full text

I Introduction

The recent advances in experimental techniques concerning the realm of cold atomic gases exp made it urgent to reach a better understanding of the behavior of closed quantum systems when driven out of equilibrium. A paradigmatic out-of-equilibrium protocol consists in the so called quantum quenchcalabrese-cardy , i.e. the system starts in the ground state, an excited state or a thermal state of its Hamiltonian and then some parameter of the Hamiltonian is abruptly changed, resulting in non-trivial dynamics. Despite the fact that the whole system is closed and so the time evolution is non dissipative, finite subsystems are expected to relax towards a steady state and immense efforts have been devoted to investigating this hypothesispssv .

Particularly appealing is the world of one dimensional systems, due to the existence of integrable modelsKorepin , i.e. those that possess infinitely many local integrals of motion in addition to the Hamiltonian which allow an exact solution of the dynamicsintegrability_review . The presence of infinitely many conserved quantities makes the steady state retain additional information of the initial state compared to the non-integrable case so that the Gibbs Ensemble needs to be modified to a Generalized Gibbs Ensemble (GGE) density matrixgge

[TABLE]

where $\mathcal{I}_{j}$ are the relevant conserved charges and $\lambda_{j}$ the associated Lagrange multipliers. The correctness of the GGE prediction has been verified in many different casesgge_work and recent works focused in understanding the locality features of the charges to be included in the GGE constructionlocalcharges .

Despite several achievements in the study of integrable models out-of-equilibrium, apart from free models analytic results have been derived only for restricted classes of observables and initial states and at the price of tremendous effortsquench_initial_states . Recent investigations probed inhomogeneous quench protocols, obtaining analytical results in the framework of transport problemstransport .

In between the complexity of quenching a truly interacting integrable model and the simplicity of a free model lies a special interaction quench in the well known Lieb Liniger modelKorepin ; quench_liebliniger . Its Hamiltonian describes a one-dimensional gas of non-relativistic bosons with mass $m$ interacting through a density-density contact potential of strength $c$

[TABLE]

Above, the scalar fields $\psi^{{\dagger}}(x)$ satisfy the canonical bosonic commutation relations $[\psi(x),\psi^{{\dagger}}(y)]=\delta(x-y)$. For arbitrary values of $c$ the model is a truly interacting integrable model, however in two special limits the dynamics is that of free particles: at $c=0$ the model reduces to free bosons, while in the $c\to\infty$ limit known as Tonks-Girardeau gas (TG) the bosons behave as free but impenetrable particles tonksgas that are equivalent to free fermions by means of a Jordan Wigner transformation. Compared with many other quenches among free theories, the modes at $c=0$ and $c=\infty$ are not connected through a simple linear relation.

This quench protocol has already been successfully addressed in the literature, choosing as initial state the BEC state and analytical expressions for the steady state correlators were obtainedkcc14 . Subsequently, these results were generalized to include a hard wall trapmckc14 . The presence of simple analytical results for an experimentally relevant quench protocol is certainly appealing, however the above mentioned works concern a rather peculiar initial state. Among the flaws of the BEC as pre-quench state are: (i) the total absence of any space dependence in the correlation functions (except when a confining potential is present); (ii) the strong dependence on the trap shapemckc14 which makes the proper thermodynamic limit difficult to define. In this sense, the BEC state is rather pathological and therefore, it is worth to investigate a different class of states which hold a more realistic behavior. On the other hand, the initial state is still required to be simple enough to allow an analytical treatment. A good compromise that meets both requirements is given by the gaussian (in terms of the bosons) ensembles, which are the initial states we consider in this paper. In particular, we present a detailed analysis of the multipoint correlation functions in the steady state attained at late times, with emphasis on the large distance decaying properties. We initially consider an homogeneous system and subsequently a thermodynamically large trap, providing closed expressions for the correlators suitable both for a numerical and analytical study. The correlation functions are notoriously difficult to be obtained in quenches in truly interacting integrable models: the simplification of quenching towards the Tonks-Girardeau limit makes these quantities more easily accessible. While our results apply to arbitrary gaussian (in terms of the bosons) initial states, we consider the prototypical example of the free thermal state, i.e. an ensemble with Bose-Einstein distribution: the density-density correlation function in the steady state is proven to decay exponentially, the decay length being a non trivial function of the temperature and chemical potential. When the zero temperature limit is taken, the decay behavior changes from exponential to algebraic: actually, in the zero temperature limit the exact correlation function is computable. This result is striking different from what is obtained quenching from a BEC kcc14 , where an exponential decay is observed, but this is attributed to the fact that the finite temperature Gibbs Ensemble is not continuously connected to the BEC state as long as one dimension is concerned. Our discussion is organized as follows. In Section II we present the general technology needed to handle this quench protocol: under the condition of initial clustering of fermionic correlators, dephasing among the free fermionic modes of the post quench dynamics leads to a gaussian steady stategaussification in terms of the fermions (i.e. the charges appearing in the GGE are quadratic in terms of the fermionic operators). Because of the free dynamics connecting the initial state with the steady state, the open technical problems reduce to (i) computing the initial fermionic correlators on a bosonic gaussian ensemble and (ii) computing back the bosonic correlators on the fermionic gaussian steady state. We derive closed expressions in terms of Fredholm determinants to fulfill this task. These results are employed in Section III to study both numerically and analytically quenches in homogeneous systems, i.e. the particles are confined in a thermodynamically large interval with periodic boundary conditions. While being an important benchmark, homogeneous systems with periodic boundary conditions are a crude approximation of realistic experimental systems, in which an inhomogeneous confining potential is always present: Section IV extends our study to the presence of a trap. We show how in the thermodynamic limit (which is necessary for the system to exhibit relaxation) a semiclassical time evolution naturally emerges and allows us to consider traps of arbitrary shapes. Our final conclusions are gathered in Section V, while Appendix A-B contain some technical details in support to the main text.

II Quench set up

In the infinite repulsive regime, the Lieb Liniger gas (2) can be described as a free model of hard-core bosons hardcore $\Psi(x)$

[TABLE]

The infinite repulsive term prevents the possibility of having two particles in the same position and this can be achieved through a definition of the hard-core bosonic fields $\Psi(x)$, obtained projecting the standard bosonic fields $\psi(x)$ on the subspace with at most one particle in position $x$

[TABLE]

Above, $P_{x}=|0_{x}\rangle\langle 0_{x}|+|1_{x}\rangle\langle 1_{x}|$ is the projector on the subspace with no more than a particle in position $x$. The constraint is evident from the commutation rules satisfied by the hard-core bosonic fields, i.e.

[TABLE]

at the same position and the standard commutation rules $[\Psi(x),\Psi(y)]=[\Psi(x),\Psi^{\dagger}(y)]=0$ whenever $x\neq y$ kcc14 . Actually, a correct definition of the hard-core bosons requires a proper lattice regularization (see Appendix A for further details): in order to consider directly the continuum limit it is convenient to introduce the fermionic field $\Phi(x)$ through a Jordan-Wigner transformation $\Phi(x)=\exp\left(i\pi\int^{x}dz\,\Psi^{{\dagger}}(z)\Psi(z)\right)\Psi(x)$. While the hard-core bosons need a lattice regularization, on the fermionic fields the continuum limit can be taken and standard commutation rules $\{\Phi(x),\Phi^{{\dagger}}(y)\}=\delta(x-y)$ are recovered. The strong-interacting Hamiltonian becoming thus quadratic in terms of those fermionic fields.

Both models, namely the non-interacting and the strong-interacting, are therefore related via a non-local transformation which nevertheless makes a well-established connection. Whenever we need to pass from one theory to the other, e.g. via a quench protocol, the problem reduces to compute the multipoint correlators of the post-quench fields in terms of the pre-quench ones. This dictionary is provided by the underling regularized hard-core bosons, whose normal ordered correlation functions can be computed as if the hard core bosons were canonical bosonic operators: the details are reported in Appendix A, thereafter we limit ourselves to quote the necessary results. For example, the fermionic multipoint correlation function can be expressed in terms of bosonic fields as

[TABLE]

where, without loss of generality, we assumed $y_{1}<\dots<y_{l}$ and $x_{1}<\dots<x_{l}$, $n(z)\equiv\psi^{{\dagger}}(z)\psi(z)$ being the density operator, and the integration domain $\mathcal{D}$ is defined via

[TABLE]

where $\chi_{j}(z)$ is the characteristic function of the interval $(y_{j},x_{j})$. The same expression can be used to compute the bosonic correlators starting from the fermionic ones, provided in Eq. (6) the bosons are replaced with fermions and vice versa. A rigorous derivation of this relation requires a lattice regularization similarly to what has been done in Ref. kcc14, and we leave these considerations to Appendix A.

Considering a quench from homogeneous free bosons towards the Tonks-Girardeau limit, the fermions evolve as free fields, therefore the system typically exhibits gaussificationgaussification , i.e. all the fermionic multipoint connected correlators vanish apart from the two point correlator. Assuming translation invariance, the two point correlator is conserved by the time evolution and it is the only information of the initial state that survives in the steady state.

Therefore, the quench problem is reduced to computing the expectation value of the two point correlator on the pre-quench state through Eq. (6): handling this expression for arbitrary pre-quench states appears as an hopeless task, since even in the computation of the fermionic two point correlator the whole set of multipoint bosonic correlators is required, as it is clear from the power expansion of the exponential in Eq. (6). Nevertheless, some special classes of states allow for a brute force computation of Eq. (6), i.e. the already studied BEC case kcc14 and the gaussian ensembles that we are going to analyze.

Notice that, being the steady state gaussian in terms of the fermionic fields, the bosonic correlators can be obtained solving the inverse problem with respect to that we analyzed so far. Thus, we discuss the computation of fermionic correlators from a gaussian bosonic ensemble: the calculations in the other direction are completely analogous. Evaluating Eq. (6) on a gaussian ensemble amounts to an extensive use of the Wick theorem, however the calculations are greatly simplified replacing Eq. (6) with an expression in terms of the connected correlators

[TABLE]

where $\langle...\rangle_{\text{c}}$ is the connected expectation value. It must be stressed that the connected part is constructed considering the particle density $n(z)$ as a unique operator and not as a composite one $n(z)=\psi^{\dagger}(z)\psi(z)$. The derivation of (8) is a simple exercise best understood through an analogy with standard Feynman diagrams where $e^{-2\int_{\mathcal{D}}dz\;n(z)}$ plays the role of the “action”. In the computation of correlation functions the expectation value of the action factorizes and we are left with the sum of connected diagrams, i.e. (8). Besides, $\langle:e^{-2\int_{\mathcal{D}}dz\;n(z)}:\rangle$ can written in terms of connected correlators thanks to the cumulant expansion weinberg

[TABLE]

The computation of the connected correlators on a gaussian ensemble is straightforward (see Fig. 1) and we reach the following expressions (where $\mathcal{C}(y,x)=\langle\psi^{\dagger}(y)\psi(x)\rangle$)

[TABLE]

and

[TABLE]

where in the first expression $z_{j+1}=y$ and $z_{0}=x$, while in the second expression $z_{j+1}=z_{1}$. Eq. (10-11) are readily resummed in compact expressions once we define $\mathcal{C}_{\mathcal{D}}$ as the linear operator obtained restricting the correlator $\mathcal{C}$ on the domain $\mathcal{D}$, i.e. $\mathcal{C}_{\mathcal{D}}(s,t)=\mathcal{C}(s,t)$ and $s,t\in\mathcal{D}$.

[TABLE]

[TABLE]

The analysis of multipoint correlators is straightforward, leading to the general expression

[TABLE]

where the sum runs over all the permutations $\sigma$.

Computing bosonic correlators in the assumption of a gaussian ensemble in the fermionic fields leads to similar expressions, with the introduction of extra minus signs due to the fermionic Wick Theorem

[TABLE]

where now $\mathcal{F}_{\mathcal{D}}$ is the fermionic two point correlator restricted to the domain $\mathcal{D}$. Fredholm-like expressions for the bosonic correlators based on a gaussian fermionic ensemble were already knownLenard ; Korepin ; kcc14 ; mckc14 , but with the inverse formula Eq. (14) it is now possible to determine the initial conditions of the free (in terms of the fermions) evolution.

III Homogeneous quenches in the Tonks gas phase

As first application of the derived expressions we consider the quench protocol in the homogeneous case towards the Tonks gas. For the seek of concreteness, we often refer to the free thermal ensemble as the pre-quench state

[TABLE]

where the chemical potential $\mu$ is fixed requiring $\mathcal{C}(x,x)=n$, with $n$ the particle density. The steady state information is completely encoded in the initial two point fermionic correlator, that is conserved through the time evolution. While the numerical evaluation of Eq. (14) can be easily implemented, its analytical solution can be obtained only in very special cases, since the presence of the finite domain $\mathcal{D}$ prevents the possibility of diagonalizing the operators in the Fourier space. However, Eq. (14) is still feasible of proper analytical approximations in several regimes.

Short distances.—

In the short distance regime or in the small density limit, a satisfactory approximation of Eq. (14) can be obtained from the series expansion in terms of $\mathcal{C}_{\mathcal{D}}$. For example, assuming low density and that the two point bosonic correlators decays fast enough to ensure the convergence of the integrals in Eq. (10-11) for infinite domain $\mathcal{D}$, we can readily write

[TABLE]

Small temperatures.—

When the two point correlator is approximately constant $\mathcal{C}(x,y)\simeq\mathcal{C}(x,x)=n$ over the whole domain $\mathcal{D}$. This is an important case, since this condition is met in the zero temperature limit, at fixed density, of the free thermal ensemble: in the low temperature limit $\mathcal{C}(x,y)\simeq ne^{-\alpha|x-y|}$ with $\alpha=m(n\beta)^{-1}$. Thus in the small temperature limit $\alpha$ approaches zero and $\mathcal{C}$ is approximatively constant on an increasing range of distances. In this limit the following two point fermionic correlator is immediately derived

[TABLE]

This correlator greatly differs from the findings in Ref. kcc14, for the BEC case, where the decaying has been found to be purely exponential rather than algebraic. The two results are not in contradiction, since the condensation of one dimensional systems is realized only exactly at zero temperature, thus the zero temperature limit of the thermal ensemble does not reproduce the BEC state.

Large distances.—

When the distance of the two fermionic operators in the two point correlator far exceeds the typical decay length of the bosonic correlator $\mathcal{C}$, we can extract the asymptotics of the two point fermionic correlator. We start considering the string contribution

[TABLE]

where we already specialized the domain $\mathcal{D}=(x,y)$. In the limit of a very large interval, most of the contribution to the integral comes from coordinates far from the edges of the domain. Thus, in first approximation, we can compute $[\log(1+2\mathcal{C}_{\mathcal{D}})](z,z)$ as if the domain $\mathcal{D}$ were the whole real axis and diagonalize the operator in the Fourier space. This approximation gives the extensive part of the integral

[TABLE]

where $\tilde{\mathcal{C}}(k)$ is the Fourier transform of the two point bosonic correlator. The same approximation can be justified on the whole two point fermionic correlator, being reliable only in the asymptotic regime

[TABLE]

In the case of the thermal ensemble, we have

[TABLE]

where the parameters $\xi$ and $\lambda$ are

[TABLE]

The analysis of the asymptotes displays an oscillating behavior due to the imaginary part of $\lambda$, however in view of our crude approximations only the overall exponential decay $\langle\Phi^{\dagger}(y)\Phi(x)\rangle\propto e^{-(\xi+\Re(\lambda))|x-y|}$ is expected to be reliable. In Fig. 2 we plot the numerical computation of the stationary two point fermionic correlator after a quench from a thermal bosonic ensemble (16) and compare it with our analytical predictions. As already stated before, we only report the pure exponential decay. Notwithstanding, we checked the predicted period $2\pi/\Im(\lambda)$ against the numerical results and, quite remarkably, we found a very good agreement. However, the asymptotic expansion does not give a correct prediction for the oscillations’ phase.

We also analyze the momentum density distribution $\tilde{n}(k)=\int dz\,{\rm e}^{-ikz}\langle\Phi^{\dagger}(z)\Phi(0)\rangle$, which turns out to be very different from the fermionic momentum distribution at thermal equilibrium. In general, for any finite temperature, the momentum density distribution is an analytic function for all momenta. Only in the zero temperature limit, the algebraic decay of the fermionic two-point function affects the short-distance behavior of the momentum distribution which reads $\tilde{n}(k)\simeq 1-\pi|k|/4n$ (for $\beta\to\infty$). On the contrary, the large momentum behavior $\tilde{n}(k)\simeq 8n^{2}/k^{2}$ it turns out to be independent of the temperature, being fixed by the non analyticity of the fermionic correlator $\langle\Phi^{\dagger}(z)\Phi(0)\rangle\simeq n-4n^{2}|z|$ in $z=0$. Let us point out that such non-analytic behavior in the fermionic two point function is somehow universal since it only depends on the properties of the initial bosonic correlator in the vicinity of $z=0$, namely $\mathcal{C}(x,x)=n$ and $\partial_{x}\mathcal{C}(x,y)|_{y=x}=-\partial_{y}\mathcal{C}(x,y)|_{y=x}=0$.

Finally, using the steady-state of the two-point fermionic correlation function, we can numerically evaluate the stationary bosonic correlators $\langle\psi^{{\dagger}}(z)\psi(0)\rangle$ via Eq. (15). The Bose gas in the steady state is fully characterized by a two-point function which exhibits a nearly exact exponential decay for all temperatures, even for $\beta=\infty$, even tough the stationary fermionic two-point function decays algebraically. Actually, the exponential decay of the bosonic correlator in the steady state can be argued by mean of the symmetric analysis performed to extract the asymptotic decay of the two-point fermionic correlator. Interestingly, for all temperatures, we can easily extract the short distance behaviour of the bosonic two-point function which, in particular, coincides with the fermionic one, i.e. $\langle\psi^{{\dagger}}(z)\psi(0)\rangle\simeq\langle\Phi^{\dagger}(0)\Phi(0)\rangle+\partial_{z}\langle\Phi^{\dagger}(z)\Phi(0)\rangle|_{0}z=n-4n^{2}|z|$.

Lastly, we remark that even though the asymptotic behavior of the two point fermionic correlator does not suffice to determine the decay of an arbitrary bosonic correlator because of the presence of the string, this is instead possible when the string is absent. This is the case, for example, in the density-density correlator, being the density equivalently expressed in terms of either the fermions or the bosons.

[TABLE]

Thus, the density-density correlator (connected part) decays twice as fast as the two point fermionic correlator (22) and it can be computed exactly in the zero temperature limit (18).

IV Quenching in a thermodynamic trap

The expression (14) for the correlators does not require translational invariance, but only gaussianity of the initial state. Thus, we can generalize the previous calculations to the more realistic situation in which a confining potential is present and replace the Hamiltonian (2) with

[TABLE]

where $V(x)$ is a trapping potential. Differently from the homogeneous case, the two point fermionic correlator exhibits now a non-trivial time evolution.

In presence of the trap, the discrete spectrum prevents the relaxation to a steady state, that can be recovered only in a proper thermodynamic limit: for this reason, we introduce a thermodynamic length-scale $L$. Parametrizing the confining potential as $V(x)=v(xL^{-1})$, where $v$ is assumed to be a smooth function kept fixed in the thermodynamic limit, the eigenvalues approach a continuum distribution when $L\to\infty$. While considering the limit of large trap, we must also rescale the number of particles consistently and keep the total density $n=L^{-1}\int dz\,\,\mathcal{C}(z,z)$ constant.

The thermodynamic limit hugely simplifies the subsequent calculations and permits to consider a generic trap. First of all, in the thermodynamic limit the thermal initial state is described by the Local Density Approximation LDA (LDA). In fact, the two point correlator decays on a typical scale (finite for any non zero temperature) much smaller than the thermodynamic length characterizing the variation of $V$. It is convenient to change variables in the two point bosonic correlator $\mathcal{C}(y,x)\to C(s,q)$, defining $s=y-x$ and $q=(x+y)/2L$. In the LDA we have

[TABLE]

As is clear from the general expression (14), the LDA on the bosonic correlators implies LDA for the fermionic correlators as well, i.e. $\mathcal{F}(x,y)\to F(s,q)$.

Following the time evolution after the quench, we can clearly distinguish two time scales.

1. On a non-thermodynamic time scale, the system behaves as if it was locally homogeneous, therefore it exhibits gaussification of the correlators on the same time scale as the homogeneous system. In this first phase, the two point correlator does not evolve: even though the ensemble is gaussian, the steady state has not been reached yet.

2. Relaxation can be realized only on a time scale sufficient for the system to explore the whole trap, thus times $t\propto L$ must be considered. The two point correlator evolves under the Schrödinger equation, that in terms of the rescaled variables reads

[TABLE]

For any finite time $t$, in the limit $L\to\infty$ the two point correlator does not evolve. However, rewriting the differential equation in terms of the rescaled time $\tau=t/L$ and Taylor expanding the potential $v$, we get a non trivial time evolution

[TABLE]

where we assumed $\tau\ll L$ (i.e. $t\ll L^{2}$). The above equation is associated with a classical evolution for the two point correlator. This is unveiled in terms of the well known Wigner distribution wigner , that amounts to a partial Fourier transform of the correlator

[TABLE]

In this new notation, Eq. (28) can be written in the form of a classical Liouville equation in terms of $q$ and its conjugated momentum $p$

[TABLE]

where $\{\,,\}_{p,q}$ are the classical Poisson brackets and $H_{\text{cl}}$ is the classical Hamiltonian $H_{\text{cl}}=p^{2}/2m+v(q)$. In fact, the thermodynamic limit we are considering is equivalent to a semiclassical limit on the time evolution of the Wigner distribution wigner .

If the frequency of the solution of the classical equation of motion does not have a trivial dependence on the energy, the correlator in the coordinate space reaches a steady state thanks to the dephasing between the different energy shells. The steady state is of course described in terms of the energy density of the initial Wigner distribution

[TABLE]

where the energy density is defined as

[TABLE]

Based on the integrability of the free fermionic model, the steady state is expected to be described by a GGE constructed out of the modes of the trap: in Appendix B we show that the quantum mechanical prediction coincides with the above semiclassical result, as it should be.

As already mentioned, relaxation is caused by the dephasing between the different classical trajectories. While this is the common scenario for a generic trap, this is no longer true in the parabolic case: being the period of the oscillations independent from the initial conditions, dephasing is absent and the two point correlator exhibits revivals. Amusingly, notice that for a parabolic potential Eq. (27) coincides with Eq. (28) without taking the thermodynamic limit, thus the Wigner distribution evolves classically even for finite harmonic traps.

As an example, we explicitly analyze the zero temperature limit. In this limit the LDA of the initial fermionic two-point correlation function reduces to

[TABLE]

and the Wigner distribution evolves as

[TABLE]

where $\tilde{f}(p_{0})=\int ds\,{\rm e}^{-ip_{0}s}n/(1+2n|s|)^{2}$, and $(q_{\tau}(p_{0}),\,p_{\tau}(p_{0}))$ is the solution of the classical equation of motion with initial condition $(0,p_{0})$. Therefore, exploiting the normalization of the Dirac delta function, we can integrate over the final momentum $p$, obtaining for the correlation function

[TABLE]

where the sum runs over all the real roots $p_{j}(q,\tau)$ of the equation $q_{\tau}(p_{j})=q$. In Fig. 3 we plot the Wigner distribution at different times $\tau$, and for two representative confining potentials, namely $v(x)=x^{2}/2$ and $v(x)=x^{4}/4$ (where, in the following, we fixed $m=1$ and $n=1$). In the harmonic case, the Wigner distribution exhibits a trivial periodic dynamics and the particle density never reaches a stationary profile, continuously oscillating between $\delta(q)$ and $\tilde{f}(q)/2\pi$.

On the contrary, in the anharmonic case, the number of roots of the equation $q_{\tau}(p_{j})=q$ increases with time and the root density approaches a continuous distribution; the system therefore relaxes to the stationary state whose correlation function is given by Eq. (31).

In Fig 4 we compare the steady-state predictions with the correlation function evaluated numerically using Eq. (IV), in the presence of the anharmonic confining potential $v(x)=x^{4}/4$. In particular, we plot the particle density versus $q$ as well as the two-point correlator as a function of $s=y-x$ for fixed $q=0$ and $q=2$. As expected, for large time the numerical data approach the stationary profiles.

V Conclusion

In this work we studied the equilibration mechanism of 1D Bose gases in the Tonks-Girardeau limit starting with initial states that are gaussian in the bosonic fields, in the homogeneous case as well as the general inhomogeneous case of an arbitrary confining trap. Despite the noninteracting (in terms of the fermions) dynamics, this problem is nontrivial due to the nonlocal nature of the Jordan-Wigner transformation which links the bosonic and fermionic fields. By solving the problem of extracting fermionic correlations from the bosonic ones in a Gaussian state and vice versa, we derived exact closed formulas for the steady state correlations and analysed them both numerically and analytically in several asymptotic limits.

When a confining potential is considered, the relaxation of local observables can be attained only in a proper thermodynamic limit, where the local density approximation is applicable. In this limit, the evolution is separated on two time scales: initially the system behaves as if it were homogeneous with a subsequent gaussianification. Thereafter, the particles explore the whole trap on a thermodynamic time scale with an emergent semiclassical evolution for the two point fermionic correlator, in terms of which all the expectation values of local observables are determined.

VI Acknowledgments

This work was supported by the European Union’s Horizon 2020 research and innovation program under the Marie Sklodowska-Curie Grant Agreement No. 701221 (M.C.). SS acknowledges support from the A*MIDEX project Hypathie (no. ANR-11-IDEX- 0001-02) and from the LIA - LYSM (Laboratoire Ypatia des Sciences Mathématiques).

Appendix A Correlators of fermions on the initial state

The key point in the analysis of our quench protocol is the mapping between bosonic and fermionic correlators, achieved by mean of the Jordan Wigner transformation (6). This mapping passes through the hard core bosons, whose rigorous treatment requires a lattice regularization. This appendix is devoted to this issue, in particular we show that even though fermionic fields are defined in terms of the hard core bosons, their correlation functions can be computed as if the hard core bosons were standard bosonic fields (after the expressions have been normal ordered). In the case when the initial state is chosen to be the Bose-Einstein condensate the explicit check has already been performed in Ref. kcc14, , however the same conclusions hold even for more general states. Through a simple generalization of the following, it is possible to show the inverse relation, i.e. that bosonic correlators are expressed in terms of fermionic correlators exchanging the role of bosons and fermions in Eq. (6).

Here we test expectation values on pure states with a fixed number of particles, but the generalization over more general density matrices (as the thermal ensemble) is trivial. Consider then a generic state $|S\rangle$ with $n$ particles

[TABLE]

where $W_{n}$ is the symmetric (normalizable) wavefunction, that we will assume to be smooth. We introduce a lattice space $\Delta$ and discretize the space $x\to\Delta j$, the continuous bosonic fields are replaced with the lattice counterpart $\psi(\Delta j)\to\Delta^{-1/2}\psi_{j}$

[TABLE]

where the lattice regularized fields follow standard commutation rules $[\psi_{j},\psi_{j^{\prime}}^{\dagger}]=\delta_{j,j^{\prime}}$. To be precise, the lattice wavefunction differs from the continuous one for terms $\mathcal{O}(\Delta)$, but we drop them since they are inessential in the limit in which we are interested in. The hard core bosons are obtained from the standard bosons by mean of a projector that excludes the presence of more than one particle at each lattice site

[TABLE]

where $P_{j}=|0_{j}\rangle\langle 0_{j}|+|1_{j}\rangle\langle 1_{j}|$. By mean of a straightforward application of this definition, the hardcore bosonic fields are immediately shown to satisfy the following commutation rules

[TABLE]

[TABLE]

With the hard core bosons we define the fermions on the lattice through the Jordan Wigner transformation

[TABLE]

and finally, in the $\Delta\to 0$ limit, the continuous fermionic fields $\Phi_{j}\to\Delta^{1/2}\,\Phi(\Delta j)$. Consider now an arbitrary normal ordered correlator of the hard core bosons

[TABLE]

where, because of normal ordering, there are not two creation (annihilation) operators with the same index. Moreover, from the sum over the indexes $j$ must be excluded all the configurations where an index $j$ is equal to an index $i$ or $i^{\prime}$. However, taking the continuous limit, this exclusion requirement can be neglected, since the forbidden region has a zero measure in the integral and we readily obtain

[TABLE]

Consider now the fermionic correlators. For the seek of simplicity we analyze the two point case (the extension to the general case is obvious)

[TABLE]

where without loss of generality we assumed $j>j^{\prime}$. Above we rely on a useful identity that permits to normal order the string (41)

[TABLE]

We now take the continuous limit of Eq. (44), letting $\Delta\to 0$ but keeping $\Delta j=x$ and $\Delta j^{\prime}=y$. Expanding the normal ordered exponential we get

[TABLE]

Aiming to the continuum limit we use (43) in each term and replace the sums with the proper integrals. The resulting expression is then simply resummed leading to

[TABLE]

The final step amounts to replace the l.h.s. with the continuous counterpart of the fermionic fields, leading to the desired identity

[TABLE]

Appendix B Quenching in the thermodynamic trap

In Section IV we considered the quench protocol in a thermodynamic trap and we showed the emergence of a classical picture for the Wigner distribution, with the subsequent classical relaxation to a steady state. On the other hand, based on the integrability of the model, we expect relaxation to a GGE constructed out of the modes of the trap. This appendix is devoted to showing the equivalence of the two descriptions, i.e. that the GGE predictions reduce, in the thermodynamical limit, to the classical steady state (31-32). This consistency check can be achieved thanks to the fact that the WKB approximation griffiths becomes exact in the TDL limit we are investigating.

We start presenting a brief review of the WKB method sufficient for our purposes, then use it to compute the GGE predictions in the thermodynamic limit.

B.1 WKB description of the thermodynamic trap

Consider the one particle Shroedinger equation in one dimension:

[TABLE]

The trapping potential is assumed in the form of Section IV, i.e. $V(x)=v(xL^{-1})$. For simplicity, we assume $v$ to be symmetric and with an unique minimum in $x=0$. In the limit of small variation of the potential the wavefunction can be approximated through the WKB solution. Let $\xi_{n}>0$ be the classical turning point identified by the equation $E_{n}=v(\xi_{n})$ ($E_{n}$ is the energy of the $n^{\text{th}}$ energy level), then the wavefunction oscillates inside the trap ($|x|<L\xi_{n}$)

[TABLE]

and is exponentially damped outside of it. For example, on the right hand side $x>L\xi_{n}$ we have

[TABLE]

In Eq. (50) $\mathcal{N}_{n}$ is a normalization constant to ensure the eigenfunction to be normalized to unity. The energy levels satisfy the WKB quantization condition

[TABLE]

As we already said, the WKB approximation is reliable when the variation of the potential is much smaller than the energy difference between $E$ and the potential, more quantitatively this amounts to ask

[TABLE]

Replacing $V$ with the rescaled potential $v$ (and using the rescaled variable $q=xL^{-1}$) we see that the WKB approximation becomes exact in the TDL

[TABLE]

with the exception of the turning points, i.e. the zeros of the denominators of the above expressions. However, calling $\delta q$ the distance from the turning point, the region of validity of the WKB approximation is readily estimated as

[TABLE]

Thus, in the rescaled variable, the region in which the WKB approximation fails is shrunk to zero as $\delta q\propto L^{-2/3}$. Some care must be used stating that the WKB becomes exact in the thermodynamic limit, since in terms of the original variable $x=Lq$ the region where WKB fails does not vanish, but rather diverges as $L^{1/3}$. Indeed, even though the single WKB eigenfunction is not reliable in a neighborhood of the turning point of increasing size, the non extensive scaling of the latter permits, ipso facto, a direct use of the WKB eigenfunctions in several expressions.

Before passing to analyze the GGE, we need to properly normalize to unity the wavefunction: this task is hugely simplified in the thermodynamic limit, since only the WKB solution inside the trap is needed. To understand this point consider the solution outside of the trap, let’s choose the right side for simplicity. Using the rescaled coordinate $q=L^{-1}x$ the wavefunction is exponentially damped as

[TABLE]

In the limit $L\to\infty$ the above exponential vanishes, meaning that the region where the wavefunction is not zero outside of the trap does not have a thermodynamic extent. Using this fact (and the exactness of the WKB approximation in the rescaled coordinate), the normalization constant $\mathcal{N}_{n}$ in (50) is readily written as

[TABLE]

where we neglected subextensive contributions to $\mathcal{N}^{-2}_{n}$, $\rho(E)$ is simply the volume of the energy shell in the classical phase space

[TABLE]

B.2 The excitation density in the thermodynamic limit

Armed with the WKB approximation, we can start constructing the thermodynamic limit of the GGE associated with the trap. As first step, we extract the excitation energy from the two point fermionic correlator in the assumption the latter satisfied the LDA approximation, i.e. $\langle\Phi^{\dagger}(x)\Phi(y)\rangle=\mathcal{F}(x,y)\to F(s,q)$ where $s=x-y$ and $q=(x+y)/2L$. $F$ does not contain any further dependence on the thermodynamic length $L$ and decays fast enough for large $s$ to justify the forthcoming approximations.

Denoting as $\eta^{\dagger}_{n},\,\eta_{n}$ the creation/annihilation operators of the modes of the trap, the number of particles in the energy level $\langle\eta^{\dagger}_{n}\eta_{n}\rangle$ can be readily extracted from the two point fermionic correlator and the eigenfunctions

[TABLE]

At the price of neglecting corrections vanishing in the thermodynamic limit, we can replace in the above the WKB eigenfunctions and the LDA approximation for the correlator. Using (50) with the truncation of the wavefunction outside of the trap, through simple trigonometric identities we get

[TABLE]

Thanks to the LDA approximation for the fermionic correlator and under the assumption that decays in the relative distance $x-y$, the above can be further simplified. In particular, the second term gives contributions suppressed in the thermodynamic limit and the first line can be approximated to

[TABLE]

The above is clearly a smooth function of the energy $\langle\eta^{\dagger}_{n}\eta_{n}\rangle=\epsilon(E_{n})$, where

[TABLE]

and $\tilde{F}$ is defined in Eq. (29). Notice that $\epsilon$ coincides with the classical expression (32).

B.3 GGE predictions in the thermodynamic limit

The GGE density matrix of the free fermion gas is defined over the modes of the trap $\rho_{\text{GGE}}=e^{-\sum_{n}\lambda_{n}\eta^{\dagger}_{n}\eta_{n}}/\mathcal{Z}$, in particular this ensemble is gaussian and it is completely determined by the two point correlator

[TABLE]

where the GGE expectation value of the occupancy $\eta^{\dagger}_{n}\eta_{n}$ is equal to the expectation value on the initial state, that reduces to (62) in the TDL. Computing (63) in the thermodynamic limit simply amounts to replace the eigenfunctions with the WKB expressions and the mode occupation with (62). In the TDL limit the energy levels become dense, thus we would like to replace the sum over the discrete levels with an integral over the energy, being the density of levels at fixed energy extracted from the quantization rule (52). However, the energy eigenfunctions are not smooth functions of the energy, since they are even (odd) if the quantum number is even (odd), as it is clear from (50). Therefore, we first split the sum over the even and odd eigenfunctions, then take the continuum approximation of each one. Doing such an operation and recollecting the terms in an unique integral we readily find

[TABLE]

Since $\epsilon(E)$ is a smooth function, the above integral decays when $x$ and $y$ are pulled far apart and in the limit of large $L$ we can safely approximate

[TABLE]

where, as usual, $s=x-y$ and $q=(x+y)/2L$. The above expression matches with the classical computation (31-32), as it should be.