One-dimensional Bose gas driven by a slow time-dependent harmonic trap

TL;DR

This paper investigates the dynamics of a one-dimensional hard-core boson cloud under slow, time-dependent harmonic trap release, providing analytical solutions and exploring effects like log-periodic modulation during expansion.

Contribution

It offers exact and approximate analytical solutions for the cloud's evolution under various trap ramp protocols, including a novel case with a time-inverse amplitude variation.

Findings

Derived a typical length scale related to the ramp

Compared perturbative and exact solutions for linear ramp

Discovered log-periodic modulation in cloud expansion

Abstract

We consider the unitary time evolution of a one-dimensional cloud of hard-core bosons loaded on a harmonic trap potential which is slowly released in time with a general ramp $g(t)$. After the identification of a typical length scale $\ell(t)$, related to the time ramp, we focus our attention on the dynamics of the density profile within a first order time-dependent perturbation scheme. In the special case of a linear ramp, we compare the first order predictions to the exact solution obtained through Ermakov-Lewis dynamical invariants. We also obtain an exact analytical solution for a cloud released from a harmonic trap with an amplitude that varies as the inverse of time. In such situation, the typical size of the cloud grows with a power law governed by an exponent that depends continuously on the initial trap frequency. At high enough initial trap amplitude, the exponent acquires an imaginary part that leads to the emergence of a log-periodic modulation of the cloud expansion.

Equations

$$H=-\frac{J}{2}\sum _{j=-L/2}^{L/2}[a_{j+1}^{†}{a}_j+h.c.]+\frac{U}{2}\sum _{j=-L/2}^{L/2}{n}_j(n_j-1)+\sum _{j=-L/2}^{L/2}{V}_j(t)n_j$$

$$H=-\frac{1}{2}\sum _{j=-L/2}^{L/2}[b_{j+1}^{†}{b}_j+h.c.]+\sum _{j=-L/2}^{L/2}{V}_j(t)n_j$$

$$c_j^{†}=\prod _{i<j}(1-2b_i^{†}{b}_i)b_j^{†}$$

$$H=\sum _{i,j=-L/2}^{L/2}{c}_i^{†}{A}_{i,j}(t)c_j,$$

$$A_{i,j}(t)\equiv V_i(t){\delta }_{i,j}-\frac{1}{2}({\delta }_{i,j+1}+{\delta }_{i+1,j}).$$

$$H=\sum _{q=0}^L{ϵ}_q(t){\eta }_q^{†}(t){\eta }_q(t),$$

$${\eta }_q^{†}(t)=-\sum _{i=-L/2}^{L/2}{\psi }_q(i,t)c_i^{†}$$

$$V_j(t)=|g(t)|j^2-\mu$$

$$A(t){\psi }_q(t)={ϵ}_q(t){\psi }_q(t).$$

$${\psi }_q(x_0\pm a,t)\simeq {\psi }_q(x_0,t)\pm a{\partial }_x{\psi }_q(x_0,t)+\frac{a^2}{2!}{\partial }_x^2{\psi }_q(x_0,t),$$

$$\frac{1}{2}{\partial }_x^2{\psi }_q(x,t)+[{ϵ}_q(t)+1-V(x,t)]{\psi }_q(x,t)=0,$$

$$\frac{1}{2}\Big{(}-\partial_{x}^{2}+\omega^{2}(t)x^{2}\Big{)}\psi_{q}(x,t)=\epsilon_{q}(t)\psi_{q}(x,t)\;,\qquad\omega(t)=\sqrt{2|g(t)|}$$

$${\psi }_q(x,t)=\sqrt{\frac{\sqrt{\omega (t)}}{2^{q}q!\sqrt{\pi }}}{e}^{-\frac{\omega (t)}{2}{x}^2}{\mathrm{He}}_q(x\sqrt{\omega (t)}),{ϵ}_q(t)=\omega (t)(q+\frac{1}{2}),$$

$${\psi }_q(x,t)={\ell }^{-1/2}{\chi }_q\left(\frac{x}{\ell }\right),{\chi }_q(u)=\frac{1}{\sqrt{2^{q}q!\sqrt{\pi }}}{e}^{-\frac{1}{2}{u}^2}{\mathrm{He}}_q(u)$$

$${ϵ}_q(t)={\ell }^{-2}(q+\frac{1}{2})$$

$$\delta \mu (x,t)\equiv \mu (x,t)-{\mu }_c\approx -|g(t)|x^w$$

$$\ell (t)\propto |\delta \mu (\ell (t),t){|}^{-\nu }$$

$$\ell (t)\propto |g(t){|}^{-{\nu }_g},{\nu }_g=\frac{\nu }{1+\nu w}$$

$$\phi (x,t)\propto {\ell }^{-x_{\phi }}\tilde{\Phi }(\frac{x}{\ell }),$$

$$\phi (x,L,t)\propto {\ell }^{-x_{\phi }}\tilde{\Phi }(\frac{x}{\ell },\frac{L}{\ell }).$$

$$|{\mathcal{S}}_0(t_0)⟩=\prod _{q=0}^{N-1}{\eta }_q^{†}(t_0)|0⟩,$$

$$E_0(t_0)=\sum _{q=0}^{N-1}{ϵ}_q(t_0).$$

$${\rho }^{ad}(i,t)$$

$${\rho }^{ad}(x,t)={\ell }^{-1}{f}^{ad}\left(\frac{x}{\ell }\right),f^{ad}(u)=\sum _{q=0}^{N-1}|{\chi }_q(u){|}^2,$$

$${\rho }^{LDA}(x,t)=\frac{1}{\pi }\arccos (V(x,t)-{ϵ}_{N-1}(t)).$$

$${\rho }^{LDA}(x,t)=\frac{1}{\pi }\arccos \left(1-\frac{1}{2{\ell }^2}\left[2N-1-{\left(\frac{x}{\ell }\right)}^2\right]\right).$$

$${\rho }^{LDA}(x,t)\simeq \frac{1}{\pi \ell }{\left(2N-1-{\left(\frac{x}{\ell }\right)}^2\right)}^{1/2}\underset{N\gg 1}{\simeq }\frac{\sqrt{2N}}{\pi \ell }{\left(1-\frac{1}{2}{\left(\frac{x}{\ell \sqrt{N}}\right)}^2\right)}^{1/2}$$

$${\ell }_N=\ell \sqrt{2N-1}\underset{N\gg 1}{\simeq }\ell \sqrt{2N}.$$

$${\rho }^{LDA}(x,t)\simeq \frac{N}{{\ell }_N}{f}^{LDA}\left(\frac{x}{{\ell }_N}\right),f^{LDA}(u)=\frac{2}{\pi }\sqrt{1-u^2}.$$

$$\ket{\mathscr{S}(t)}\simeq e^{-i\int_{t_{0}}^{t}dt^{\prime}\,E_{0}(t^{\prime})}\Big{[}\ket{\mathscr{S}_{0}(t)}+\sum_{p=N}^{L}\sum_{k=0}^{N-1}a_{p,k}(t)\,\ket{\mathscr{S}_{0}[\not{k},p\,](t)}\Big{]}$$

Full text

1 Introduction

Advances in ultracold atomic gases have led to the possibility of realizing low-dimensional model systems, e.g.[1, 2, 3, 4, 5, 6, 7]. This has opened the road to probe experimentally theoretical predictions on the dynamical aspects related to non-equilibrium effects, see [8, 9] for recent reviews. In particular, quantum gases parametrically driven through a quantum phase transition have played a central role [8] for testing new ideas related to the thermalization (or the absence of it) of such gases. A paradigmatic model in this context is the Bose-Hubbard (BH) model [10, 11], describing interacting bosons on a lattice, and exhibiting Superfluid-to-Mott Insulating (SF-MI) quantum phase transitions [11, 12]. To mention just a few examples of the wide literature dedicated to such so called quantum quenches, we may quote J. M. Zhang et al. [13, 14] that have studied the quantum quench driven by the on-site interaction or J. Dziarmaga et al. [15, 16, 17] that have investigated the case of a time-dependent hopping magnitude, focusing on the loss of adiabaticity and on Kibble-Zurek scaling regims [22, 23].

In this context, our aim here is to investigate close to the SF-MI transition the dynamics of a one dimensional cloud of bosons driven by the slow release of an inhomogeneous confining potential. We focus our attention on the case of a harmonic potential $V(x)\sim g(t)x^{2}$ for a given time ramp $g(t)$. We consider in particular the low-density regime of the Bose gas in the limit of high repulsive interactions for which the Tonks-Girardeau model [18, 19, 20, 21] is a good effective description. On one hand this problem can be handled within the framework of quantum quenches ideas where off-equilibrium behaviors are related to the breakdown of adiabaticity near the phase transition. The dynamical behavior after a sudden release of the trap has been studied so far quite extensively, see for example [24, 25, 26, 27]. In the case of a slow driving, which is the situation considered in this work, the system can be investigated perturbatively around the adiabatic evolution [9]. The departure from the equilibrium can be quantified studying the first off-equilibrium corrections of the density profile or counting, à la Kibble-Zurek, the number of excitations generated during the quench [28]. On the other hand, another way to approach the problem is to find dynamical invariants associated to the time-dependent system, which basically after diagonalization reduces to a set of time-dependent harmonic oscillators, see e.g. [29, 30]. For the special case of a linear ramp $g(t)$ one is then able to solve explicitly the associated non-linear Pinney differential equation [31] and rebuild explicitly the many-body wave function from which exact solutions for the physical observables and in particular for the particle density are available [32, 33].

The paper is organized as follows: in the next section we present the model and its mapping to a Fermi system in the limit of hard-core bosons. After the explicit diagonalization of the instantaneous Hamiltonian which is performed in section 3 we present the results obtained for the density profile in the first order time-dependent perturbation theory framework in section 4. The dynamical invariants approach is presented in section 5 where the exact result for the density profile, obtained in the case of a linear ramp, is compared to the first order perturbative one. The release of the harmonic trap with a frequency that varies as the inverse of the time is also considered there and solved explicitly. It is shown in particular that the bosonic cloud expands with a power-law behaviour which exponent is a continuous function of the initial amplitude of the trap. At high initial amplitudes, the exponent becomes complex and leads to the appearance in the expansion of the cloud to a log-periodic modulation in time of square-root growing law. Finally, a brief summary and conclusions are given in the last section.

2 The model

The Hamiltonian of the one dimensional Bose-Hubbard model, which describes a set of bosons living on a lattice with repulsive on-site interaction $U$ and submitted to an external time-dependent potential $V(t)$ is given by

[TABLE]

where the creation and annihilation operators $a^{\dagger},a$ satisfy the usual canonical bosonic algebra $[a_{j},a^{\dagger}_{k}]=\delta_{j,k}$, $[a_{j},a_{k}]=[a^{\dagger}_{j},a^{\dagger}_{k}]=0$ and where $n_{j}=a^{\dagger}_{j}a_{j}$ stands for the occupation number at site $j$. The first term proportional to $J$ describes the kinetic part of the system and in the following we will set the hopping amplitude $J=1$. In the hard core limit, that is for a very large repulsive interaction $U/J\gg 1$, the Bose-Hubbard Hamiltonian reduces to

[TABLE]

with a new set of operators $b^{\dagger},b$ that still satisfy the bosonic algebra for different sites but fulfills the on-site anti-commutation relations $\{b^{\dagger}_{j},b_{j}\}=1$, $\{b_{j},b_{j}\}=\{b^{\dagger}_{j},b^{\dagger}_{j}\}=0$, which prevent a double occupancy of a given site (the occupation operator is $n_{j}=b^{\dagger}_{j}b_{j}$). Obviously, the Pauli raising and lowering operators $\sigma_{j}^{+}$ and $\sigma_{j}^{-}$ realize the algebra generated by the operators $b_{j}^{\dagger}$ and $b_{j}$.

The standard procedure to diagonalize the Tonks-Girardeau (TG) Hamiltonian (2) is first to fermionize it through the Jordan-Wigner transformation [34], mapping the operators $b$ with mixed bosonic and fermionic characters to simple fermionic $c$ operators:

[TABLE]

and the associated mapping for the adjoint annihilation operators $c_{j}$. Under this transformation, the TG Hamiltonian reduces to a spinless tight-binding Fermi system. Explicitly, for a finite size lattice with open boundary conditions one has

[TABLE]

where we have introduced the matrix $A(t)$:

[TABLE]

Notice that the occupation number operator $n_{j}=b^{\dagger}_{j}b_{j}=c^{\dagger}_{j}c_{j}$. At a given time $t$ the quadratic form (4) is readily diagonalized through a unitary transformation reducing $(\ref{A})$ into a free theory:

[TABLE]

where the $\epsilon_{q}(t)$ are the single particle energies and where the diagonal Fermi operators $\eta^{\dagger}$, $\eta$ are defined through

[TABLE]

and the associated relation for the adjoints $\eta_{q}(t)$. Notice here that the time $t$ appears in these expressions as a simple parameter. The Bogoliubov coefficients $\psi_{q}(i,t)$ satisfy the orthonormality condition $\sum_{i}\psi^{*}_{q}(i,t)\psi_{p}(i,t)=\delta_{qp}$ which in turn implies the canonical anti-commutation algebra $\{\eta^{\dagger}_{q}(t),\eta_{p}(t)\}=\delta_{qp}$, $\{\eta^{\dagger}_{q}(t),\eta^{\dagger}_{p}(t)\}=\{\eta_{q}(t),\eta_{p}(t)\}=0$. The minus sign in (7) is irrelevant and set for further conveniences.

The time-dependent potential $V(t)$ confining the Bose gas is

[TABLE]

where the shift $\mu$, with $-1<\mu<1$, can be interpreted as a chemical potential and the time-dependent amplitude $g(t)$ is assumed to be a slowly-varying function (this will be made precise later).

3 Instantaneous diagonalization and adiabatic evolution

3.1 Instantaneous diagonalization

The instantaneous single particle energy spectrum and the corresponding eigenvectors are derived through the diagonalization of the matrix $A(t)$ (notice that here the time variable $t$ is just a parameter):

[TABLE]

Numerical exact diagonalization are easily performed for such a problem (see figure 1). However, in the thermodynamic limit, in which the lattice site $i$ is replaced by a continuous variable $x=ai$, and expanding

[TABLE]

the eigenvalue problem $(\ref{eigendiscrete})$ reduces to

[TABLE]

where $V(x,t)=|g(t)|x^{2}-\mu$ is the continuum limit of the lattice potential of the Eq. $\eqref{trap}$. Choosing $\mu=-1$ (which corresponds in the absence of the trapping potential to the transition point between the trivial Mott phase with zero density for $\mu<-1$ and the superfluid phase for $|\mu|<1$) one obtains the stationary Schrödinger equation for the $1d$ quantum harmonic oscillator

[TABLE]

with energies $\epsilon_{q}$ and corresponding eigenfunctions $\psi_{q}$. Explicitly the solution is

[TABLE]

where $q\in\mathbb{N}$ and ${\rm He}_{q}$ denotes the $q^{th}$ Hermite polynomial with physical normalization implied from the normalization of the single particle wave functions $\psi_{q}(x,t)$. In figure 1 we show the exact low-lying single particle rescaled energies as a function of the trap amplitude $|g(t)|$ and compare them to the continuum limit prediction given above. As the system size is increased the agreement gets better and better.

The associated lowest eigenvectors for different sizes are shown in figure 2.

The frequency $\omega(t)$ defines an instantaneous length scale in the problem which is given by $\ell(t)\equiv(\omega(t))^{-1/2}$. In terms of this length scale, the solution is

[TABLE]

for the eigenstates, and

[TABLE]

for the energy spectrum, where the time-dependence of $\ell(t)$ has been implicitly considered. This scaling form for the eigenfunctions and energies is in agreement with general time-dependent trap-size scaling arguments. Indeed, for a driven inhomogeneous quantum quench through the MI-SF critical point, the inhomogeneous control parameter

[TABLE]

induces a finite-length scale $\ell(t)$ [36, 35, 37, 28, 38]. This typical length scale can be derived self-consistently from the space dependence of the quantum deviation parameter by

[TABLE]

where $\nu$ is the correlation critical exponent. The solution of this equation is

[TABLE]

The scaling behavior of a local quantity $\varphi(x,t)$ with scaling dimension $x_{\varphi}$ is expected to be

[TABLE]

where $\tilde{\Phi}$ is a scaling function. At the MI-SF transition, where the critical exponents are $\nu=1/2$ for the correlation length and $z=2$ for the dynamics, fixing $w=2$ for a parabolic trap one recovers the scaling forms (14) and (15) with $\nu_{g}=1/4$ and $x_{\psi}=1/2$ for the scaling dimension of the single particle wave function. Notice that for a finite size system, the size $L$ of the system itself becomes a scaling field and one expects that a given quantity depends on the two length scales $\ell(t)$ and $L$, such that

[TABLE]

For $L\ll\ell(t)$ it is expected that the scaling relation $(\ref{scalL})$ matches the ordinary finite-size scaling behavior $\varphi\propto L^{-x_{\varphi}}$ while for $L\gg\ell(t)$ the system becomes independent on the lattice size $L$ and matches the infinite volume behaviour. The thermodynamic limit is therefore taken as the limit $L\rightarrow\infty$, $\ell\rightarrow\infty$ with the ratio $L/\ell^{2}$ (or equivalently $L^{2}g$) fixed such that finite size corrections are avoided.

3.2 Adiabatic evolution of the density

At the initial time $t_{0}$ the system is prepared in the $N$-particles ground state of (2) which is given by

[TABLE]

where $|0\rangle$ is the vacuum state such that $\eta_{q}(t_{0})|0\rangle=0$ $\forall q$, since all the single particle energies are positive. The energy associated to the initial state $|\mathscr{S}(t_{0})\rangle$ is thus simply given by

[TABLE]

For a very slow variation of the confining potential, it is expected that the system adapts itself to the instantaneous Hamiltonian $(\ref{TGH})$ and evolves remaining in the instantaneous ground state $|\mathscr{S}_{0}(t)\rangle=\prod_{q=0}^{N-1}\eta^{\dagger}_{q}(t)|0\rangle$. The evolution of the particle density is thus expected to be given by the adiabatic density

[TABLE]

or in the thermodynamic limit by the scaling form

[TABLE]

where the functions $\chi_{q}$ are defined in (14). In figure 3 we show the convergence of the adiabatic density profile obtained from exact numerical diagonalization toward the thermodynamical limit expression for small particle numbers $N$. In the large size limit the exact numerical results match perfectly the analytical expression (24). At smaller sizes the finite size corrections to the excitation spectrum and to the corresponding eigenvectors lead to quite a discrepancy between the thermodynamic limit result and the actual finite size density profile.

In figure 4 we show the adiabatic density for $N=3$ bosons for three different values of the amplitude $g(t)$.

For a large number of bosons the adiabatic density profile matches its Local Density Approximation (LDA) limit. The LDA is obtained assuming that for each instant of time $t$ and around each (coarse-grained) point $x$ there is a local flat band of excitations with dispersion $\varepsilon(x)=-\cos q_{F}(x)+V(x,t)$. Locally the single particle band is filled up to the global Fermi level given here by $\varepsilon_{F}=\varepsilon_{N-1}(t)$ and the LDA $\rho^{LDA}(x)$ is deduced from the associated local Fermi momentum $q_{F}(x)=\pi\rho^{LDA}(x)$ leading to [24, 39, 40]

[TABLE]

This is shown in figure 5 for three different values of $N$. The exact numerical profile, as obtained from exact diagonalization, is compared to the thermodynamic limit (24) and LDA (25) results. One can recover the scaling form (24) from (25) in the limit $N\ll\ell^{2}$. Indeed, using $\ell=(\omega(t))^{-1/2}=(2|g(t)|)^{-1/4}$ and the expression of the potential $V(x,t)=1+|g(t)|x^{2}$ one has

[TABLE]

For $N/\ell^{2}\ll 1$, expanding the $\arccos$ function to the leading order in $N^{1/2}$, one obtains a semi-circle law

[TABLE]

in agreement with the scaling form (24). The support of the adiabatic density profile is in $[-\ell_{N},\ell_{N}]$ with a number of particles typical length scale

[TABLE]

In terms of that typical length scale $\ell_{N}$, the adiabatic profil takes the scaling form

[TABLE]

4 First order correction to the adiabatic evolution

4.1 First order correction to adiabaticity

The first order correction to the adiabatic evolution, starting from an initial $N$ particles ground state $\ket{\mathscr{S}_{0}(t_{0})}$ as defined in (21), is given by 111see the Appendix A of the work in Ref. [28]:

[TABLE]

where $\ket{\mathscr{S}_{0}[\not{k},p\,](t)}\equiv\eta^{\dagger}_{p}(t)\,\eta_{k}(t)\ket{\mathscr{S}_{0}(t)}$ is the instantaneous ground state in which a particle has been promoted from the lowest levels $k=0,\dots,N-1$ ($\not{k}$ denotes a vacancy in the position $k$) toward higher ones $p=N,\dots,L$ (since the particles are fermionic in nature double occupancy of a state is forbiden). The transition amplitude $a_{p,k}(t)$ is given by

[TABLE]

where

[TABLE]

is the instantaneous transition amplitude of the perturbation

[TABLE]

In the scaling limit $g\rightarrow 0$, $L\rightarrow\infty$, the first order transition amplitude $a_{p,k}(t)$ can be computed analytically using the thermodynamic limit expressions (14) and (15). After a straightforward computation one obtains

[TABLE]

Starting from the initial $N$ particles ground state $\ket{\mathscr{S}_{0}(t_{0})}$, as seen from (30) at the first order in perturbation only the lowest energy levels are activated: $(p,k)=(N+1,N-1)$ and for $N>1$, $(p,k)=(N,N-2)$. A crude approximation $p\simeq k\simeq N\gg 1$ shows that the transition amplitude $(\ref{a})$ is of order $N\,\ln\frac{\ell(t)}{\ell(t_{0})}$ . Consequently, the order of the approximation (30) is

[TABLE]

with $\delta g(t)\equiv g(t)-g(t_{0})$.

4.2 Particle density

The particle density at site $i$ and time $t$ is given by

[TABLE]

where the two-point function

[TABLE]

With the expansion (34) at the leading order in $\varepsilon$ the two-point function is given by

[TABLE]

The first term of the two-point function gives the adiabatic contribution (23). The deviation to the adiabatic density is thus expressed as

[TABLE]

In the scaling limit $g\rightarrow 0$, $L\rightarrow\infty$, using (14) and (34) we obtain

[TABLE]

with $C_{1,-1}=0$ by convention. Using recursion relations of the Hermite polynomials the density deviation $\delta\rho(x,t)$ can be rewritten as

[TABLE]

where $F_{N}(u)$ is given by

[TABLE]

A plot of the density profile in the quasi-adiabatic evolution for low values of $N$ is shown in figure 6.

At large particle number $N$, the Hilbert-Hermite functions $\chi_{N}(u)$ take significant values only in the region $|u|<\sqrt{2N}$ where the zeros of the Hermite polynomials are located. Outside that region the Hilbert-Hermite functions decay exponentially fast. In the limit $N\gg 1$, it has been shown in [41] that the Hermite polynomials have the asymptotic representation for $\theta\in]-\pi/2,\pi/2[$ given by

[TABLE]

with the phase

[TABLE]

With Stirling formula $N!\simeq\sqrt{2\pi N}(\frac{N}{e})^{N}$ and the asymptotic representation given above, the Hilbert-Hermite functions, in the limit $N\gg 1$, take the form

[TABLE]

Using this, one has for (42) at large $N$ the asymptotic expression

[TABLE]

The $\sin^{2}h_{N}(u)$ gives a widely oscillating term and taking its average, $\sin^{2}h_{N}(u)\sim 1/2$, one finally obtains a scaling form for the deviation $\delta\rho(x,t)$ as a function of the scaling variable $x/\ell_{N}$:

[TABLE]

where we have set the small parameter $\varepsilon$ associated to the first order correction to

[TABLE]

This behavior is shown in figure 7 and figure 8 for the associated scaling function at large number $N$.

Notice that using the asymptotic expression (45) one has after averaging the $\cos^{2}$ term

[TABLE]

where $\theta(q)$ is the Heaviside function. The adiabatic density is thus

[TABLE]

which is nothing but the LDA semi-circle law (29).

5 Exact Ermakov-Lewis evolution

5.1 Dynamical invariant approach

The Ermakov-Lewis approach [29, 30] based on the identification of dynamical invariants is an alternative way that leads to exact results for particular ramps $g(t)$ . Let us first reconsider the single-particle Schrödinger equation:

[TABLE]

with $\omega(t)=\sqrt{2|g(t)|}$ and suppose that the initial condition at $t_{0}$ is an eigenstate $\phi_{q}(x,t_{0})=\psi_{q}(x,t_{0})$ (see (13)) of the harmonic oscillator with initial pulsation $\omega(t_{0})\equiv\omega_{0}$. According to the dynamical invariant approach, the time-evolved single-particle wave function $\phi_{q}(x,t)$ can be expressed as [29, 30, 31, 32]

[TABLE]

where $\zeta(t)$ is the solution of the non-linear differential equation

[TABLE]

with initial conditions $\zeta(t_{0})=1$ and $\overset{\centerdot}{\zeta}(t_{0})=0$. The problem is thus reduced to solving this differential equation given the time-dependent protocol $g(t)$. For instance, setting a linear ramp

[TABLE]

with time-scale $1/\alpha$, a solution for $\zeta$ is explicitly known in terms of Airy functions, see for example [39, 42], and is plotted in figure 9 as an illustration. The time-evolution of the density profile for the linear ramp protocol $\eqref{g}$ is shown in figure 10.

From the knowledge of the single-particle wave functions $\eqref{one}$ we can write down the $N$-particles state as

[TABLE]

with $\vec{x}\equiv(x_{0},x_{1},...,x_{N-1})$. The Vandermonde determinant $\Delta$ symmetrizes the Slater determinant under particle exchange giving us the wave function of $N$ hard-core bosons since

[TABLE]

Introducing the generating functional

[TABLE]

the time-dependent particle density can be expressed as a functional derivative of $\mathcal{Z}[a]$:

[TABLE]

Using the random matrix approach [43, 44], we end up with the explicit result [33, 39]

[TABLE]

where $\ell_{0}\equiv\ell(t_{0})$. The effect of the dynamics is completely absorbed in the definition of a non-trivial length scale

[TABLE]

The adiabatic limit is recovered by $\lim_{t\rightarrow t_{0}}\zeta(t)=1$, for which $\xi=\ell_{0}$.

5.2 Comparison with the quasi-adiabatic case

Let us consider the protocol (54) in the quasi-adiabatic regime, i.e. when the quench rate $\alpha\rightarrow 0$. Expanding $\zeta(t)$ and $\omega^{2}(t)$ to the leading order in $\alpha$

[TABLE]

we can solve (53) perturbatively. At the zeroth order, (53) leads to

[TABLE]

which has the trivial solution $\zeta_{0}=1$, which is necessary for continuity at $t_{0}$. The first order equation is

[TABLE]

with initial condition $\zeta_{1}(0)=0$, $\overset{\centerdot}{\zeta}_{1}(0)=0$ and its solution is

[TABLE]

In the scaling limit we considered so far ($g_{0}\sim 1/L^{2}$) the function $\zeta(t)$, to the leading order in $\alpha$, shows a cubic growth in time:

[TABLE]

The quasi-adiabatic density profile (59) can be characterized through the length scale

[TABLE]

instead of the adiabatic length scale

[TABLE]

identified in the previous section (see (24)). In the perturbative regime considered here, $g_{0}t^{3}\ll t$ and the expansion/contraction of the cloud is always slower than what would have been expected from a naive adiabatic guess. This feature is related to the freezing out of the dynamics close to the critical point with the consequent breakdown of the adiabatic behavior. In figure 11 we show a comparison between the density profiles obtained from this Ermakov Lewis approach and the quasi-adiabatic approach developed in the previous section. The comparison is made using $\varepsilon=\nu_{g}\,N\,\alpha\,t$, which follows from (35) and (54), and we see an excellent agreement between both aproaches.

5.3 The special case of the decrease of the frequency as the inverse of time

A special case of interest is the situation where the harmonic trap frequency decreases as the inverse of time

[TABLE]

from an initial time $t_{0}=1$ set to one in the following such that the initial frequency is $\omega_{0}\equiv\omega(t_{0})=\lambda$. This time dependence is generated by the ramp

[TABLE]

In a recent work [45] it has been shown that the release of a scale-invariant Fermi gas confined within a harmonic trap with this type of $1/t$ time dependence leads to the appearance of a discrete scaling symmetry in time. Such a discrete scale invariance is known to produce log-periodic modulations of the physical quantities [46]. This behaviour has been observed and reported in [45] where the size of the expanding Fermi gas grows through a sequence of plateaus which are distributed log-periodically. As the free Fermi gas is closely related to the Tonks Girardeau gas one expects that such a phenomenon also exist in that case [45]. Indeed, here we prove that there is a regime where this log-periodic modulation of the expansion appears. Considering the Pinney equation with (68)

[TABLE]

one can derive an equivalent Pinney equation with a time independent frequency. Indeed, with the substitution

[TABLE]

we arrive at

[TABLE]

with boundary conditions

[TABLE]

Therefore, we have two distincts regimes, one with a high initial frequency $\omega_{0}=\lambda>1/2$ for which

[TABLE]

and an other one at low initial frequency, $\omega_{0}=\lambda<1/2$, for which

[TABLE]

The solution of (72) with the boundary conditions (73) is in the high initial frequency regime

[TABLE]

and

[TABLE]

in the low initial frequency one. The typical size of the bosonic cloud $2\ell_{N}(t)=2\sqrt{2N}\xi(t)$ is extracted from the exact density profile (59) and it is given through

[TABLE]

Clearly, we see in the high initial frequency case, that is for $\lambda>1/2$, a square root expansion with a typical log-periodic modulation. On the contrary, at lower initial frequencies, for $\lambda<1/2$, the expansion of the cloud at long times is governed by a pure power-law

[TABLE]

with a dynamical exponent $z(\lambda)=1/x(\lambda)$ that varies continuously with the initial frequency $\lambda$ and which is given through

[TABLE]

Such a behaviour is reminiscent of marginal perturbations, such as the Hilhorst-van Leeuwen ones [47, 48, 49], affecting the equilibrium critical exponents at a second order phase transition continuously. Notice that the log-periodic modulation can be seen as the emergence of a complex exponent $x(\lambda)$. Indeed, in the high frequency regime, the parameter $\kappa(\lambda>1/2)=is(\lambda)$ and the power law behaviour $t^{1/2+\kappa(\lambda)\lambda}=t^{1/2(1+i2s(\lambda)\lambda)}$. Taking its real part gives the log-periodic modulation:

[TABLE]

In figure 12, we report the typical length-scale $\eqref{xiEfi}$ for $\lambda>1/2$ and the behavior of the dynamical exponent $\eqref{Xexp}$ as a function of $\lambda$.

6 Summary and conclusions.

We have investigated the unitary dynamics of a low-density one-dimensional gas of impenetrable bosons at zero temperature under a harmonic trap release. In particular, we have studied the off-equilibrium aspects emerging from a slow general time ramp of the trapping potential. The presence of the trap leads to the emergence of a typical length scale $\ell(t)$ (see $\eqref{ell}$) which fully characterizes the time evolution of the cloud for an adiabatic process (see figure 3). We have computed the first off-equilibrium corrections arising away from adiabaticity using time-dependent perturbation theory. The departure from the equilibrium density profile has been obtained and shown to exhibit a scaling form (47) in the large particle number limit. We have also provided an exact solution using dynamical invariants in the case of a linear time ramp. In that case, we have identified the exact typical length scale $\xi(t)$ governing the process. This length scale is the product of the instantaneous initial typical length scale $\ell(t_{0})$ with the solution $\zeta(t)$ of the Pinney non-linear differential equation $\eqref{Pinney}$. The connection with the first order perturbative result is made through a series expansion of the exact scale $\xi(t)$ in the limit of a small linear quench rate. Thanks to that the density profile, obtained from first order perturbation theory, matches perfectly the exact one. When the trap frequency decreases as the inverse of time, we have shown by an exact solution of the Pinney equation that the cloud expansion is modulated by a log-periodic function reminiscent of a discret scale invariance for high initial frequencies. If the initial frequency is low enough the expansion follows a power law with an exponent that continuously depends on the initial frequency. The proper identification of a typical length scale for these slow out-of-equilibrium processes may be useful in experimental contexts as for example in cold atoms setup. Indeed, the characteristic scales that we pointed out may give a theoretical hint on the fluctuations of the particle density due to smooth modification of the optical cavities.

We would like to thank J. Unterberger and J. Dubail for useful discussions. We also gratefully acknowledge R. Qi for putting to our attention the nice work reported in [45] who gave us the hint to consider the $1/t$ case treated in the last part of this work.