Finite-temperature dynamics of a Tonks-Girardeau gas in a frequency-modulated harmonic trap

TL;DR

This paper investigates the out-of-equilibrium finite-temperature dynamics of a harmonically trapped Tonks-Girardeau gas under periodic trap frequency modulation, providing exact solutions and stability analysis through Mathieu's equation.

Contribution

It offers the first exact solutions for the density and momentum distributions of this system and demonstrates the equivalence of many-body and hydrodynamic stability diagrams via Mathieu's equation.

Findings

Exact solutions for density and momentum distributions.

Identification of parametric resonances in the system.

Equivalence of many-body and hydrodynamic stability diagrams.

Abstract

We study the out-of-equilibrium dynamics of a finite-temperature harmonically trapped Tonks-Girardeau gas induced by periodic modulation of the trap frequency. We give explicit exact solutions for the real-space density and momentum distributions of this interacting many-body system and characterize the stability diagram of the dynamics by mapping the many-body solution to the solution and stability diagram of Mathieu's differential equation. The mapping allows one to deduce the exact structure of parametric resonances in the parameter space characterized by the driving amplitude and frequency of the modulation. Furthermore, we analyze the same problem within the finite-temperature hydrodynamic approach and show that the respective solutions to the hydrodynamic equations can be mapped to the same Mathieu equation. Accordingly, the stability diagram and the structure of resonances following from the hydrodynamic approach is exactly the same as those obtained from the exact many-body solution.

Equations

$${\Psi }_{N,k}^B(x_1,...,x_N;t)=A(x_1,...,x_N){\Psi }_{N,k}^F(x_1,...,x_N;t),$$

$$i\hslash \frac{\partial {\varphi }_{k_i}(x,t)}{\partial t}=-\frac{{\hslash }^2}{2m}\frac{{\partial }^2{\varphi }_{k_i}(x,t)}{\partial x^2}+V(x,t){\varphi }_{k_i}(x,t),$$

$$i\hslash \frac{\partial {\Psi }_{N,k}^F}{\partial t}=\hat{H}{\Psi }_{N,k}^F,$$

$$\hat{H}=\sum _{i=1}^N\left[-\frac{{\hslash }^2}{2m}\frac{{\partial }^2}{\partial x_i^2}+V(x_i,t)\right].$$

$$\rho (x,t)=\rho (x,x;t),$$

$$n(k,t)=\iint dxdy{e}^{-ik(x-y)}\rho (x,y;t).$$

$$\rho (x,y;t)=\sum _{i,j=0}^{\infty }\sqrt{f_i}{\varphi }_i(x,t)Q_{ij}(x,y;t)\sqrt{f_j}{\varphi }_j^{\ast }(y,t),$$

$$P_{ij}(x,y;t)={\delta }_{ij}-2\text{mbox}sgn(y-x)\sqrt{f_i{f}_j}{\int }_x^{y}d{x}^{\prime }{\varphi }_i(x^{\prime },t){\varphi }_j^{\ast }(x^{\prime },t),$$

$$V(x,t)=m{\omega }^2(t)x^2/2,$$

$${\varphi }_j(x,t)=\frac{1}{\sqrt{\lambda }}{\varphi }_j\left(\frac{x}{\lambda },0\right)\exp \left[i\frac{m{x}^2}{2\hslash }\frac{\dot{\lambda }}{\lambda }-i{E}_j(t)t\right],$$

$$\ddot{\lambda }+{\omega }^2(t)\lambda =\frac{{\omega }^2(0)}{{\lambda }^3},$$

$$\rho (x,y;t)=\frac{1}{\lambda }{\rho }_0\left(x/\lambda ,y/\lambda \right)e^{im\dot{\lambda }(x^2-y^2)/2\hslash \lambda },$$

$$\ddot{\lambda }+{\omega }^2(t)\lambda =0,$$

$$\lambda (t)=\sqrt{\mathcal{A}{\lambda }_1^2(t)+\mathcal{B}{\lambda }_2^2(t)+2\mathcal{C}{\lambda }_1(t){\lambda }_2(t)}$$

$${\partial }_t\rho +{\partial }_x(\rho v)$$

$${\partial }_{t}v+v{\partial }_{x}v=-\frac{1}{m}{\partial }_{x}V(x,t)$$

$${\partial }_{t}s+v{\partial }_{x}s$$

$$\rho (x,t)$$

$$v(x,t)$$

$$T(t)$$

$$\omega^{2}(t)=\omega_{0}^{2}\big{(}1-\alpha\sin(\Omega t)\big{)}.$$

$$\ddot{z}(\tau)+\big{(}a-2q\cos(2\tau)\big{)}\,z(\tau)=0,$$

$$a={\left(\frac{2{\omega }_0}{\Omega }\right)}^2\text{ and }q=\frac{2{\omega }_0^2\alpha }{{\Omega }^2}.$$

$${\lambda }_1(t)=C\left(\frac{4{\omega }_0^2}{{\Omega }^2},\frac{2\alpha {\omega }_0^2}{{\Omega }^2},\frac{\pi }{4}-\frac{\Omega t}{2}\right),$$

$${\lambda }_2(t)=S\left(\frac{4{\omega }_0^2}{{\Omega }^2},\frac{2\alpha {\omega }_0^2}{{\Omega }^2},\frac{\pi }{4}-\frac{\Omega t}{2}\right).$$

$$\lambda (t)=\frac{1}{|W|}$$

$$+{{\left({\dot{\lambda }}_2(0){\lambda }_1(t)-{\dot{\lambda }}_1(0){\lambda }_2(t)\right)}^2]}^{1/2}.$$

$${\omega }^2(t)={\omega }_0^2(1+\alpha \cos (\Omega t)).$$

$$z_1(\tau )=e^{i\nu \tau }f(\tau ),$$

$$f(\tau )=\sum _{n=-\infty }^{\infty }{c}_n{e}^{2in\tau },$$

$$\nu =\frac{2}{\pi }\arcsin \left(\sqrt{\Delta (0){\sin }^2(\sqrt{a}\pi /2)}\right),$$

$$\nu=\frac{1}{\pi}\arccos\big{(}2\Delta(1)-1\big{)}.$$

$$\Delta (\nu )\equiv \det \left({\delta }_{nm}+\frac{({\delta }_{n+1,m}+{\delta }_{n-1,m})q}{(2n+\nu {)}^2-a}\right),$$

Full text

I Introduction

Characterising and understanding the behaviour of quantum many-body systems driven out-of-equilibrium is one of the grand challenges of modern physics. The problem is particularly challenging when treating the dynamics of realistic finite-temperature systems, rather than systems evolving from the zero-temperature ground state. Here we solve this problem for a paradigmatic model in many-body physics—a Tonks-Girardeau (TG) gas Girardeau (1960) of bosons interacting in one dimension (1D) via hard-core repulsion. The specific dynamical protocol that we consider is periodic modulation of the frequency of the harmonic trapping potential, which continuously drives the system out-of-equilibrium. Periodically driven systems are rather common in nature and form an important class of problems in quantum dynamics. In the present context, a periodically driven TG gas can be realised in ultra-cold atom experiments Paredes et al. (2004); Kinoshita et al. (2004, 2006); Haller et al. (2009); Meinert et al. (2015), and therefore our results are of direct relevance to the ongoing experimental progress in ultra-cold atom physics, aimed at developing a fundamental understanding of out-of-equilibrium phenomena in many-body physics Bloch et al. (2008); Cazalilla and Rigol (2010); Polkovnikov et al. (2011); Cazalilla et al. (2011); lam .

Previously, a harmonically trapped single particle has been widely studied and the effect of modulating either the position of the trap center or the strength of the trapping potential has been thoroughly investigated Husimi (1953); Breuer and Holthaus (1989); He et al. (2014); Brown (1991); Hänggi and Zerbe (1993); Kohler et al. (1997); Zerbe et al. (1994); Thorwart et al. (2000). More recently, in the context of many-body out-of-equilibrium dynamics, the periodic modulation of the trapping frequency of a one-dimensional system of interacting bosons at zero temperature was considered by Quinn and Haque Quinn and Haque (2014). In particular, they focussed primarily on the energy absorption of the system and the structure of energy resonances for interaction strengths varying from no interactions through to the strongly interacting TG regime.

In this work, we extend the analysis of Quinn and Haque in the TG regime to characterise the non-equilibrium dynamics of the gas at finite temperatures and show how the solution to this many-body problem can be mapped to the solution of Mathieu’s differential equation Mathieu (1868); McLachlan (1964). The approach we use has been developed recently in Ref. Atas et al. (2017a) and provides an exact finite-temperature theory that can in general model the dynamics of a TG gas in arbitrary trapping potentials. So far this treatment has been applied to a harmonically trapped TG gas driven out of equilibrium via a quench of the trap frequency from an initial value to a fixed final value Atas et al. (2017b). Such a quench protocol invokes familiar breathing-mode oscillations of the density profile, however, an unexpected many-body effect—dubbed many-body bounce—was also found in the dynamics of the momentum distribution. The many-body bounce manifests itself as an additional narrowing of the momentum distribution during the breathing mode oscillation cycle and is similar to the phenomenon of frequency doubling first observed in a weakly interacting 1D Bose gas Fang et al. (2014); Bouchoule et al. (2016).

The mapping to Mathieu’s equation is facilitated through a scaling solution that exist for the nonequilibrium dynamics of the TG gas in a time-varying harmonic trap. Minguzzi and Gangardt (2005); Quinn and Haque (2014); Atas et al. (2017a, b). Such a scaling solution allows the evolution of the density and momentum distributions of the gas to be determined through a single scaling parameter, which itself satisfies an ordinary second-order nonlinear differential equation, the Ermakov-Pinney equation Erm ; Pinney (1950). For a specific case of sinusoidal modulation of the trap frequency, the solutions to the Ermakov-Pinney equation can be found through the mapping to Mathieu’s equation. Mathieu’s equation itself possesses stable and unstable solutions, with a nontrivial structure of parametric resonances, which ultimately determine the dynamics and the stability diagram of the many-body TG gas that we analyse here in detail.

In addition to the exact many-body analysis, we study the dynamics of the TG gas using a finite-temperature hydrodynamic approach developed in Ref. Bouchoule et al. (2016). Solutions to the hydrodynamic equations are facilitated through the same scaling parameter (satisfying the same Ermakov-Pinney equation) as that in the exact many-body theory Erm ; Pinney (1950). Accordingly, the mapping to Mathieu’s equation, the stability diagram and the structure of parametric resonances that follow from the hydrodynamic approach are also the same. Thus, the frequency-modulation of a harmonically trapped Tonks-Girardeau gas represents yet another dynamical scenario, in addition to a simple quench protocol studied in Ref. Atas et al. (2017b), which can be adequately described by such a theory.

II The model and scaling solutions

II.1 Exact many-body solution in a time-dependent trapping potential

The model system we consider is a 1D gas of $N$ impenetrable (hard-core) bosons Girardeau (1960); Lieb and Liniger (1963) of mass $m$ confined by a one-body time-dependent trapping potential $V(x,t)$. (To keep the discussion general, we do not specify the form of $V(x,t)$ until the specifics of harmonic trapping come into play.) Due to the Fermi-Bose gas mapping Girardeau (1960); Girardeau and Wright (2000); Das et al. (2002), to every $N$-particle wavefunction of the bosonic problem, $\Psi^{B}_{N,k}(x_{1},...,x_{N};t)$ (where the index $k$ enumerates the different states that may occur in the given $N$-particle sector), there corresponds an antisymmetric wavefunction $\Psi^{F}_{N,k}(x_{1},...,x_{N};t)$ for a similarly trapped system of noninteracting spinless fermions, and vice versa:

[TABLE]

where the unit antisymmetric function $A(x_{1},...,x_{N})\!=\!\prod_{1\leq j<i\leq N}\mbox{sgn}(x_{i}-x_{j})$ ensures the symmetrization of the bosonic wavefunction. Accordingly, observables that do not depend on the sign of the many-body wave functions are readily given by their fermionic counterparts Das et al. (2002), whereas those that depend on the sign behave differently.

The free-fermion $N$-particle wavefunction itself is constructed as a Slater determinant $\Psi_{N,k}^{F}(x_{1},...,x_{N};t)\!=\!\mathrm{det}_{i,j=1}^{N}\left[\phi_{k_{i}}(x_{j},t)\right]/{\sqrt{N!}}$ of single-particle wavefunctions $\phi_{k_{i}}(x,t)$ evolving according to the time-dependent Schrödinger equation,

[TABLE]

with the initial wavefunctions $\phi_{k_{i}}(x,0)$ being the eigenstates of the trapping potential $V(x,0)$, with eigenenergies $E_{k_{i}}$ such that the total energy $E_{k}\!=\!\sum_{i=1}^{N}E_{k_{i}}$. Thus, the index $k\!=\!\{k_{1},...,k_{N}\}$ in $\Psi^{B/F}_{N,k}$ labels the different sets of single-particle quantum numbers $k_{i}$ that may occur in the initial $N$-particle state. The initial $N$-particle state is described by $\Psi_{N,k}^{F}(x_{1},...,x_{N};0)$, which itself evolves according to the $N$-particle Schrödinger equation,

[TABLE]

where

[TABLE]

We next assume preparation of the system such that its initial state at time $t\!=\!0$, in the trapping potential $V(x,0)$, is described by a grand-canonical ensemble at temperature $T_{0}$ and chemical potential $\mu_{0}$. This defines the statistical weights $P_{N,k}\!=\!\frac{1}{\mathcal{Z}}e^{(\mu_{0}N-E_{k})/k_{B}T_{0}}\geq 0$ attached to the occurrence of states $\Psi^{B}_{N,k}(x_{1},...,x_{N};0)$ in the ensemble, where $\mathcal{Z}=\!\sum_{N,k}e^{(\mu_{0}N-E_{k})/k_{B}T_{0}}$ is the grand-canonical partition function. With the initial state preparation specified, the Fermi-Bose mapping allows one to reduce the problem of treating finite-temperature dynamics of a trapped TG gas to a single-particle basis. In particular, combining such mapping with the Fredholm determinant approach of Ref. Atas et al. (2017a) leads to a computationally efficient way of calculating exactly the reduced one-body density matrix of the system, $\rho(x,y,t)$, and hence important physical observables of the gas, such as its real-space density

[TABLE]

and the momentum distribution

[TABLE]

More specifically, the reduced one-body density matrix can be computed as a simple double sum Atas et al. (2017a),

[TABLE]

where $Q_{ij}$ are the matrix elements of the operator $\mathbf{Q}(x,y;t)=(\mathbf{P}^{-1})^{\mathsf{T}}\mathrm{det}\;\!\mathbf{P}$, with

[TABLE]

and $f_{i}$ are the Fermi-Dirac occupancies $f_{i}\!=\![e^{(E_{i}-\mu)/k_{B}T_{0}}+1]^{-1}$ of single-particle orbitals $\phi_{i}(x,0)$.

II.2 Scaling solution in a harmonic trap

In general, the calculation of the one-body density matrix using Eq. (7) requires the evaluation of the overlap matrix elements $P_{ij}(x,y;t)$ between the time-evolved single-particle wavefunctions $\phi_{j}(x,t)$, starting from the initial wavefunctions $\phi_{j}(x,0)$. For the case of a harmonic trap,

[TABLE]

the initial wavefunctions $\phi_{j}(x,0)$ are given by the well-known Hermite-Gauss orbitals for the initial trap of frequency $\omega_{0}\!=\!\omega(0)$, with the 1D harmonic oscillator energy eigenvalues given by $E_{j}=\hbar\omega_{0}(j+\frac{1}{2})$, whereas the evolved wavefunctions, for arbitrary $\omega(t)$, can be found exactly using a scaling transformation Popov and Perelomov (1969); Perelomov and Zel’dovich (1998); Minguzzi and Gangardt (2005),

[TABLE]

where $E_{j}(t)\!=\!\hbar\omega_{0}(j+\frac{1}{2})\frac{1}{t}\int_{0}^{t}\mathrm{d}t^{\prime}/\lambda^{2}(t^{\prime})$, and the scaling function $\lambda(t)$ satisfies a second-order nonlinear differential equation known as the Ermakov-Pinney equation Erm ; Pinney (1950),

[TABLE]

with the initial conditions $\lambda(0)\!=\!1$ and $\dot{\lambda}(0)\!=\!0$.

The scaling transformation (10) leads to a simple scaling solution for the evolution of the reduced one-body density matrix Atas et al. (2017a, b),

[TABLE]

where $\rho_{0}(x,y)\!=\!\rho(x,y;0)$ is the one-body density matrix for the initial thermal equilibrium state. Such a scaling solution simplifies the analysis enormously as the one-body density matrix needs to be calculated only once, at time $t\!=\!0$, while the specifics of the dynamics is reduced to the solution of the Ermakov-Pinney equation for a single scaling parameter $\lambda(t)$

As has been shown by Pinney Pinney (1950), the general solution to the Ermakov-Pinney equation (11) can be constructed by combining two independent solutions \big{(}\lambda_{1}(t),\lambda_{2}(t)\big{)} of the respective homogeneous equation

[TABLE]

which describes a simple harmonic oscillator with time-dependent frequency. More specifically, the general solution is constructed as

[TABLE]

where the constants $\mathcal{A},~{}\mathcal{B}$, and $\mathcal{C}$ satisfy the constraint equation $\mathcal{A}\mathcal{B}-\mathcal{C}^{2}=\omega^{2}(0)/W^{2}$, and where $W=\lambda_{1}\dot{\lambda}_{2}-\lambda_{2}\dot{\lambda}_{1}$ is the Wronskian (which is a constant, in accordance with Abel’s identity) of the two independent solutions. This constraint equation, together with the initial conditions on the scaling function, $\lambda(0)\!=\!1$ and $\dot{\lambda}(0)\!=\!0$, uniquely fix the coefficients $\mathcal{A}$, $\mathcal{B}$ and $\mathcal{C}$. Note that the function $\lambda$ must be real and positive in order to keep the probability densities $|\phi_{j}(x,t)|^{2}$ positive, and as such the positive branch of the square root is assumed in (14).

In Sec. III below we will consider a specific dynamical protocol in which the TG gas evolves under sinusoidal modulation of the harmonic trap frequency $\omega(t)$. In this case, the homogeneous differential equation (13) can be mapped to Mathieu’s equation. Solutions to Mathieu’s equation and their stability properties are known and will be presented in Sec. IV, hence allowing us to construct the solutions to the Ermakov-Pinney equation and ultimately analyse the dynamics of the density and momentum distributions of the TG gas. Before doing so, however, we momentarily pause to point out that evolution of a harmonically trapped TG gas under the same dynamical protocol can be also analysed using a finite-temperature hydrodynamic approach developed in Ref. Bouchoule et al. (2016); Atas et al. (2017b).

II.3 Hydrodynamic treatment

In the hydrodynamic approach, the TG gas evolves according to the following equations for the local 1D density, $\rho(x,t)$, the hydrodynamic velocity, $v(x,t)$, and the entropy per particle, $s(x,t)$ Bouchoule et al. (2016); Atas et al. (2017b):

[TABLE]

where $V(x,t)$ is the trapping potential (given by Eq. (9 in the present case), and $P(x,t)$ is the local pressure that follows from the thermodynamic equation of the state.

As the equation of state for the TG gas is the same as that for an ideal Fermi gas, the solutions to Eqs. (15) through (17) are also the same as those for the ideal Fermi gas; they are given by the following scaling transformations Bouchoule et al. (2016); Atas et al. (2017b),

[TABLE]

Here, $\rho(x,t=0)=\rho_{0}(x)$ is the initial density profile and $T_{0}$ is the initial temperature of the gas, whereas $\lambda(t)$ is a scaling parameter which satisfies the same Ermakov-Pinney equation (11) as in the exact many-body treatment Bouchoule et al. (2016); Atas et al. (2017b). Therefore, under sinusoidal modulation of $\omega(t)$ as in Eq. (21) below, one can use the mapping to Mathieu’s equation and the same solutions for the scaling parameter $\lambda(t)$ that follow from them. Accordingly, the stability properties of these solutions (see Sec. IV.2) have the same implications on the hydrodynamics of the TG gas as in the exact the many-body treatment, implying that the ensuing stability diagram and the structure of parametric resonances of the TG in a frequency-modulated harmonic trap can be reproduced exactly from the hydrodynamic approach.

Specific examples of calculations of the dynamics of the TG gas that follow from the hydrodynamic approach will be presented in Sec. VI. These calculations are carried out as in Refs. Bouchoule et al. (2016); Atas et al. (2017b) and include the dynamics of the momentum distribution of the gas. The latter does not trivially follow from the above solutions for the real-space density distribution. Despite this, the momentum distribution of the gas can be constructed from the solutions for the density distribution, owing to the local density approximation that is intrinsic to the hydrodynamic approach. For ease of reference we outline this construction in Appendix A.

III Mapping to Mathieu’s equation for sinusoidal modulation

From here on we focus on the dynamics of the harmonically trapped TG gas in response to a sinusoidal modulation of the trap frequency,

[TABLE]

Here, $\omega_{0}$ is the trap frequency in the preparation stage of the initial ($t\leq 0$) thermal equilibrium state, whereas $\Omega$ and $\alpha$ ($0\leq\alpha\leq 1$), which define the parameter space of the problem, characterise the frequency and amplitude of subsequent modulation of $\omega(t)$ after time $t>0$.

Using $\omega^{2}(t)$ given by Eq. (21), along with the change of variable $\Omega t=\pi/2-2\tau$, one can show that the homogeneous differential equation (13) takes the canonical form of Mathieu’s equation Mathieu (1868) for the function z(\tau)=\lambda\big{(}(\pi/2-2\tau)/\Omega\big{)},

[TABLE]

which is parametrized in terms of

[TABLE]

The two independent solutions of Mathieu’s equation are given by the even and odd Mathieu functions denoted, respectively, by $C(a,q,\tau)$ and $S(a,q,\tau)$, and therefore the two independent solutions to the homogeneous Ermakov-Pinney equation can be written as

[TABLE]

[TABLE]

By combining these two solutions and solving the constraint equations for the coefficients $\mathcal{A}$, $\mathcal{B}$ and $\mathcal{C}$ in Eq. (14), the general solution to the Ermakov-Pinney equation can be written as

[TABLE]

This function indeed satisfies the initial conditions $\lambda(0)\!=\!1$ and $\dot{\lambda}(0)\!=\!0$, and we note that since the Wronskian of two independent solutions of a differential equation cannot be equal to zero the function $\lambda(t)$ is well defined.

We mention here that physically a slightly different modulation scenario to Eq. (21) is to use cosine modulation of the form

[TABLE]

In this case, the modulation of frequency begins slowly and without discontinuity in the derivative (i.e., with a zero slope in the time derivative, whereas in the former case of sine modulation the modulation begins abruptly). Nevertheless, the case of cosine modulation can still be solved using the solutions we present in this work for sine modulation; however, one must be careful in making the appropriate variable transformations. Specifically, the solutions for cosine modulation can be obtained from the solutions (24) and (25) by a simple time translation $t\rightarrow t-\pi/2\Omega$ and setting $q\rightarrow-q$, where $q$ is defined in Eq. (23).

IV Solutions to Mathieu’s equation

In what follows, we discuss the explicit form of the Mathieu functions, the method of calculating them, and their stability properties in the $(a,q)$ parameter space, which ultimately determine the physical properties of the periodically modulated TG gas. We consider the frequency modulation given by Eq. (21) and restrict ourselves to the case of $q>0$, or equivalently $\alpha>0$ (from Eq. (23)), with the recognition that the transformation $\tau\rightarrow\tau+\pi/2$ simply changes the sign of $q$ in Eq. (22). Therefore the case $q<0$ ($\alpha<0$) can be treated by simply translating the time $t\rightarrow t+\pi/\Omega$ in the solutions for $q>0$,(i.e., take $t\rightarrow t+\pi/\Omega$ with $q>0$ in Eqs. (24) and (25) and propagate these solutions through Eq. (26) to obtain $\lambda(t)$ for the equivalent $q<0$ case).

IV.1 Construction of the two independent solutions to Mathieu’s equation

The differential equation (22) represents the motion of a classical particle in a $\pi$-periodic potential, and therefore from Floquet’s theorem we have that a first fundamental solution can be written in the form

[TABLE]

where $\nu$ is the Floquet exponent (also referred to as the Mathieu characteristic exponent), which depends on the parameters $a$ and $q$, and $f(\tau)$ is a periodic function with period $\pi$, i.e., $f(\tau+\pi)=f(\tau)$. By virtue of its $\pi$-periodicity, the function $f$ can be represented using the Fourier expansion

[TABLE]

where $c_{n}$ are the expansion coefficients.

Inserting this form of $f(\tau)$ into Mathieu’s equation (22) one obtains a system of linear equations for the coefficients $c_{n}$, which can be written in matrix form. The condition of having a nontrivial solution to this set of equations results in a matrix determinant equation involving the Floquet exponent $\nu$ (see Appendix B for details). For given values of $a$ and $q$ this equation can be solved to obtain $\nu$. In particular, for $a\neq 4j^{2}$ (where $j$ is an integer) one finds that

[TABLE]

whereas for $a=4j^{2}$ ($j\in\mathbb{Z}$) one finds that

[TABLE]

Here, $\Delta(0)$ and $\Delta(1)$ are numerical coefficients to be obtained from

[TABLE]

at $\nu=0$ and $\nu=1$, respectively. In the above equation, the integer indices $n$ and $m$ in the infinite determinant run from $-\infty$ to $\infty$, but in practice they must be truncated at a high absolute value that is sufficient for numerical convergence.

From the properties of the matrix determinant involving $\nu$, one can further show (see Appendix B) that if $\nu$ is a solution to the determinant equation, then $\pm\nu\pm 2n$ (with $n$ being an integer) is also a solution. Therefore, from Eqs. (30) and (31) one can deduce that when the arguments of the inverse trigonometric functions lie within the interval $[-1,1]$ (depending on the values of $a$ and $q$) it is sufficient to only consider solutions for $\nu$ in the interval $\nu\in[0,1]$. The cases with $\nu=0$ and $\nu=1$ correspond, respectively, to even and odd integer values of the Floquet exponent. When, on the other hand, the arguments of the inverse trigonometric functions in (30) or (31) lie outside the interval $[-1,1]$, then the Floquet exponent becomes complex, in which case it can be either a pure imaginary number $\nu=i\eta$ or have the form $\nu=1+i\eta$, with $\eta<0$.

Once the Floquet exponent $\nu$ is known from Eq. (30) or (31), one can solve the system of equations for the Fourier coefficients $c_{n}$ (see Appendix C) and finally construct the first fundamental solution $z_{1}(\tau)$, Eq. (28), in explicit form. For constructing the second fundamental solution, $z_{2}(\tau)$, we note that Mathieu’s equation is unchanged by the time reversal transformation $\tau\!\rightarrow\!-\tau$. This means that

[TABLE]

is also an independent solution of Mathieu’s equation (22), provided that the Floquet exponent is either complex or real but not an integer. Therefore, the function $z_{2}(\tau)$ given by Eq. (33) can be taken as the second independent solution of Mathieu’s equation, in these cases.

When $\nu$ is an integer (we recall that the only two integer values of $\nu$ that we need to consider are $\nu=0,1$), the first fundamental solution $z_{1}(\tau)$ becomes $\pi$ or $2\pi$ periodic (according to the parity of $\nu$), and by Ince’s theorem McLachlan (1964) the second linearly independent solution cannot have either a period $\pi$ or $2\pi$. This means that $z_{2}(\tau)$, defined by Eq. (33), no longer constitutes an independent solution of Mathieu’s equation (22) as it becomes simply proportional to $z_{1}(\tau)$. Instead, the second independent solution in this case is given by $z_{2}(\tau)=\beta\tau z_{1}(\tau)+g(\tau)$ McLachlan (1964); Arscott (2014), where $\beta$ is a constant and the function $g(\tau)$ has the same periodicity as $z_{1}(\tau)$. An explicit expression for $z_{2}(\tau)$ in terms of a series of products of Bessel functions is given in Appendix C.

In conclusion, except for the special case when $\nu=0$ or $1$, the two independent solutions $z_{1}(\tau)$ and $z_{2}(\tau)$ can always be combined to form the even and odd solutions of Mathieu’s equation,

[TABLE]

such that $C(a,q,0)=1$ and $\dot{C}(a,q,0)=0$, whereas $S(a,q,0)=0$ and $\dot{S}(a,q,0)=1$. Note that the normalizations used here is quite arbitrary and has no consequences in the construction of the scaling function $\lambda(t)$ as the normalization constants are absorbed by the coefficients $\mathcal{A}$, $\mathcal{B}$ and $\mathcal{C}$ in Eq. (14).

IV.2 Stability of solutions to Mathieu’s equation

The long-time asymptotic behaviour of the scaling function $\lambda(t)$, and hence the behaviour of physical observables such as the density and momentum distributions of the trapped TG gas, are in direct bijection with the stability properties of the solutions to Mathieu’s equation. As such, to characterise the long-time dynamical behaviour of the harmonically trapped TG gas, one need only to examine the stability properties of these solutions. The solutions are said to be stable if they remain bounded when $\tau\rightarrow\infty$, and they are unstable if they tend to $\pm\infty$ when $\tau\rightarrow\infty$. From Eqs. (28) and (33), we can see see that stable solutions correspond to $\nu$ being real, whereas unstable ones correspond to $\nu$ being complex. In addition to their stability, the solutions can be classified according to their periodicity: when $\nu=p/r$ is a rational fraction less than unity, with $p$ and $r$ being mutually prime integers, the solutions have a period of $2\pi r$, and when $\nu$ is an irrational number the solutions have no specific periodicity.

In Figs. 1 (a) and (b) we show the density plots of, respectively, the real and imaginary parts of the Floquet exponent $\nu$ in the $(a,q)$ plane, calculated using Eqs. (30) and (31). On these plots, the dashed (magenta) lines show the boundaries between the stable and unstable regions. The Floquet (or the Mathieu characteristic) exponent on these boundaries is an integer (modulo the periodicity), $\nu(a,q)=n\,\mathrm{mod}\,2$ (i.e., $\nu=0$ or $1$), with $n\in\mathbb{Z}$ indexing the band. For a given $q$, the values of the parameter $a$, for which the Floquet exponent $\nu(a,q)$ lies on the stability boundaries are called the characteristic values of Mathieu’s equation and are denoted as $a_{n}(q)$ for even (cosine elliptic) solutions and $b_{n}(q)$ for odd (sine elliptic) solutions. As can be seen from the Floquet solution (28) and the Fourier expansion (29), on the even band the fundamental solution $z_{1}(\tau)$ is $\pi$ periodic, whereas on the odd band the fundamental solution has period $2\pi$. We recall that, for $\nu=0$ or $1$, only the fundamental solution $z_{1}(\tau)$ has a period of $\pi$ or $2\pi$, whereas the second independent solution is not periodic and is given by a series involving products of Bessel functions (see Appendix C).

The parameter space between the different bands corresponding to unstable solutions is shown in Fig. 1 (a) by black (where $\mathrm{Re}\,\nu=0$) and white (where $\mathrm{Re}\,\nu=1$) regions. On the other hand, the parameter space corresponding to stable solutions is shown by the colour-mapped regions, located in $a_{n}(q)<a(q)<b_{n+1}(q)$ (with $n=0,1,2,\dots$ and $a_{0}<0$). The Floquet exponent is real in these regions, with $0<\nu<1$, and the solutions to Mathieu’s equation are given by (28) and (33). The Fourier expansion coefficients in Eq. (29) are also real in these regions, and lead to real-valued solutions.

Respectively, in Fig. 1 (b), the black regions, where $\mathrm{Im}\,\nu=0$, correspond to stable solutions, whereas the colour-mapped regions correspond to unstable solutions. In the unstable regions, the Floquet exponent is either pure imaginary $\nu=i\eta$, with $\eta>0$, which is the case in the regions $b_{2n}(q)<a(q)<a_{2n}(q)$ [where $\mathrm{Re}\,\nu=0$ in Fig. 1 (a)], or it is given by $\nu=1+i\eta$, with $\eta<0$, which is the case in the regions $b_{2n+1}(q)<a(q)<a_{2n+1}(q)$ [where $\mathrm{Re}\,\nu=1$ in Fig. 1 (a)]. In these unstable regions, the Fourier expansion coefficients are generally complex, which ultimately leads to complex solutions of Mathieu’s equation. In practice, however, it is simpler to work with real valued solutions, and in Appendix C we show how one can use the symmetry of the Fourier expansion coefficients in order to absorb their complex phase factor in the definition of the first and second fundamental solutions, after which these solutions become real valued.

Discarding the actual values of the Floquet exponent in the $(a,q)$ parameter space and concentrating merely on whether $\nu$ is real or complex leads to the familiar stability digram (see, e.g., McLachlan (1964)) of solutions to Mathieu’s equation, which is shown in Fig. 2 in the positive quadrant of the $(a,q)$ plane.

V Stability diagram of the TG gas

We now return to our physical problem of the TG gas in the frequency-modulated harmonic trap. The frequency modulation is parametrised in terms of two independent parameters, $\Omega$ and $\alpha$, characterising, respectively, the frequency and amplitude of modulation of the trap frequency through Eq. (21). Mathieu’s equation, on the other hand, is parametrised in terms of $a$ and $q$, and the mapping between these pairs of parameters is given by Eq. (23). This particular form of mapping means that Mathieu’s parameters $a$ and $q$ are not independent, but are related by

[TABLE]

Moreover, the physical restriction $0\leq\alpha\leq 1$ on the amplitude parameter $\alpha$ restricts the values of $q$ to $0\leq q\leq a/2$, hence the physically unattainable region of $q>a/2$ or $a<2q$, shown in Fig. 2.

For given values of the physical parameters $(\Omega,\alpha)$ (where $\Omega$ is measured in units of $\omega_{0}$), the point $(a,q)$ is uniquely specified, and gives the stability of the solution to Mathieu’s equation through the Floquet exponent $\nu(a,q)$. However, from an experimental point of view it is more practical to fix, for instance, the modulation amplitude $\alpha$ and search for modulation frequencies $\Omega$ (or values of $a$) that lead to the desired stable or unstable dynamics. It is therefore more convenient to work in the $(a,\alpha)$ parameters space, rather than in the $(a,q)$ space.

Using $(a,\alpha)$ as our parameters, we can now map the stability diagram of Fig. 2 to the parameter space of $(a,\alpha)$ by first fixing the amplitude $\alpha$, or equivalently the slope $2/\alpha$ of any straight line $a=(2/\alpha)q$ in Fig. 2, and then following the trace of this line through stable and unstable regions. The intersections of this line with consecutive characteristic lines $a_{n}(q)$ and $b_{n+1}(q)$, where $n=0,1,2,\dots$, define a set of values of $q$ and hence $a$ in between which the solution is stable. In this way, we can arrive at the new stability diagram in the $(a,\alpha)$ parameter space shown in Fig. 3, which is equivalent to the physically attainable part of the diagram of Fig. 2.

V.1 Parametric Resonances

The stability diagram of Fig. 3 can also be read in terms of an experimental scenario in which one fixes the modulation frequency $\Omega$ (or the parameter $a$) and then searches for the modulation amplitude $\alpha$ for stable or unstable dynamics. This corresponds to following a straight horizontal line at a fixed value of $a$ in Fig. 3; the values of $\alpha$ for which such a line remains in the white (grey) region give the range of modulation amplitudes for which the dynamics are unstable (stable). Examples of such lines, corresponding to driving frequencies $\Omega$ which are equal to fractions of $2\omega_{0}$,

[TABLE]

are shown in Fig. 3 (a) for $j=1,2,..,5$.

These special values of frequencies are the same as in the phenomenon of parametric resonance, well-known in the problem of a classical oscillator subject to (an externally driven) time-variation of the parameters Landau and Lifshitz (1976), as well as in the problem of a (single-particle) quantum oscillator with time-varying frequency Popov and Perelomov (1969). In these systems, as well as in the present many-body problem of the TG gas, the oscillatory response of the system can be resonantly amplified (and is hence unstable) at certain driving frequencies.

The primary parametric resonance corresponds to $a=1$, or to $j=1$ in $\Omega_{j}=2\omega_{0}/j$, which is the frequency of the natural breathing mode oscillations Atas et al. (2017b) of the TG gas. It occurs for arbitrary values of the modulation amplitude $\alpha$ as indicated by the fact that the respective horizontal line in Fig. 3 (a) remains in the white (unstable) region for all $0<\alpha\leq 1$. The second resonance, corresponding to $a=4$ (or $j=2$), displays a similar behaviour.

In contrast, for $j\geq 3$ the lines of fixed $a=j^{2}$ (or $a=9,16,25,...$) start (for $\alpha>0$) in the gray (stable) regions. This is because the lower boundaries of the unstable regions (dashed lines in Fig. 3) corresponding to $j\geq 3$ all have positive curvature at small $\alpha$ and remain above the respective values of $a=9,16,25,...$ for all $\alpha$. As such, the resonant enhancement (i.e., unstable behavior) occurs at values of $a$ that are shifted upwards relative to $a=9,16,25,...$ (in terms of frequencies, $\Omega$ is shifted downward relative to the exact fractional values of $\Omega_{j}=2\omega_{0}/j$). We note that this upward shift, as well as the width of unstable regions, increases with $j$; in fact, for relatively small modulation amplitudes, the width of the resonance region $\Delta a_{j}$ scales as $\Delta a_{j}\approx\frac{j^{2j}\alpha^{j}}{2^{3(j-1)}[(j-1)!]^{2}}$ Bell (1957). Moreover, for large enough $j$ (namely, $j\geq 8$) it is possible for the upper bound of the unstable region to cross the next resonant line of constant $a=64,81,...$, or in other words, it becomes possible for the horizontal line of fixed $a=64,81,...$ to re-enter the white (unstable) region of the previous (lower $j$) resonance as $\alpha$ is increased from 0 to 1.

The structure of resonances of the parametrically driven TG gas has been previously discussed in Ref. Quinn and Haque (2014), however, the exact boundaries of stable and unstable regions, as obtained here in Fig. 3 from the stability properties of Mathieu’s equation, have not been identified to the best of our knowledge.

VI Evolution of the density and momentum distributions

We first reiterate that the long-time dynamical behaviour of the trapped TG gas can be determined through the stability diagram presented in Fig. 3. For analysing the transient behaviour, however, one needs to resort to explicit time-dependent calculations of the dynamics. Such calculations provide insight on how the choice of modulation parameters, as well as the finite temperature of the gas, can affect dynamical features at short to intermediate times. In what follows, we thus examine both the real-space density and momentum distributions of the gas and present two typical examples which exemplify common behaviours in both the stable and unstable regions of the parameter space, using both the exact many-body approach and the hydrodynamic approach outlined in Sec. II. We also show how finite-temperature effects alter the behaviour of the gas in these two examples.

VI.1 Unstable dynamics

In Fig. 4 we show the evolution of the scaling function $\lambda(t)$ for modulation parameters $(a,\alpha)=(1,0.6)$. The choice $(a,\alpha)=(1,0.6)$ lies on the primary parametric resonance, corresponding to $j=1$ in Eq. (37) and hence the driving frequency of $\Omega=2\omega_{0}$, and therefore the dynamics are unstable. Accordingly, the sequence of the peaks in the scaling solution $\lambda(t)$ grows exponentially (on a long-time scale), whereas the oscillations that accompany the exponential growth are due to the natural breathing mode behaviour Atas et al. (2017b).

Figure 5 shows the resulting dynamics of a finite temperature TG gas containing $N=16$ particles at a relatively low initial temperature $T_{0}$. The temperature is parametrized in terms of the dimensionless parameter $\theta_{0}=k_{B}T_{0}/N\hbar\omega_{0}$, where $N\hbar\omega_{0}/k_{B}$ is the temperature of quantum degeneracy. We show the evolution of the particle number density $\rho(x,t)$, the momentum distribution $n(k,t)$, and the respective widths of the gas determined from both the exact quantum theory and the hydrodynamic approach. According to Eq. (12), the magnitude of the density $\rho(x,t)$ depends inversely on $\lambda(t)$. Thus, the sequential minima of $\rho(0,t)$, when the the density profiles are the broadest, diminish as $\lambda(t)$ grows, whereas the respective peak widths (when the density profile is the narrowest) grow in proportion to $\lambda(t)$. The behaviour of the momentum distribution $n(k,t)$ is reciprocal to that of the real-space density. However, the sequence of narrow and broad distributions is interrupted by an additional narrowing of the momentum distribution at time instances when the density distribution is also narrow. This is a manifestation of the phenomenon of quantum many-body bounce, studied for the breathing mode oscillations of the Tonk-Girardeau gas in Ref. Atas et al. (2017b), which in turn is similar to the phenomenon of frequency doubling in the weakly interacting 1D quasicondensate Fang et al. (2014); Bouchoule et al. (2016).

We note that, for the lower temperature dynamics which we present in Fig. 5, the hydrodynamic theory predicts very well the features obtained from the exact theory.

In Fig. 6, we present the dynamics resulting from the same modulation parameters as in Fig. 5, but for a gas at a higher initial temperature of $\theta_{0}=0.5$. As can be seen from the width of the momentum distribution, the higher temperature example displays thermal broadening. The thermal broadening, however, is apparent only for a short amount of time, and eventually it is overwhelmed by the large-scale dynamical broadening resulting from the unstable nature of the scaling parameter $\lambda(t)$. Accordingly, most of the dynamical features remain largely unaltered by the temperature of the gas at longer times. However, the thermal broadening does significantly affect the many-body bounce Atas et al. (2017b), which is far less pronounced compared to the low temperature case.

Within the hydrodynamic theory, which is derived here for temperatures lower than the temperature of quantum degeneracy (so that the density profile is well approximated by an inverted semicircle), the real space density does not depend on temperature. As such, the density profile cannot capture any alterations from thermal effects. On the other hand, the momentum distribution in the hydrodynamic theory is explicitly dependent on temperature through the width of the Lorentzian distribution (see Appendix A). In our higher temperature example, this is evident in the overall broadening of the momentum distribution, which is in agreement with the exact theory. The hydrodynamic theory also accurately captures the blurring effect of the increased temperature on the momentum width at the inner turning points where the many-body bounce occurs, even though the quantitative agreement with the exact theory is poorer in the immediate vicinity of these inner turning points Atas et al. (2017b).

VI.2 Stable dynamics

In Fig. 7 we show the evolution of the scaling function $\lambda(t)$, for modulation parameters $a=16$ (or $\Omega=\omega_{0}/2$) and $\alpha=0.8$. This parameter choice corresponds to stable or bounded dynamics, which generally can be periodic or aperiodic. More specifically, the Floquet exponent in this example is $\nu\approx 0.1991$, for which the solution is aperiodic. However, this value of $\nu$ is close to $\nu=1/5$, which results in strictly periodic dynamics (with a period of $10\pi$, following from the periodicity conditions given in Section IV.2); accordingly, the example under consideration is nearly periodic and hence we only show the behaviour of $\lambda(t)$ within $\Omega t\in[0,10\pi]$.

The high-frequency oscillations in $\lambda(t)$ within this period are a result of beating arising from the contribution of only a few most significant Fourier coefficients in the expansion of Eq. (29). As $\rho(x,t)=\rho(x,0)/\lambda(t)$, the same high-frequency oscillations show up in the dynamics of the density profile and its width shown in Fig. 8. In contrast, the dynamics of the momentum distribution and its width show oscillations that occur at approximately twice the frequency of oscillations of the density profile due to the phenomenon of quantum many-body bounce. The thermal broadening effect in the higher temperature example in Fig. 8 is qualitatively the same as in the example of unstable dynamics (Sec. VI.1) and is more clearly seen in the width of the momentum distribution of the gas. Similarly, the comparison between the hydrodynamic and exact theories is qualitatively the same as before, and we restrict ourselves to showing the hydrodynamic results only for the widths of the density and momentum distributions [dashed lines in Fig. 8 (b), (d), (f), and (h)].

VII Summary

In conclusion, we have characterised the out-of-equilibrium dynamics of a finite-temperature harmonically trapped TG gas in response to periodic modulation of the trap frequency. The analysis was performed using exact quantum many-body theory as well as a finite-temperature hydrodynamic approach. Due to the scaling transformations, which exist for the single-particle orbitals of the system as well as in the hydrodynamic theory, this otherwise complex many-body problem can be reduced to the solution of a single differential equation, Eq. (11). For periodic modulations of the form given by Eq. (21), we have shown that the solution to Eq. (11) can be mapped to the solutions of Mathieu’s differential equation.

We have provided a detailed analysis of how to construct the solutions to Mathieu’s equation and identified the physical and unphysical parameter spaces of these solutions when they are mapped to describe the TG many-body problem. Using the stability properties of Mathieu’s equation in terms of Floquet theory, we have constructed the stability diagram of the TG system as a function of the modulation frequency and amplitude. The stability diagrams following from the exact quantum theory and the hydrodynamic approach are identical. We have also identified the structure of parametric resonances within this stability diagram.

In order to explore the transient dynamics of the TG system, we have provided examples of the modulation-induced dynamics in both the unstable and stable regions of the parameter space. We find that the collective many-body bounce effect reported in Ref. Atas et al. (2017b) persists in both of these dynamical regimes. Furthermore, in each regime we include examples at both a low and high temperature so that we can comment on thermal effects within this system. We note that higher temperatures cause significant thermal broadening to the width of the momentum distribution and also cause blurring of the collective many-body bounce. In the unstable regime, however, the overall thermal broadening of the momentum distribution is eventually suppressed and drowned out by the exponentially increasing nature of the dynamics.

Our results open the path to analysing periodic modulation protocols in the context of nonequilibrium thermodynamics of interacting quantum many-body systems. As an example, they provide the tools for calculating quantum work distributions and designing periodic cycles suitable for quantum many-body heat engines.

In view of applicability of our results to realistic physical systems, we note that the TG model corresponds to the limit of a trapped 1D Bose gas, or the Lieb-Liniger model Lieb and Liniger (1963), with a two-body contact interaction potential, $g\delta(x-x^{\prime})$, in which the 1D coupling strength $g\rightarrow\infty$. The relevant dimensionless interaction parameter in the Lieb-Liniger model is $\gamma(x)=mg/\hbar^{2}\rho(x)$, where $\rho(x)$ is the local 1D density Kheruntsyan et al. (2005). In the TG limit, this local dimensionless interaction parameter must tend to infinity at all $x$. Physically realisable strongly interacting 1D Bose gases are characterised, however, by a very large ($\gamma(x)\gg 1$) but finite interaction strength Paredes et al. (2004); Kinoshita et al. (2004, 2006); Haller et al. (2009); Meinert et al. (2015). Accordingly, the results following from the strictly TG model serve only as a good approximation for such strongly interacting gases. As such, our dynamical results in the unstable regime of modulation (where the peaks of the density $\rho(x,t)$ grow exponentially with time) will eventually become physically invalid as the dimensionless local interaction strength $\gamma(x,t)$ can become of the order of, or much smaller than, unity, for which the TG model becomes inapplicable.

Acknowledgements.

The authors acknowledge stimulating discussions with D. M. Gangardt and support by the Australian Research Council Discovery Project grant DP170101423.

Appendix A Finite-temperature hydrodynamics

For specific dynamical calculations using the scaling solutions to the hydrodynamic equations, the initial density profile in Eq. (18) can be approximated by the Thomas-Fermi (TF) semicircle Atas et al. (2017b),

[TABLE]

which is valid when the thermal energy scale is much less than the kinetic energy scale, i.e., when $k_{B}T_{0}\ll\hbar^{2}\rho_{0}(0)^{2}/m$. Here, $R_{\text{TF}}=2N/\pi\rho_{0}(0)=\sqrt{2N}l_{ho}$ is the TF size of the cloud, with harmonic oscillator length $l_{ho}=\sqrt{\hbar/m\omega_{0}}$. Hence, Eq. (38) together with Eq. (18) provide the low-temperature solution for the density profile of the TG gas within the hydrodynamic theory.

The momentum distribution of the TG gas is somewhat more involved to obtain. One first employs the local density approximation (LDA) and defines $\bar{n}(k;\rho,T)$ to be the equilibrium momentum distribution of a uniform TG gas of density $\rho$ at temperature $T$, normalised to $\rho$ (i.e., $\int dk~{}\bar{n}(k;\rho,T)=\rho$). In this approach, the momentum distribution of each uniform slice of gas can be added together to obtain the full momentum distribution for the trapped TG gas,

[TABLE]

Next, one models the momentum distribution of a uniform TG gas using a Lorentzian $\bar{n}(k;\rho,T)=(2\rho l_{\phi}/\pi)/[1+(2l_{\phi}k)^{2}]$ Bouchoule et al. (2016); Atas et al. (2017b), where $l_{\phi}(x,t)=\hbar^{2}\rho(x,t)/mk_{B}T(t)$ is the phase coherence length of the gas. Although this model is only valid for small momenta and low temperatures ($|k|\ll 1/l_{\phi}\ll\rho$), it captures well the bulk of $\bar{n}$ and provides the dominant contribution to the bulk of $n(k,t)$. Using this Lorentzian approximation for $\bar{n}(k;\rho,T)$, the full momentum distribution of the TG gas is given by Bouchoule et al. (2016)

[TABLE]

Defining the initial phase coherence length in the trap centre to be $l_{\phi}^{(0)}=\hbar^{2}\rho_{0}(0)/mk_{B}T(0)$ and applying the transformation $u=x/\lambda R_{\text{TF}}$ allows the momentum distribution to be written in dimensionless form as Bouchoule et al. (2016)

[TABLE]

where $\tilde{k}=kl_{\phi}^{(0)}$ and $\tilde{T}_{0}=T_{0}/T_{d}$, with $T_{d}=\hbar^{2}\rho_{0}(0)^{2}/2mk_{B}$ being the initial temperature of quantum degeneracy of a uniform 1D gas at density $\rho_{0}(0)$. In addition, $\tilde{\lambda}(\tau)\equiv\lambda(\tau/\Omega)$, where $\tau=\Omega t$ is the dimensionless time so that $\dot{\tilde{\lambda}}=d\tilde{\lambda}/d\tau=\dot{\lambda}(t)/\Omega$ Atas et al. (2017b).

For the purpose of performing numerical calculations, one can use the relations $l_{\phi}^{(0)}=\hbar^{2}\rho_{0}(0)/mk_{B}T(0)=l_{ho}\sqrt{2/N}/\pi\theta_{0}$ and $\tilde{T}_{0}=\pi^{2}\theta_{0}$ (where $\theta_{0}=k_{B}T_{0}/N\hbar\omega_{0}$ is the initial dimensionless temperature of the gas) to rewrite Eq. (41) as

[TABLE]

Here $C=4\sqrt{2N}/\pi^{3}\theta_{0}$, $B=8/\pi^{2}\theta_{0}^{2}N$, and $\bar{k}=kl_{ho}$. Note that, in order to construct the hydrodynamic momentum distribution given by Eq. (42), only $N$, $\theta_{0}$ and the ratio $\Omega/\omega_{0}$ need to be specified.

Appendix B Exact expression of the Floquet exponent

In this Appendix, we show how to obtain the exact expression for the Floquet exponent of Mathieu’s equation, following the same line of derivation as in Ref. McLachlan (1964). First, by inserting the Fourier expanded form of the periodic function $f(\tau)$, Eq. (29), into Mathieu’s equation (22), one can find that the expansion coefficients $c_{n}$ satisfy the following three-term recurrence relation:

[TABLE]

This system of linear equations can be written in a matrix form,

[TABLE]

where the elements of the tridiagonal matrix $M_{\nu}$ are given by

[TABLE]

and the vector $\mathbf{c}=(\dots,c_{-n},\dots,c_{-1},c_{0},c_{1},\dots,c_{n},\dots)$ contains the Fourier coefficients.

The system of equations (44) has a non-trivial solution if

[TABLE]

This condition actually constitutes an equation for the Floquet exponent $\nu$, and requires computation of the infinite determinant. In practice, however, the infinite determinant must be truncated, so in order to ensure the convergence of the determinant, let us divide (43) by $(2n+\nu)^{2}-a$ and define

[TABLE]

The determinant which appears in Eq. (46) then becomes

[TABLE]

and setting $\Delta(\nu)=0$ will provide non-trivial solutions for the Floquet exponent $\nu$.

Even in this form, the equation for $\nu$ remains rather intricate and difficult to solve. However, a closed form expression for the Floquet exponent can be obtained by making the following observations. First, since the value of a determinant does not change when rows are interchanged, it follows that the function $\Delta(\nu)$ is even [i.e. $\Delta(\nu)=\Delta(-\nu)$]. Next, the transformation $\nu\rightarrow\nu\pm 2$ leaves the determinant unchanged since rows of the matrix are just shifted by $\pm 1$ in that case. Hence, the function $\Delta(\nu)$ is periodic in $\nu$ with period $2$. The only singularities of $\Delta$ are simple poles located at $\nu_{n}=\pm\sqrt{a}-2n$, with $n\in\mathbb{Z}$. The function

[TABLE]

presents all the above mentioned properties, and by Liouville’s theorem there exists a constant $\beta$ such that $\Delta(\nu)-\beta\xi(\nu)$ is also a constant McLachlan (1964); Arscott (2014). A closed form expression for $\nu$ can now be obtained by determining this constant, as well as $\beta$ by extension. These two constants will provide an equation for $\Delta(\nu)$ which can be set equal to zero and rearranged to obtain an expression for $\nu$.

Since $\Delta(\nu)-\beta\xi(\nu)$ is a constant, it can be evaluated at any $\nu$. Taking the limit $i\nu\rightarrow\infty$, we see that $\Delta(i\nu\rightarrow\infty)=1$, as all the off-diagonal elements $\zeta_{n}$ of the matrix in (48) tend to zero in this limit, whereas the diagonal elements are equal to $1$, and $\xi(i\nu\rightarrow\infty)=0$. Consequently,

[TABLE]

and hence the constant $\beta$ is given by $\beta=(\Delta(\nu)-1)/\xi(\nu)$.

In particular, $\beta$ can be evaluated at $\nu=0$, provided that $a\neq 4j^{2}$ ($j\in\mathbb{Z}$) as to avoid the matrix elements $\zeta_{n}(0)$ becoming infinite for $n=j$, which results in $\beta=[1-\cos(\sqrt{a}\pi)][\Delta(0)-1]$. Substituting this into Eq. (50) gives an explicit expression for the determinant,

[TABLE]

In this expression, $\Delta(0)$ is the determinant (48) with the off-diagonal elements given by $\zeta_{n}(0)=q/(4n^{2}-a)$. Setting $\Delta(\nu)=0$ (to obtain a non-trivial solution for the Floquet exponent) and rearranging for $\nu$ results in the following closed expression for the Floquet exponent:

[TABLE]

Similarly, by evaluating $\beta$ at $\nu=1$, it can be shown that, when $a=4j^{2}$,

[TABLE]

Here, $\Delta(1)$ is the determinant (48) with the off-diagonal elements given by $\zeta_{n}(1)=q/((2n+1)^{2}-a)$. In practice, one fixes the point in the parameter space $(a,q)$ and then computes the Floquet exponent $\nu$ using Eq. (52) or (53), which are the same equations as Eqs. (30) and (31) of the main text.

Appendix C Symmetry of the Fourier coefficients and form of the solutions in the different regions

Once the Floquet exponent is obtained, the Fourier coefficients can be computed from the eigenvector of the matrix $(M_{\nu}-a)$ which corresponds to the zero eigenvalue, where $M_{\nu}$ is given by Eq. (45). From the exact expression of the Floquet exponent, (30) or (31), and the periodicity property of the determinant equation, Eq. (46), the value of $\nu$ can always be taken in the interval $[0,1]$ when it is real. In particular, the case $\nu=0,1$ corresponds to the integer values of the Floquet exponent. The Floquet exponent can also be complex and it is easy to see that in that case the Fourier coefficients (43) are also complex and lead to complex fundamental solutions (28) and (33). In practice, it is easier to work with real representations of the fundamental solutions. From the properties of the Floquet exponent outlined in the previous section, one can deduce symmetries for the Fourier coefficients and simplify the expression for the fundamental solutions in the different regions of the phase diagram (Fig. 1) in order to get a real representation. In what follows, we derive the exact expressions used in the numerical calculation of the two fundamental solutions in the different regimes of stability.

C.1 Case when $\nu$ is complex ($\nu=1-i\eta$)

When $(a,q)$ lies in between $b_{2n+1}(q)$ and $a_{2n+1}(q)$, the Floquet exponent has the form $\nu=1-i\eta$, which ultimately results in complex Fourier coefficients $c_{n}$. Consequently, the fundamental solution (28) is not necessarily a real function. We can however use the properties of the Fourier coefficients to get a real solution. For that purpose, let us insert the form of the Floquet exponent $\nu=1-i\eta$ into the recurrence relation (43). Setting $n\rightarrow-(n+1)$ and taking the complex conjugate of the resulting recurrence equation shows that the coefficients $c_{-n-1}^{\ast}$ satisfy the same recurrence relation as the coefficients $c_{n}$. This means that they must be proportional to each other, i.e., $c_{-n-1}^{\ast}=(c_{-1}^{\ast}/c_{0})\,c_{n}$, where the proportionality constant is found by taking $n=0$, and that they must have the same complex amplitude, i.e., $|c_{-n-1}^{\ast}|^{2}=|c_{n}|^{2}$. Hence, the proportionality constant is a pure phase factor $c_{-1}^{\ast}/c_{0}=\exp(-2i\theta)$, and we finally have

[TABLE]

Using this result in the Fourier expansion (29), we thus find that the fundamental solution (28) takes the form

[TABLE]

where we have used the polar representation of the Fourier coefficients $c_{n}=\rho_{n}\exp(i\phi_{n})$ and where $2\theta=\phi_{0}+\phi_{-1}$. The solution is complex due to the presence of the phase factor $e^{i\theta}$. For real $\tau$ however, it is possible to absorb the phase factor into the definition of the fundamental solution $z_{1}\rightarrow z_{1}e^{-i\theta}$ (which remains a solution) in order to get a real solution

[TABLE]

The second fundamental solution in this region of the phase diagram is given by

[TABLE]

These two solutions are then combined according to (34) and (35) to obtain the even and odd solutions, respectively, of Mathieu’s equation.

C.2 Case when $\nu$ is pure imaginary ($\nu=i\eta$)

When $(a,q)$ lies between $b_{2n}(q)$ and $a_{2n}(q)$, the Floquet exponent is a pure imaginary number $\nu=i\eta$ with $\eta>0$ and therefore the Fourier coefficients are complex in this case as well. By inserting the Floquet exponent into the recurrence relation (43) and setting $n\rightarrow-n$, along with the complex conjugate transformation, we deduce that the coefficients $c_{-n}^{\ast}$ are proportional to $c_{n}$. More precisely, we obtain $c_{n}=(c_{0}/c_{0}^{\ast})c_{-n}^{\ast}=\exp(2i\theta)c_{-n}^{\ast}$ with $\theta=\phi_{0}$. This property of the Fourier coefficients leads to the following form for the first fundamental solution:

[TABLE]

where, as before, we have absorbed the extra phase coefficient $e^{i\theta}$ into the definition of the solution.

The second solution is obtained by making the substitution $\tau\rightarrow-\tau$,

[TABLE]

C.3 Case when $\nu$ is real but not integer ($\nu\in\left]0,1\right[$)

When the pair $(a,q)$ lies in the stable region between $a_{n}(q)$ and $b_{n+1}(q)$, the Fourier coefficients are real and do not present any particular symmetry. The real valued even and odd solutions take the following form:

[TABLE]

with $0<\nu<1$. Notice that the second independent solution is obtained by replacing the cosine function by a sine function and that the normalization constant has been omitted in both representations.

C.4 Case when $\nu$ is integer ($\nu=0,1$)

As discussed above, $\nu=0,1$ correspond to real integer values of the Floquet exponent and in this case only one solution can be periodic, in accordance with Ince’s theorem McLachlan (1964). The second independent solution is constructed following a different scheme compared to the ones discussed previously. When $(a,q)$ lies exactly on a characteristic line (dashed magenta line in Fig. 1), i.e., $a=a_{n}(q)$ or $a=b_{n}(q)$, the Fourier coefficients are all real. The fundamental $\pi$ or $2\pi$ periodic solution is even on $a_{n}$ and odd on $b_{n}$.

C.4.1 On the characteristic lines of type $a_{n}(q)$

When $a=a_{2n}(q)$, the $\pi$ periodic fundamental solution has the form

[TABLE]

which is proportional to the function commonly denoted in the literature as $ce_{2n}(\tau,q)$ McLachlan (1964). The second fundamental non-periodic solution is given by

[TABLE]

where $J_{n}$ and $Y_{n}$ are the Bessel functions of first and second kind, respectively. This function is commonly denoted as $fe_{2n}(\tau,q)$ in the literature McLachlan (1964).

When $a=a_{2n+1}(q)$, the first fundamental solution has period $2\pi$ and can be written as

[TABLE]

whereas the second solution takes the form

[TABLE]

In the literature, one can find the notation $ce_{2n+1}(\tau,q)$ and $fe_{2n+1}(\tau,q)$ for these two solutions, respectively.

C.4.2 On the characteristic lines of type $b_{n}(q)$

When $a=b_{2n+1}$, the is odd and one has

[TABLE]

which is proportional to the commonly known function $se_{2n+1}(\tau,q)$ McLachlan (1964). The second solution has the form

[TABLE]

Finally, when $a=b_{2n+2}$ the solutions have the form

[TABLE]

and

[TABLE]