Extreme dynamics and relaxation of quantum gases: A hydrodynamic approach

TL;DR

This paper uses hydrodynamics to analyze the extreme out-of-equilibrium behavior of quantum gases, revealing self-similar solutions, shock formation, and relaxation dynamics that align with experimental observations.

Contribution

It introduces a hydrodynamic framework for studying the dynamics of quantum gases, deriving scaling laws and demonstrating shock formation and relaxation phenomena.

Findings

Existence of self-similar solutions with specific scaling exponents.

Development of shock fronts during gas expansion.

Steady states in relaxation dynamics match analytical predictions.

Abstract

The evolution of quantum gases, released from traps, are studied through hydrodynamics, both analytically and numerically, in one and two dimensions. In particular, we demonstrate the existence of long time self-similar solutions of the Euler equations, for the density and velocity fields, and derive the scaling exponents as well as the scaling functions. We find that the expanding gas develops a shock front and the size of the cloud grows in time as a powerlaw. We relate the associated exponent to that appearing in the corresponding equation of state of the quantum gas. Furthermore, we study the relaxation dynamics of a trapped quantum gas and show that the resulting steady state is in excellent agreement with that derived analytically. Our hydrodynamic approach is versatile and can be used to unravel several other far-from-equilibrium collective phenomenon of extreme nature, relevant to the growing experimental interests in quantum gases.

Equations

$$\frac{\partial \rho }{\partial t}+\nabla \cdot (\rho \mathbf{v})=0,$$

$$\frac{\partial \mathbf{v}}{\partial t}+\mathbf{v}\cdot \nabla \mathbf{v}+\nabla w=-\nabla U_{\mathrm{trap}},$$

$$w(\rho )=D{\rho }^{\delta }.$$

$$\frac{\partial \rho }{\partial t}+\frac{\partial \rho v}{\partial x}$$

$$\frac{\partial v}{\partial t}+\frac{\partial }{\partial x}\left(\frac{v^2}{2}+D{\rho }^{\delta }\right)$$

$$\rho /{\rho }_0\to \rho ,x/\sigma \to x,\text{leavevmode}\frac{v}{(D{\rho }_0^{\delta }{)}^{1/2}}\to v,\text{leavevmode}\frac{t(D{\rho }_0^{\delta }{)}^{1/2}}{\sigma }\to t.$$

$$\rho (x,t)=t^{-b}f(\xi ),v(x,t)=t^{-c}g(\xi ),\mathrm{where}\text{leavevmode}\text{leavevmode}\xi =\frac{x}{t^b},$$

$$b=\frac{2}{\delta +2},c=\frac{\delta }{b+\delta }.$$

$$\frac{d(fg)}{d\xi }-b\xi \frac{df}{d\xi }-bf$$

$$\frac{d{f}^{\delta }}{d\xi }+\frac{1}{2}\frac{d{g}^2}{d\xi }-b\xi \frac{dg}{d\xi }-cg$$

$$R(t)={\xi }_f{t}^b,$$

$$U=\dot{R},$$

$$v_{-}=U,\text{leavevmode}\text{leavevmode}\text{leavevmode}\frac{v_{-}^2/2+{\rho }_{-}^{\delta }}{v_{-}}=U,$$

$$f_{-}={\left(\frac{{\xi }_f^2{b}^2}{2}\right)}^{1/\delta },\text{leavevmode}\text{leavevmode}\text{leavevmode}{g}_{-}={\xi }_{f}b.$$

$$f(\xi )$$

$$g(\xi )$$

$${\int }_0^{{\xi }_f}d\xi f(\xi )=M,$$

$${\xi }_f={\left(\frac{(2-\delta {)}^{-1/\delta }(\delta +2{)}^{2/\delta }M}{{}_2{F}_1\left(\frac{1}{2},-\frac{1}{\delta };\frac{3}{2};\frac{\delta }{\delta -2}\right)}\right)}^{\frac{\delta }{\delta +2}},$$

$$\rho (x,0)$$

$$v(x,0)$$

$$\frac{\partial \Delta \rho }{\partial t}+{\rho }_B\frac{\partial \Delta v}{\partial x}=0,\frac{\partial \Delta v}{\partial t}+\delta {\rho }_B^{\delta -1}\frac{\partial }{\partial x}\left(\Delta \rho \right)=0,$$

$$c_s=\sqrt{\delta {\rho }_B^{\delta }}.$$

$$\rho ={\rho }_B+f(x/t),\text{leavevmode}\text{leavevmode}v=g(x/t),$$

$$-\xi \frac{df}{d\xi }+{\rho }_B\frac{dg}{d\xi }+\frac{dfg}{d\xi }$$

$$-\xi \frac{dg}{d\xi }+g\frac{dg}{d\xi }+\delta ({\rho }_B+f{)}^{\delta -1}{f}^{\prime }$$

$$(g-\xi )\frac{df}{d\xi }+(f+{\rho }_B)\frac{dg}{d\xi }$$

$$\delta (f+{\rho }_B{)}^{\delta -1}\frac{df}{d\xi }+(g-\xi )\frac{dg}{d\xi }$$

$$g=\xi \pm \sqrt{\delta }(f+{\rho }_B{)}^{\delta /2}.$$

$$\pm (1+\delta /2)\sqrt{\delta }(f+{\rho }_B{)}^{\delta /2}{f}^{\prime }+(f+{\rho }_B)$$

$$f(\xi )$$

$$g(\xi )$$

$$\Phi=\left(\frac{2(\delta+1)\text{{\hbox{\rho_B}}}^{\delta/2}-\delta^{1/2}\xi}{\delta+2}\right)^{2/\delta}.$$

$$f(\xi =c_s)=0,g(\xi =c_s)=0.$$

$$\displaystyle f(\xi-c_{s})=A(\xi-c_{s})+O\big{(}(\xi-c_{s})^{2}\big{)}\,,$$

Full text

Turbulence

I Introduction

The non-equilibrium dynamics of cold atomic gases has been a subject of great interest recently both experimentally and theoretically [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]. Fascinating experiments have been performed [11, 12, 13, 14, 15] to observe the free expansion and relaxation of trapped quantum gases and there has been significant progress in their theoretical understanding through the framework of hydrodynamics [16, 17] and generalized hydrodynamics [18]. While many of the observed experimental features can be explained by Euler hydrodynamics, the role of viscosity and its quantum aspects has been explored in experiments [19, 20, 21] and some features are now understood from theory [22, 23, 24]. However, the exact form of the viscous dissipation in the hydrodynamic equation is far from known.

In the context of dynamics in quantum gases, there has been immense activity in bosons [25], fermions [26] and various mixtures [27, 28, 29]. Such dynamics shines light on the consequences of the underlying interactions that lead to exotic collective phenomenon. There has been studies in understanding dynamics near an equilibrium background (such as sounds waves [30, 31, 32, 33, 34, 35, 36]) and dynamics very far from equilibrium (for e.g., resulting in shock waves [14, 17, 37, 38] or energy dynamics [39]). However, in these systems dynamics of extreme nature is far less understood. One such specific problem is to study an evolution that mimics a blast that can arise from a sudden release of a substantial amount of particles or energy from a localized region. Some natural challenges while studying these problems include: (i) establishing rigorously the collective hydrodynamic description (ii) unavailability of analytical solutions, especially very far from equilibrium even if a phenomenological hydrodynamic theory is available.

At the level of Euler hydrodynamics, one of the most interesting results are the so-called self-similar scaling solutions [40], first discussed by Taylor, von Neumann and Sedov [41, 42, 43, 44] in the context of atomic explosions. The main question there is to understand the time evolution of the hydrodynamic fields, namely density, momentum and energy, for an initial state where a huge amount of localized energy is released in an otherwise cold ambient gas — this is referred to as the blast problem. Following the classic work of TvNS, there has been much recent interest in the blast problem, in particular looking at direct comparisions of the Euler predictions with results from microscopic molecular dynamics simulations and understanding the effects of dissipation [45, 46, 47, 48, 49, 50, 51, 52]. So far, all work has focused on classical systems and the main goal of the present paper is to examine self similar solutions for blasts in quantum gases. Furthermore we consider quantum gases confined in traps and study the approach to equilibrium using the framework of hydrodynamics.

While the physics of blasts is rooted in the fascinating, collective dynamics of the constituent bosons or fermions, it is possible to capture the essential, experimentally measurable dynamics through a coarse-grained hydrodynamic description. Here we consider the zero temperature case where it is sufficient to consider the equations for the density field, $\rho({\bf x},t)$, and velocity field, ${\bf v}({\bf x},t)$. The hydrodynamic Euler equations, in the presence of an external trapping potential, $U_{\rm trap}({\bf x})$, are given by

[TABLE]

where the specific enthalpy $w$ is related to the pressure, $P$, via the relation $\partial_{\rho}P=\rho\,\partial_{\rho}w$. A knowledge of the zero-temperature equation of state of the quantum gas thus enables us to write the hydrodynamics equations. Interestingly, a large number of quantum gases, both bosonic as well as fermionic, have a simple power-law form for the specific enthalpy

[TABLE]

As examples in one-dimensional ($d=1$) systems, we have: $\delta=1$ and $\delta=2$ for the Lieb-Liniger gas of bosons at weak and strong coupling respectively; $\delta={2}/{5}$ is the limiting case of a quasi-1D (cigar-shaped) unitary Fermi gas [14, 17]; $\delta=2/3$ for unitary Fermi gas in plane-wave limit. In dimensions $d>1$ examples include: $\delta=2/d$ for weakly interacting fermions at zero temperature. On the other hand, the unitary Fermi gas in $d$ dimensions also shares the same dependence on $\delta$. However, there is an overall renormalization in comparison to weakly interacting case [14, 37]. Remarkably, the unitary fermi gas, being a strongly correlated quantum fluid, shares deep resemblance to very different systems such as the quark–gluon plasma [53, 54, 2] thereby enabling a hydrodynamic description.

It is worth noting that, beyond cold gases, the well-known Riesz family of gases [55, 56, 57, 58, 59] provides an important and distinctly complementary example described by Eqs. (1,2). This family comprises of particles with all-to-all interactions of power-law type characterized by an exponent, say $s$. The collective theory of the Riesz family of gases in a broad range of parameter space falls under the general structure of Eqs. (1a) and (1b). In particular, a class of these models are of finite-ranged type [57]. It was shown in Ref 57 that a finite-ranged version of the Riesz gas in one dimension characterised by exponent $s$, would correspond to $\delta=s$ [Eq. (2)] for any $s>0$. This immediately offers us a plethora of possibility of realization of the same collective dynamics for a wide window of $\delta$.

The precise set-up where we study blasts in quantum gases is motivated by experiments in cold atoms using, e.g, time-of-flight (TOF) techniques, a standard tool in atomic physics for probing momentum and coherence in quantum gases. In TOF experiments, atoms are cooled in a (typically harmonic) trap, abruptly released, and allowed to expand for a set time, and then imaged via absorption imaging techniques [60, 61, 62]. In the present work we study, through Euler hydrodynamics (Eqs. 1), the dynamical evolution of such gases after their release from the trap. We show that, at long times, they evolve into self-similar forms and provide a complete analytic and numerical understanding. Secondly, we also consider the problem of a quantum gas in a trap ($U_{\rm trap}\neq 0$) and study the relaxation dynamics to the steady state distribution. In this case, we add the Navier-Stokes corrections to Eqs. (1). For reasons explained later, we restrict ourselves to $0<\delta<2$ for the one dimensional case and $0<\delta<1$ for the two dimensional case.

This paper is organized as follows. In Sec. II, we discuss the one-dimensional gas without any external trap. Both evolution in zero (Sec. II.1) and finite background (Sec. II.2) are discussed. We then discuss the relaxation dynamics of one-dimensional gas in a trap in Sec. III. The two dimensional case is discussed in Sec. IV. We end with conclusions along with an outlook in Sec. V. Details of direct numerical simulations are provided in Appendix. A.

II Free evolution in one-dimensional case

We begin with the one-dimensional (1D) problem with the trap switched off and consider the Euler equations for the density and velocity fields:

[TABLE]

We consider initial condition where we start with a density profile on top of a uniform density background $\rho_{B}$ and zero velocity field. The cases with background density either zero and non-zero will be considered separately in the following sections. Let us assume that the density profile has a spatial spread over a characteristic length $\sigma$ and a typical size $\rho_{0}$. We can then transform to dimensionless variables

[TABLE]

This corresponds to setting $D=1$ in Eq. (3b) and we will henceforth assume this.

We will consider below the two case where an excess of particles is introduced in a localized region, either in an empty region (Sec. II.1) or one with a finite density of particles (Sec. II.2).

II.1 Evolution in a vacuum

In this section, the initial condition is a spatially localized mass (with centre of mass at the origin) of a gas with $\rho_{B}=0$ and net zero momentum. It is expected that the Euler equations should, at long times, develop a shock and the hydrodynamic fields have scaling forms. At long times the fields $\rho,v$ will therefore take the forms

[TABLE]

$f$ and $g$ are scaling functions and $b$ and $c$ are scaling exponents. Plugging the forms in Eq. (5) into Eq. (3b) and requiring that the resulting equations in $f(\xi)$ and $g(\xi)$ have no explicit time-dependence fixes the scaling exponents $b$ and $c$ as

[TABLE]

We then get the following ordinary differential equations (ODEs) for the scaling functions $f$ and $g$:

[TABLE]

Let us indicate the location of the shock in scaled variable as $\xi_{f}$. Therefore, it follows that the position of the shock is given by

[TABLE]

Let us denote the shock velocity

[TABLE]

and the fields $(\rho_{-},v_{-})$ and $(\rho_{+}=0,v_{+}=0)$ just behind and in front of the shock, respectively. The Rankine–Hugoniot [50, 51] boundary conditions then give:

[TABLE]

and hence

[TABLE]

With these boundary conditions, we can solve Eq. (7) to get

[TABLE]

The factor $\xi_{f}$ is fixed by using the mass conservation equation

[TABLE]

which, for $\delta<2$, gives

[TABLE]

where ${}_{2}F_{1}$ represents the Hypergeometric function. Equation. (14) shows that we can get an explicit form of shock location $\xi_{f}$ for $\delta\leq 2$. One can immediately see that the case $\delta=1$ gives $\xi_{f}=3(M/4)^{1/3}$.

We now test the scaling solutions in Eqs. (12a) and (12b) through direct numerical solution of the hydrodynamics equations. The details of the pseudospectral numerical method employed is provided in Appendix. A. We chose extremely localised yet generic, initial conditions. For a system of length $L$,

[TABLE]

In Fig. 1 we show results from direct numerical simulation for the time evolution of density and velocity for $\delta=2/5$. In panel (a), we show relatively early time evolution and the subsequent emergence of the shock. In panels (b) and (c) we show the late time evolution of density and velocity and their anticipated scaling and fit to the analytical predictions for the scaling functions given in Eqs. (12a) and (12b).

II.2 Evolution in a finite density background

We now consider the case where we add mass to an existing finite background density $\rho_{B}$. The localized mass density will then propagate and simultaneously spread. At sufficiently large times, the excess density is a perturbation and one can then do a linearized theory, which would predict the propagation to be ballistic characterized by a speed of sound $v_{s}$. If the initial excess density is large then the initial propagation will be sub-ballistic. However, after a considerable time has passed, the excess density will become small compared to the background thereby validating a linearized theory. The linearized theory can be obtained from Eq. (3) as (setting $D=1$)

[TABLE]

where $\Delta\rho$ and $\Delta v$ denote small perturbations to justify linearization. Equation (17) immediately gives the sound speed

[TABLE]

We now discuss the scaling solutions at long times. We expect ballistic behaviour at long times and hence we assume the form

[TABLE]

where $f(x/t)$ and $g(x/t)$ are scaling functions that we seek to compute analytically. Plugging Eq. (19) into Eq. (3) yields

[TABLE]

where we recall that $\xi=x/t$. Equation (20) can be recasted as

[TABLE]

Equation (21) is a linear homogeneous equation in $f^{\prime}$ and $g^{\prime}$ and therefore a solution exists only if the determinant $(g-\xi)^{2}-\delta(f+\rho_{B})^{\delta}=0$. This immediately gives us a relation between $g$ and $f$

[TABLE]

Plugging Eq. (22) into Eq (21a) one gets

[TABLE]

This can be integrated to yield (the physical solution)

[TABLE]

where,

[TABLE]

Note that, as requried, the following conditions are satisfied

[TABLE]

In Fig. 2(a) we show the density (upper panel) and velocity (lower panel) profiles for different times with the background density $\rho_{B}=1$ with $\delta=2/5$. Clearly, both fields evolve ballistically. After some time, the shock front $R(t)$ moves with the sound speed, i.e., $R(t)=c_{s}t$ where we recall that $c_{s}$ is given in Eq. (18). Given the symmetric nature of the evolution, we track the shock moving towards the right. In Fig. 2(b) we show the plot in scaled and shifted variables ($\xi-c_{s}$) where $\xi=x/t$. In Fig. 2(c) we plot the scaling exponent $b$ as a function of background density $\rho_{B}$. In the large system size limit, we observe exponent $b=1$ independent of $\delta$. Clearly, as $\rho_{B}\to 0$, we recover the freely-evolving gas limit where, as seen earlier $b=\frac{2}{2+\delta}$. In the inset of Fig. 2(c) we plot the position of shock and the ramp as a function of time. The thick black line, as a guide to the eye, confirms the ballistic motion of $R(t)$. We further analyze the shock front in Fig. 3. We demonstrate that as one increases the system size, the numerical profile approaches the asymptotic limit and converges to the analytical expression given in Eq. (24).

We next discuss some analytical estimates for the location of the shock and its scaling with time. For this, we revisit the expression of $f(\xi)$ in Eq. (24a). Note that for relatively small $\xi-c_{s}$, we have

[TABLE]

where

[TABLE]

Next we require by normalization,

[TABLE]

where we recall that $M$ is the excess mass and the factor of $2$ is because we are considering one side of the profile. Simplifying Eq. (29) by using Eq. (27) gives

[TABLE]

We next evaluate the height of the shock at $\xi=c_{s}+x_{\text{max}}/t$. To do so, we plug the value of $x_{\text{max}}/t$ in Eq. (27) to give

[TABLE]

We have numerically that verified that $f(\xi-c_{s}=x_{\rm max}/t)\sim 1/\sqrt{t}$. In fact, our numerical results are also consistent with the precise coefficient $\sqrt{M\,A}$ given in Eq. (31).

III Relaxation of a one-dimensional trapped gas

So far, we have discussed cases when there is no external confinement. In this section, we discuss the case when we have a confining potential $U_{\rm trap}\neq 0$. Starting from some nonequilibrium initial condition we study how the gas evolves and eventually reaches a steady state. The hydrodynamic equations are given by

[TABLE]

where we are using dimensionless variables as in Eq. (4) and $\omega$ and $\nu$ are also made dimensionless via

[TABLE]

Note that a more realistic form for the dissipative term in Eq. (33) would be $\partial_{x}(\rho\partial_{x}v)/\rho$. For simplicity, we replace it with the form $\nu\partial^{2}_{x}v$ which is expected to be a good approximation for shallow traps and avoids numerical issues related to divergences associated with the $1/\rho$ factor.

It is easy to see that the steady state solutions of Eq. (32) takes the form

[TABLE]

where the constant $a$ is fixed by the constraint

[TABLE]

In Fig. 4, we show results from direct numerical simulations of Eq. (32) for $\delta=2/5$. Since we are using a pseudospectral method with Fourier basis which is suited for periodic boundaries we smoothen the quadratic potential near the boundary and make it periodic to avoid numerical errors (see Appendix. A). In other words, without the smoothening, $\frac{\partial U_{trap}}{\partial x}$ is discontinuous at the boundary making the pseudospectral method ill-suited. We show in Fig. 4(a) density and (b) velocity profiles at several different times. Clearly at late time both the density and velocity fields converge to analytical expressions given in Eq. (34). Figure 4(c) displays an exponential decay in $\rho$ and $v$ at a fixed spatial location. We analytically derive the corresponding decay timescale later in this section.

As can be noticed in Fig. 4 (a) and (b), both the density and velocity profiles show damped oscillatory relaxation to the steady state. We now provide an in-depth understanding of this transient behaviour and, in particular, compute the associated damping rate and frequency. A full solution of the time-dependent hydrodynamic equations is challenging. However, for initial conditions that are close to the steady state one can consider a linearized dynamics. As we shall now show, such equations are amenable to analytic treatment and in fact also provide insights into the late time dynamics of even far-from equilibrium initial conditions.

To proceed further, we make the following substitution into Eq. (32)

[TABLE]

where recall that $\rho_{\rm ss}$ is given in Eq. (34). Retaining terms linear in the perturbation we obtain

[TABLE]

Note that Eq. (37) is linear in $\Delta\rho$ and $\Delta v$ but with space dependent coefficients making them still challenging to solve. Nevertheless we can in fact obtain the spectrum completely. Let us look for solutions of Eq. (37) of the form

[TABLE]

Plugging this into Eq. (37) gives us the following equations for $\psi$ and $\phi$:

[TABLE]

Plugging in the explicit form of $\rho_{\rm ss}$ in Eq. (34) into Eq (39) leads to the form:

[TABLE]

where,

[TABLE]

These equations have to be solved with the boundary conditions

[TABLE]

This will give us to a set of eigenfunctions $\{\psi_{n},\phi_{n}\}$ and corresponding allowed eignevalues $\{\lambda_{n}\}$ for our linear equations. From the structure of the eigenvalue equations [Eq. (40)], it is clear that the eigenvalues and eigenfunctions come in complex conjugate pairs. We now note that Eq. (41) is in the form of the Legendre equation whose general solution can be expressed in terms of the associated Legendre functions of the first and second kinds:

[TABLE]

where

[TABLE]

and $A_{1},A_{2}$ are arbitrary constants. Without loss of generality we can set $A_{1}=1$ since this can be fixed by normalization. We then have two complex unknown parameters, namely $A_{2}$ and $\lambda$, These are determined by imposing the two boundary conditions $\psi_{n}(\pm a)=0$. We solve for these unknown parameters numerically using Mathematica.

The general solution of Eq. (37) is given by

[TABLE]

where the constants $p_{n}$ are determined from the initial conditions. At large times, the eigenvalue with the smallest value for the real part will dominate. There are two such eigenvalues which we denote as

[TABLE]

where $\gamma>0$ and $\Omega$ are real. These provide us with the time scale for decay and the frequency of the oscillations, namely $\gamma^{-1}$ and $\Omega$ respectively. The expected forms for the density and velocity fields at late times can also be determined. Writing the dominant eigenfunctions [Eq. 49] in the form

[TABLE]

we have the following late time forms:

[TABLE]

where $C_{1},C_{2},\beta_{1},\beta_{2}$ are constants that are fixed by the initial conditions. In the numerics for the nonlinear case these constants can be fixed by using the value of $\Delta\rho$ and $\Delta v$ at some late time snapshot because in the early time there might be contribution from other eigen-functions which may lead to an incorrect estimation of the constants. In Fig. (5), we compare the analytical predictions [Eqs. (52), (53)] with direct numerical solutions of the linearized equations [Eq. (37)] and find excellent agreement. To further validate the calculations with the linearized equations, in Fig. 6 we explicitly compare numerical simulations of the full nonlinear system [Eq. (32)] with those of the perturbation equations [Eq. (37)]. For small perturbations, the linearized equations accurately approximate the full nonlinear dynamics. However, as the perturbation strength increases, deviations become apparent, indicating the breakdown of the linear approximation.

IV Two-dimensional case

In this section, we discuss the evolution, in two dimensions, of an initially localized mass of a quantum gas with zero total momentum. We only consider expansion into a vacuum. The initial mass distribution is assumed to have radial symmetry and so we look for solutions of the hydrodynamic equations which have this symmetry. In particular, we need to only consider the density field $\rho$, and the radial component of the velocity vector field $v$. Clearly, the space-time dependence of the fields will have the form $\rho=\rho(r,t),v=v(r,t)$ and the corresponding Euler equations in radial coordinates are

[TABLE]

As in Sec. II.1 we expect that asymptotically the dynamics develops a scaling form as

[TABLE]

where the factor, $2b$, for the density scaling appears since $\int_{0}^{\infty}drr\rho(r,t)$ has to be conserved. Plugging Eq. (55) into Eq. (54) and requiring that any explicit time dependence goes away, we get

[TABLE]

and the corresponding ordinary differential equations for the scaling functions

[TABLE]

We now discuss the Rankine-Hugoniot conditions which turn out to be identical to the 1D case given by Eqs. (10) or (11). To derive the first Rankine–Hugoniot condition, we integrate the first conserved quantity, $\int r\rho\,{\rm d}r$, across the shock at position $R$, i.e., over the interval $[R-\Delta,R+\Delta]$ and then take the limit as $\Delta\to 0$. This yields

[TABLE]

which implies

[TABLE]

where $\rho_{-}$ and $v_{-}$ denote the radial density and velocity before the shock, and $\rho_{+}$ and $v_{+}$ denote the density and velocity after the shock. Noting that $\rho_{+}=v_{+}=0$ gives the first Rankine–Hugoniot condition

[TABLE]

To derive the second Rankine–Hugoniot condition, we use the second conserved quantity, $\int rv\,\mathrm{d}r$, and integrate it around the shock

[TABLE]

which implies

[TABLE]

We now simplify the second term

[TABLE]

Plugging Eq. (IV) back in Eq. (62) and taking the limit $\Delta\to 0$ gives us the second Rankine-Hugoniot condition

[TABLE]

By using Eq. (60) in Eq. (64), we get

[TABLE]

In terms of scaled variables given in Eq. (55), we get

[TABLE]

which takes into account that the shock position is given by $R=\xi_{f}t^{b}$ and shock velocity $U=\dot{R}$.

We now want to find the solutions to Eq. (56). We take an Ansätz $g(\xi)=k\,\xi$ and substitute it in Eq. (57a) yielding

[TABLE]

To ensure a non-zero solution for Eq. (67), we need $k=b$. So, $g(\xi)=b\,\xi$. We verify via direct numerical simulations that this is indeed the case. Now we substitute $g(\xi)=b\,\xi$ in Eq. (57b) and get

[TABLE]

which simplifies to

[TABLE]

Integrating Eq. (69) from $\xi_{f}$ to $\xi$ we get,

[TABLE]

Finally, the analytical solution for $f$ and $g$ can be summarized as

[TABLE]

where we recall that $b$ is given in Eq. (56) and $f(\xi_{f})$ can be obtained from the Rankine-Hugoniot condition [Eq. (66)]. What remains is to find $\xi_{f}$ from mass conservation

[TABLE]

which yields (for $\delta<1$)

[TABLE]

Next we test our analytical prediction via direct numerical simulations. First we solve the full two-dimensional equations i.e. Eqs (1a) and (1b) for $\delta=2/5$ and $U_{trap}=0$ starting from a localised density profile and zero velocity. Figure. 7 shows the density evolution for early and late times demonstrating blast-like phenomenon and the emergence of a scaling form. In Fig. 8, we show the numerical results obtained from solving Eq. (54) numerically for $d=2$ and $\delta=2/5$. (a) Figure 8(a) shows the early time behavior and upper half of panels (b) and (c) show the radial density and velocity for relatively late times. In the lower panels pf Fig. 8 (b) and (c) we show excellent agreement between rescaled numerical profiles and analytical scaling functions given in Eq. (71).

V Conclusions and Outlook

To summarize, we have investigated quantum dynamics in highly out-of-equilibrium situations, both in one and two dimension, via hydrodynamics and using both analytical calculations and direct numerical simulations. We demonstrated that at long times, the hydrodynamic fields evolve to self-similar forms, that are independent of the details of the localized initial conditions. We derived analytically the scaling exponents and scaling functions and thus characterized completely the features of this highly nonequilibrium process. In particular, for a gas with the equation of state $P\sim\rho^{\delta+1}$ and evolving in a vaccuum, the shock front grows with time as $R(t)\sim t^{b}$, where $b=2/(2+\delta)$ and $1/(1+\delta)$ in one and two dimensions respectively. We also investigated the relaxation dynamics of a confined quantum gas. At late times we see a damped oscillatory relaxation, where the damping coefficient and frequency can be computed analytically from a linearized theory. We validated our analytical findings by extensive numerical simulations and found excellent agreement between the two.

Our methodology is highly adaptive and expected to be of relevance in understanding far from equilibrium dynamics in a variety of interacting quantum systems. It will be interesting to explore quantum gases in higher spatial dimensions, relaxing the assumptions of symmetry [64, 65, 66], which would make the effective one dimensional description ill-suited. It would also be important and interesting to extend our techniques to open quantum systems such as dissipative quantum fluids [67, 68]. Beyond cold quantum gases, our collective formalism and subsequent analysis of self-similar solutions is also relevant in several other contexts such as finite-ranged Riesz gas [57] and hydrodynamics of a plasma of strongly interacting quarks and gluons produced in relativistic heavy ion collisions [53, 54, 69, 2]. As noted earlier, our protocol is well-suited to TOF measurements, where atoms are cooled, released, and imaged via absorption [60, 61, 62]. Testing our predictions is thus highly feasible and the observation of self-similar scaling solutions seems like a fascinating possibility.

Acknowledgements

We thank P. L. Krapivsky for very helpful discussions. AD and MK thank the VAJRA faculty scheme (No. VJR/2019/000079) from the Science and Engineering Research Board (SERB), Department of Science and Technology, Government of India. AD acknowledges the J.C. Bose Fellowship (JCB/2022/000014) of the Science and Engineering Research Board of the Department of Science and Technology, Government of India. AD and MK thank the hospitality of Yukawa Institute for Theoretical Physics, Kyoto University during the workshop on “Hydrodynamics of low-dimensional interacting systems: Advances, challenges, and future direction”. The authors acknowledge the support of the DAE, Government of India, under projects nos. 12-R&D-TFR-5.10-1100 and RTI4001.

Appendix A Direct Numerical Simulations

We perform direct numerical simulations using a fully-dealiased pseudospectral method with periodic boundary conditions in both one and two dimensions. Spatial discretization is carried out using a Fourier basis, and time integration is performed via either the classical fourth-order Runge–Kutta (RK4) method or the integrating factor RK4 (IFRK4) scheme, depending on the stiffness of the equations [70].

To simulate the inviscid equations, we introduce a small artificial viscosity in both the density and velocity equations to regularize sharp gradients and suppress numerical instabilities. More precisely, in the simulations we replace the Euler equations in Eqs. (1a)-(1b) and work with the following dissipative equations

[TABLE]

For the range of values that we consider in this paper, i.e, $0<\delta<2$ in 1D and $0<\delta<1$ in 2D, one can easily check that, for the case of free evolution, the dissipation terms become irrelevant in the long time scaling regime. The viscosity is chosen to ensure numerical stability while preserving the key features of the inviscid dynamics. Convergence is verified by simultaneously increasing the resolution and decreasing the viscosity [71]; in the asymptotic regime, we confirm that the features of the solution remain unaffected for the case of free evolution. Note however, that for the 1D trapped gas discussed in Sec. III, the viscosity plays a significant role in the long time relaxation to the steady state.

For one-dimensional simulations, we use $2^{12}\leq N\leq 2^{16}$ collocation points. In two dimensions, we use $N\times N$ grids with $2^{9}\leq N\leq 2^{11}$. The time step, $dt$, is chosen adaptively based on the system size, subject to the Courant–Friedrichs–Lewy (CFL) stability condition.

For freely evolving systems (i.e., without background density), we set the system size to $L=N$ ($dx=1$), and find convergence in the scaling regime for $L\sim 2^{14}$. For cases with a finite background density, convergence is slower, and we obtain reliable results for $L\sim 2^{16}$.

In the presence of a harmonic trapping potential, we solve Eq. (32) without adding artificial viscosity to the continuity equation. While the Fourier pseudospectral method remains applicable, the non-periodic nature of the potential is incompatible with strict periodic boundary conditions. To suppress boundary artifacts, we apply a mild smoothing of the potential near the edges of the domain. All simulations are performed with a fixed computational domain of length $L=2\pi$. This issue does not arise in the linearized perturbation equations [Eq (37)], where the quadratic trap does not appear explicitly and its effects are encoded through the steady-state background.

To solve the two-dimensional equations in polar coordinates [Eq. (54)], we impose non-periodic boundary conditions appropriate for the radial geometry. Specifically, we set $v=0$ and $\partial_{r}\rho=0$ at $r=0$, following Ref 72. In contrast to the Cartesian case, periodic boundary conditions hence pseudospectral methods are no longer applicable. We used a second-order finite difference scheme to approximate the spatial derivatives, while time integration was carried out using the classical fourth-order Runge–Kutta (RK4) method.

Most of the codes used in the numerical results reported here can be found in Ref 63.