Nonequilibrium quantum thermodynamics of determinantal many-body systems: Application to the Tonks-Girardeau and ideal Fermi gases

The study of thermalization in isolated quantum systems has fuelled the recent development of quantum thermodynamics Gemmer et al. (2009); Binder et al. (2018); Vinjanampathy and Anders (2016); Kosloff (2013) and has highlighted its connection to the dynamics of quantum correlations, scrambling of quantum information, and entanglement dynamics Kim et al. (2011); Diaz de la Cruz and Martin-Delgado (2014); Nandkishore and Huse (2015); D’Alessio et al. (2016); Jaramillo et al. (2016). Furthermore quantum thermodynamics plays an essential role in designing quantum heat engines. Here the unique quantum properties of the working fluid of the engine can be exploited to achieve quantum supremacy in terms of the engine efficiency and output power. In this work we consider a system amenable to analytical and numerical exploration in the context of quantum thermodynamics.

Recent advances in quantum simulation platforms utilizing a variety of atomic, molecular, and optical platforms have facilitated the observation of fundamental topics in quantum thermodynamics in the presence of different particle statistics, flexible dimensionality, and short-range and long-range interactions Zheng and Poletti (2015); Gärttner et al. (2017); Schweigler et al. (2017); Friis et al. (2018); Brydges et al. (2019); Landsman et al. (2019). More recently, a variety of experimental platforms, including trapped-ions and nitrogen vacancy (NV) centers, have observed quantum effects in the absence of interactions von Lindenfels et al. (2019); Klatzow et al. (2019). This underscores the need for theoretical studies which consider the role of quantum many-body interactions. The next generation of experiments will feature quantum many-body interactions, allowing one to use quantum correlations and multi-particle entanglement as an additional resource to control the performance of quantum heat engines and quantum refrigerators. Trapped ultracold atomic gases lend themselves as a particularly promising platform in this respect Deffner and Lutz (2008); Abah and Lutz (2016); Sotiriadis et al. (2013); Chiara et al. (2015); Zheng and Poletti (2015); Jaramillo et al. (2016); Beau et al. (2016); Li et al. (2018); Rylands and Andrei (2019); Perfetto et al. (2019); Niedenzu et al. (2019); Roy and Eckardt (2020); Fogarty and Busch (2020) owing to a high degree of control over system parameters such as inter-atomic interactions and trapping potentials.

The quantum nature of the system complicates the calculation of physical quantities that are often considered in classical thermodynamics such as the mean work $\langle W(t)\rangle$. First, $W(t)$ is a randomly distributed quantity requiring two projective energy measurements at initial time and at time $t$ Vinjanampathy and Anders (2016); Kurchan (2000); Talkner et al. (2007); Yi et al. (2012); Campisi et al. (2011). Hence, quantum work is not represented by a Hermitian operator and hence cannot be regarded as an ordinary quantum observable Talkner et al. (2007) (see also Ref. Roncaglia et al. (2014); Cerisola et al. (2017)). Furthermore the outcome of these projective measurements depend on the transition probabilities between different quantum states. This probabilistic nature implies that in order to fully characterize and understand the work performed during time $t$ one needs to know the full work probability distribution $P(W)$ or its corresponding characteristic function $G_{\beta}(\vartheta)$. Finally, since these calculations involve many-body expectation values, their computational complexity may grow exponentially with system size. This is a major obstacle for a theoretical optimization of the engine performance as it limits any such study to small systems.

In this work we show that quantum many-body systems whose many-body wave functions can be represented as Slater determinants provide ideal platforms for exploring ideas of quantum thermodynamics. As we show here, we are able to recast the quantum thermodynamics of paradigmatic interacting many-body systems such as a 1D gas of hard-core bosons and anyons in terms of finite temperature dynamics, which itself can be reduced to the dynamics of single-particle wave functions. We illustrate this by using the Tonks-Girardeau (TG) gas Girardeau (1960, 1965) which is realizable experimentally Paredes et al. (2004); Kinoshita et al. (2004); Haller et al. (2009); Wilson et al. (2020), is amenable to analytical treatment for its finite-temperature dynamics Atas et al. (2017a), and as we have shown in our previous work Atas et al. (2019), the system dynamics in a periodically modulated harmonic trap Minguzzi and Gangardt (2005); Quinn and Haque (2014); Atas et al. (2017a, b); Scopa and Karevski (2017); Scopa et al. (2018) can be stable or unstable at a given modulation frequency depending on the amplitude of the modulation.

In particular, we show that, for an arbitrary modulation of the trapping potential, thermodynamic quantities of interest, such as the characteristic function $G_{\beta}(\vartheta)$, the moments $\langle W^{n}\rangle$ of the work distribution $P(W)$, and the nonadiabaticity parameter $Q^{*}(t)$, can all be evaluated using Fredholm determinants associated solely with single-particle wave functions. This is similar to the functional determinant approach used in the treatment of time-dependent perturbations of Fermi gases Abanin and Levitov (2005); d’Ambrumenil and Muzykantskii (2005); Knap et al. (2012); Schmidt et al. (2018), that relies on Levitov-Lesovik formula Levitov and Lesovik (1993) for the characteristic function of charge distribution in electron transport, even though our derivations do not rely on the said formula.

We then apply our general results to study nonequilibrium thermodynamics of a periodically (sinusoidally) modulated harmonically trapped TG gas. Periodically driven systems are rather common in nature and form an important class of problems in quantum dynamics, which motivates our choice and distinguishes it from sudden quenches studied before. For periodic modulation of a harmonically trapped TG gas, we are able to present explicit analytical expressions for the aforementioned thermodynamic quantities, and we illustrate as an example that by simply tuning the amplitude of the modulation one can significantly increase the amount of energy pumped into the system (equivalent to the amount of work done on an isolated system) in a given time. Our results can serve as a guide for designing quantum engine cycles based on periodic modulation rather than the more traditional sudden quenches (see, e.g., Jaramillo et al. (2016)) or other quantum control approaches such as shortcut to adiabaticity Deffner et al. (2014); Guéry-Odelin et al. (2019); Beau and del Campo (2020); Xu et al. (2020). Finally, our simple determinantal expression for the characteristic function allows us to explicitly calculate the expectation value $\langle e^{-\beta W}\rangle$ (where $\beta=1/k_{B}T$ is the inverse temperature) and hence demonstrate the validity of the quantum version of the Jarzynski equality Jarzynski (1997); Kurchan (2000); Esposito et al. (2009); Yi et al. (2012), which relates quantum work associated with a nonequilibrium process with the equilibrium free energy difference of the initial and final states of this quantum many-body system.

As an example of an interacting many-body system with determinantal structure of the quantum many-body function we consider the Tonks-Girardeau gas of $N$ bosons of mass $m$ interacting in one dimension via two-body hard-core interaction potential. Even though we are specifying the TG gas as the primary example for explicit illustration and discussion of our results, we note that all thermodynamic quantities calculated in this work are equally applicable to a 1D gas of spinless (or spin-polarized) noninteracting fermions Gong et al. (2014); Vicari (2019); Fre , as well as to a 1D gas of hard-core anyons Girardeau (2006); Pâţu et al. (2008); Pâţu (2020). In the purely bosonic or anyonic hard-core cases, the underlying statistics of the particles is governed by the constraints on the commutation relations of the field operators imposed by the hard-core diameter (with additional phase slips for anyons, compared to bosons), which at the same time allow for the Bose-Fermi or anyon-Fermi mapping to a pure ideal Fermi gas Girardeau (1960, 1965); Girardeau and Wright (2000); Pezer and Buljan (2007); Yukalov and Girardeau (2005); Girardeau (2006); Pâţu et al. (2008); Pâţu (2020). In all cases the evolution of the system in an external time-dependent one-body trapping potential $V(x,t)$ is governed by the free-fermion Hamiltonian

$$\hat{H}(t)=\sum_{j=1}^{N}\left[-\frac{\hbar^{2}}{2m}\frac{\partial^{2}}{\partial x_{j}^{2}}+V(x_{j},t)\right].$$

(1)

The TG or anyonic gas can equivalently be described as the limiting case of the Lieb-Lininger model with infinitely strong two-body $\delta$-function interaction potential Lieb and Liniger (1963).

The thermodynamics scenarios that we consider correspond to the system prepared at time $t=0$ in a thermal equilibrium with a reservoir at temperature $T$ and chemical potential $\mu$. It is then decoupled from the reservoir and evolves unitarily in an arbitrary time-dependent trapping potential $V(x,t)$. According to Bose-Fermi or anyon-Fermi mapping, the evolving many-body wave function of the system can be expressed in terms of the fermionic many-body wave function

$$\Psi_{s}(x_{1},...,x_{N};t)\!=\!A(x_{1},...,x_{N})\Psi_{s}^{F}(x_{1},...,x_{N};t),$$

(2)

where $A(x_{1},...,x_{N})$ is a statistical factor that ensures the right symmetry of the wave function under particle exchange, whereas the pure fermionic wavefucntion $\Psi_{s}^{F}(x_{1},\dots,x_{N};t)$ can be written as Slater determinant of the single-particle wave functions $\phi_{s_{i}}(x,t)$, where $\{s_{j}\}_{j=1,2,\dots,N}$ are a set of relevant quantum numbers. For a TG gas of hard-core bosons, $A(x_{1},...,x_{N})$ is a unit antisymmetric function given by $A(x_{1},\ldots,x_{N})=\prod_{1\leq j\leq i\leq N}\text{sgn}(x_{i}-x_{j})$, with $\text{sgn}(0)=0$, whereas for hard-core anyones it has an additional phase factor Girardeau (2006); Pâţu et al. (2008); Pâţu (2020). The relevant property that ultimately makes the calculation of our thermodynamic quantities identical to those of free fermions is that $|A(x_{1},...,x_{N})|^{2}=1$ and hence can be dropped from the final results; for the cases where $|A(x_{1},...,x_{N})|^{2}$ returns $0$, this corresponds to particles sharing the same position and is hence already taken into account by the fermionic Slater determinant which respects the Pauli exclusion principle.

With these introductory remarks in mind, we now turn to the discussion of thermodynamic work performed on the system by the external trapping potential $V(x,t)$. We can evaluate the quantum work using two projective energy measurements, one performed at time $t=0$ and the other at time $t$. The work in each realization of the experiment is given by the difference in the measured energies, and the associated work probability distribution function in an ensemble of realizations is given by Kurchan (2000); Talkner et al. (2007); Yi et al. (2012); Campisi et al. (2011); Jarzynski et al. (2015)

$\displaystyle P(W,t)\!=\!\sum_{N=0}^{\infty}\!\sum_{s}\!\sum_{{s^{\prime}}}p_{N,s}^{(0)}\,p_{{s^{\prime}}|s}^{(t)}\,\delta(W\!-\!E^{(t)}_{N,{s^{\prime}}}\!+\!E_{N,s}^{(0)}).$

(3)

In the following, we will explicitly omit the time dependence of the work distribution and conveniently write $P(W)\equiv P(W,t)$. Here, $p_{N,s}^{(0)}$ represents the probability of the first measurement returning an energy eigenvalue $E_{N,s}^{(0)}$ corresponding to the $N$-particle eigenstate $|\Psi_{s}(0)\rangle$ (where we omit the index $N$ for notational simplicity) of the initial Hamiltonian $\hat{H}(t\!=\!0)$, satisfying $\hat{H}(t\!=\!0)|\Psi_{s}(0)\rangle=E_{N,s}^{(0)}|\Psi_{s}(0)\rangle$. We consider initial thermal equilibrium states described by the grand-canonical ensemble, and therefore $p_{N,s}^{(0)}$ is given by the normalized Gibbs factor, $p_{N,s}^{(0)}=\frac{1}{\mathcal{Z}_{0}}e^{-\beta(E_{N,s}^{(0)}-\mu N)}$. Here, $\beta=1/k_{B}T$ is the initial inverse temperature (with $k_{B}$ the Boltzmann constant), $\mu$ is the initial chemical potential, and $\mathcal{Z}_{0}$ is the corresponding grand-canonical partition function. The system is then isolated from the reservoir and undergoes unitary evolution generated by $\hat{U}(t)=\mathcal{T}e^{-i\int_{0}^{t}dt^{\prime}\hat{H}(t^{\prime})/\hbar}$. The resulting state is $|\Psi_{s}(t)\rangle=\hat{U}(t)|\Psi_{s}(0)\rangle$. Next, a second projective energy measurement at time $t$ returns one of the instantaneous energy eigenvalues $E_{N,s^{\prime}}^{(t)}$, with the corresponding instantaneous eigenstate denoted via $|\Phi_{s^{\prime}}^{(t)}\rangle$, such that $\hat{H}(t)|\Phi_{s^{\prime}}^{(t)}\rangle=E^{(t)}_{N,{s^{\prime}}}|\Phi_{s^{\prime}}^{(t)}\rangle$. Therefore, the second probability entering into Eq. (3), $p_{s^{\prime}|s}^{(t)}$, is given by the transition probability from $|\Psi_{s}(t)\rangle$ to $|\Phi_{s^{\prime}}^{(t)}\rangle$, namely $p_{s^{\prime}|s}^{(t)}=|\braket{\Phi_{{s^{\prime}}}^{(t)}}{\Psi_{s}(t)}|^{2}=|\braket{\Phi_{{s^{\prime}}}^{(t)}}{\hat{U}(t)\Psi_{s}(0)}|^{2}$. Finally, the $\delta$-function in Eq. (3) ensures the conservation of energy, such that in each realization of the protocol the work is given by $W(t)=E_{N,s^{\prime}}^{(t)}-E_{N,s}^{(0)}$.

However, in practice, it is often more convenient to work with the characteristic function of the work distribution, which is given by the Fourier transform of $P(W)$,

$$G_{\beta}(\vartheta)=\int\mathrm{d}We^{i\vartheta W/\hbar}P(W)=\langle e^{i\vartheta W/\hbar}\rangle.$$

(4)

The characteristic function is the generating function of the moments of the distribution of work through successive differentiation,

$$\langle W(t)^{n}\rangle=(-i\hbar)^{n}\frac{d^{n}G_{\beta}(\vartheta)}{d\vartheta^{n}}\Big{|}_{\vartheta=0},$$

(5)

with the first moment corresponding to the mean work $\langle W(t)\rangle$.

Evaluating $G_{\beta}(\theta)$ for a generic quantum many-body system is a daunting task. However, for systems with the determinantal structure of the many-body wave-function as in Eq. (2), one can find numerically tractable expressions for $G_{\beta}(\theta)$ for a general $V(x,t)$. Since the mapping (2) holds true for the instantaneous eigenfunctions $\Phi_{s}^{(t)}(x_{1},...,x_{N})$ as well, the determinantal structure of the fermionic many-body wave function reduces the evaluation of Eq. (4) to that of evaluating integrals involving time-evolved single-particle wave functions $\phi_{s_{i}}(x,t)$ and the instantaneous single-particle eigenfunctions $\phi_{s_{i}}^{(t)}(x)$ with respective instantaneous eigenenergies $E_{s_{i}}^{(t)}$. Here, $\phi_{s_{i}}(x,0)$ are the energy eigenstates for the initial trapping potential $V(x,0)$ with eigenenergies $E_{s_{i}}^{(0)}$, such that the total energy of the $N$-particle system system is given by $E_{N,s}^{(0)}=\sum_{i=1}^{N}E_{s_{i}}^{(0)}$.

Using the properties of Fredholm integral equations and determinants Lenard (1966); Bornemann (2010), the characteristic function can be written as (see Appendix A for details),

$$G_{\beta}(\vartheta)=\frac{\mathcal{D}_{\hat{K}}^{k}(\vartheta)}{\mathcal{D}_{\hat{F}_{0}}^{f_{0}}}=\frac{\mathrm{det}(1+\hat{K})}{\mathrm{det}(1+\hat{F}_{0})},$$

(6)

where the superscripts $k$ and $f_{0}$ indicate the kernels of the respective Fredholm determinants, $\mathcal{D}_{\hat{K}}^{k}$ and $\mathcal{D}_{\hat{F}_{0}}^{f_{0}}$, given by

	$\displaystyle k(x,y)$	$\displaystyle=\sum_{i,j}\phi_{i}(x,t)k_{ij}(t)\left(\phi_{j}^{(t)}(y)\right)^{*},$		(7)
	$\displaystyle f_{0}(x,y)$	$\displaystyle=z\sum_{i}e^{-\beta E_{i}^{(0)}}\phi_{i}(x,0)\left(\phi_{i}^{(0)}(y)\right)^{*},$		(8)

with

$$k_{ij}(t)\!=\!ze^{-E_{i}^{(0)}\left(\beta+i\vartheta/\hbar\right)+i\vartheta E_{j}^{(t)}/\hbar}\!\int\!dw\phi_{i}^{*}(\omega,t)\phi_{j}^{(t)}(\omega),$$

(9)

and $z\!=\!e^{\beta\mu}$ denoting the fugacity. Here, the operators $\hat{K}$ and $\hat{F_{0}}$ are integral operators with kernels given by $k(x,y)$ and $f_{0}(x,y)$, respectively. For instance, the action of $\hat{K}$ on an arbitrary function $\xi(r)$ is given by $(\hat{K}\xi)(w)=\int_{\mathbb{R}}k(w,v)\xi(v)dv$.

Equation (6) is a significant result of this work as it allows us to express the characteristic function of the work distribution function of the quantum many-body system exclusively in terms of the single-particle quantities. This can be made more explicit if one rewrites Eq. (6) in terms of the eigenvalues of $\hat{K}$ and $\hat{F}_{0}$, which we denote by $\{\Lambda_{i}^{(\hat{K})}\}_{i=0,1,\dots}$ and $\{\Lambda_{i}^{(\hat{F}_{0})}\}_{i=0,1,\dots}$, respectively. The eigenvalues for the operator $\hat{K}$ (and similarly for $\hat{F}_{0}$) and their corresponding eigenfunctions $\{\theta_{i}^{(\hat{K})}(w)\}_{i=0,1,\dots}$ satisfy the equation

$$\int_{\mathbb{R}}dv~{}k(w,v)\theta_{i}^{(\hat{K})}(v)=\Lambda_{i}^{(\hat{K})}\theta_{i}^{(\hat{K})}(w),$$

(10)

with a similar equation for the operator $\hat{F}_{0}$, where $k(w,v)$ is replaced by $f_{0}(w,v)$. Using the expansion of the determinant of an operator as a product over its eigenvalues, Eq. (6) can be rewritten as

$$G_{\beta}(\vartheta)=\prod_{i}\left(\frac{1+\Lambda_{i}^{(\hat{K})}}{1+\Lambda_{i}^{(\hat{F}_{0})}}\right),$$

(11)

which expresses the characteristic function of the work distribution of the TG gas under an arbitrary driving protocol as an infinite product over the eigenvalues of two integral operators with kernels given by (7) and (8).

While the exact analytical evaluation of the eigenvalues is not always possible, the Fredholm determinant representation of the characteristic function, Eq. (6), offers a compact and a computationally practical way to efficiently evaluate many-body thermodynamic quantities of interest numerically, in terms of single-particle wave functions. For sufficiently smooth kernels, a simple tabulation of the kernel on a finite grid and the use of Nyström classical quadrature routine enables one to numerically evaluate the Fredholm determinant with small absolute errors Bornemann (2010). In Appendix B, we show that the eigenvalues $\Lambda_{i}^{(\hat{K})}$ (and $\Lambda_{i}^{(\hat{F}_{0})}$ ) and the determinant of the Fredholm integral equation can also be efficiently obtained through direct diagonalization of a matrix with elements proportional to the overlap between the time evolved and instantaneous eigenfunctions.

The above discussion is written generally for a thermal initial state. However, in many quantum simulation platforms, the relevant initial state is a pure state such as the many-body ground state of the initial Hamiltonian $|\Psi_{0}(0)\rangle$. In this regime, the sensitivity of the system to the driving protocol can be probed via the Loschmidt echo, whose amplitude is equivalent to the zero-temperature characteristic function $\mathcal{G}(t)=G_{\beta\to\infty}(t)$ Silva (2008). The Loschmidt echo amplitude is given by (see, e.g., Ref. Gorin et al. (2006) for a review)

$$\mathcal{G}(t)=\langle\Psi_{0}(0)|e^{i\hat{H}(0)t/\hbar}\mathcal{T}e^{-i\int_{0}^{t}dt^{\prime}\hat{H}(t^{\prime})/\hbar}|\Psi_{0}(0)\rangle,$$

(12)

with the Loschmidt echo itself given by

$$L(t)=|\mathcal{G}(t)|^{2}.$$

(13)

The Loschmidt echo gives the survival probability of the initial eigenstate $|\Psi_{0}(0)\rangle$ first evolved forward in time according to $\hat{H}(0)$ and then backward in time according to the dynamics generated by $\hat{H}(t)$. Its utility is in characterizing the survival of a quantum state when an imperfect time-reversal is applied. Since $\hat{H}(0)|\Psi_{0}(0)\rangle=E_{0}|\Psi_{0}(0)\rangle$ and $\mathcal{T}e^{-i\int_{0}^{t}dt^{\prime}\hat{H}(t^{\prime})/\hbar}|\Psi_{0}(0)\rangle=\hat{U}(t)|\Psi_{0}(0)\rangle=|\Psi_{0}(t)\rangle$ is the time-evolved state of the system, one can rewrite the Loschmidt echo as

$$L(t)=|\langle\Psi_{0}(0)|\Psi_{0}(t)\rangle|^{2}.$$

(14)

We now recognize this expression as the dynamical fidelity, $\mathcal{F}(t)=|\langle\Psi_{0}(0)|\Psi_{0}(t)\rangle|^{2}$, which is simply the squared overlap between the initial state and the time-evolved state of the system.

In order to derive a compact and computationally tractable expression for the Loschmidt echo $L(t)$, we use a similar procedure to that used in the derivation of Eq. (6) (see Appendix E). The echo amplitude $\mathcal{G}(t)$ can then be written as the determinant of an $N\times N$ matrix $\mathbf{A}$ (where $N$ is the number of particles in the system) containing the overlaps between the initial and time evolved single-particle eigenfunctions $\phi_{n}(x,t)$ [see Eq. (17) below], and hence

$$L(t)=|\mathrm{det}\mathbf{A}(t)|^{2},$$

(15)

where

$$A_{mn}(t)=\int_{-\infty}^{\infty}\phi_{m}^{\ast}(x,0)\phi_{n}(x,t)\,dx,$$

(16)

and $m,n=0,1,\dots,N-1$. This determinantal result for the Loschmidt echo reproduces the expression derived previously for the TG gas in Refs. del Campo (2011); Lelas et al. (2011, 2012); Pons et al. (2012). An equivalent determinantal expression for a system of 1D free fermionic chains has recently been derived in Ref. Gamayun et al. (2020).

In this section, we apply the Fredholm determinant formalism to a time-varying harmonic potential $V(x,t)=m\omega^{2}(t)x^{2}/2$. This model offers the advantage of being analytically solvable, in addition to being typical of ultracold atom experiments, and thus provides an ideal platform for studying the work distribution and other thermodynamic quantities of interest. We point out that the thermodynamic quantities for the TG in a harmonic trap with arbitrary time-variation of the trap frequency $\omega(t)$ were previously derived by Jaramillo et al. in Ref. Jaramillo et al. (2016) using the scale invariance of the corresponding many-body wave function. However, our determinantal approach allows us to derive a more general result for the mean work $\langle W(t)\rangle$ and the nonadiabaticity parameter $Q^{*}(t)$ that is valid for an arbitrary spatial shape of $V(x,t)$ beyond the simple harmonic trapping. An immediate application of our general results to this limit—as is done in this section—reproduces their results for the harmonic trapping.

Our numerical results in the harmonic trapping regime are further distinguished from previous studies due to the use of a nontrivial driving protocol wherein the trap potential is modulated sinusoidally in time  with $V(x,t)=m\omega^{2}(t)x^{2}/2$ and $\omega(t)^{2}=\omega(0)^{2}[1-\alpha\sin(\Omega t)]$. Here, $\Omega$ is the modulation frequency and $\alpha$ characterizes the modulation amplitude. Under such periodic modulation, the TG trap displays a rich phase diagram in the $(\Omega,\alpha)$ parameter space Atas et al. (2019), including regions of stable (bounded) and unstable (exponentially growing due to parametric resonance) dynamics Quinn and Haque (2014). The latter regime is particularly advantageous for creating highly nonequilibrium states, with very large values of the nonadiabaticity parameter $Q^{*}(t)\gg 1$. This in turn corresponds to large amount of work that can be done on the system accompanied with modest changes in the trap size, which may prove useful in optimising the performance of a quantum refrigerator with the TG gas serving as the working medium.

With an arbitrary time-varying harmonic potential, the single-particle Schrödinger equation is exactly solvable and the time evolved single-particle wave functions $\phi_{n}(x,t)$ can be obtained through a simple scaling transformation Perelomov and Zel’dovich (1998); Minguzzi and Gangardt (2005) given by

$$\phi_{n}(x,t)=\frac{1}{\sqrt{\lambda}}\phi_{n}\left(\frac{x}{\lambda},0\right)\exp\left[i\frac{mx^{2}}{2\hbar}\frac{\dot{\lambda}}{\lambda}-i\frac{\mathcal{E}_{n}(t)\,t}{\hbar}\right].$$

(17)

Here, the initial wave functions $\phi_{n}(x,0)$ are the Hermite-Gauss polynomials of a quantum mechanical $1$D harmonic oscillator,

$$\phi_{n}(x,0)=\frac{\exp\left(-x^{2}/2l^{2}_{\mathrm{ho}}(0)\right)H_{n}(x/l_{\mathrm{ho}}(0))}{\pi^{1/4}\sqrt{2^{n}n!l_{\mathrm{ho}}(0)}},$$

(18)

with frequency $\omega(0)$, energy eigenvalues $E_{n}^{(0)}=\hbar\omega(0)(n+1/2)$, and harmonic oscillator length $l_{\mathrm{ho}}(0)=\sqrt{\hbar/m\omega(0)}$. Furthermore, the scaling function $\lambda(t)$ is a solution to the Ermakov-Pinney equation erm ; Atas et al. (2019),

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

(19)

with initial conditions $\lambda(0)$=1 and $\dot{\lambda}(0)=0$, whereas the time-dependent phase factor $\mathcal{E}_{n}(t)$ in Eq. (17) is defined in terms of $\lambda(t)$ and reads as

$$\mathcal{E}_{n}(t)=\hbar\omega(0)\left(n+\tfrac{1}{2}\right)\,\frac{1}{t}\int_{0}^{t}\frac{dt^{\prime}}{\lambda^{2}(t^{\prime})}.$$

(20)