Weakly interacting disordered Bose gases out of equilibrium: From multiple scattering to superfluidity^(a)

In quantum gases, the interplay between disorder and interactions leads to a rich variety of phenomena. At low temperature, for instance, a one-dimensional Bose gas in equilibrium behaves either as an insulator or a superfluid depending on the relative strength of disorder and interactions [1]. The insulating phase has been predicted to be robust at finite temperature, up to a critical point known as the many-body localization transition [2]. Initially described for interacting electrons in dirty conductors [3,4], many-body localization has, since then, triggered a considerable interest in various fields of physics [5]. Another intriguing problem is the dynamical evolution of interacting disordered gases brought out of equilibrium by an initial quench. In the homogeneous case, out-of-equilibrium interacting gases generically thermalize to a universal Gibbs ensemble [6]. For sufficiently strong disorder and interactions, on the other hand, many-body localization may prevent thermalization from occurring [5]. In between these two limits, the quench dynamics of disordered quantum gases remains largely unexplored, in particular in dimensions larger than one.

In this context, a yet poorly explored problem is the out-of-equilibrium dynamics of weakly interacting disordered gases. In this regime, mostly states in the ergodic phase of the many-body localization transition contribute to the dynamics, so that thermalization is the rule. For bosons, the weakly interacting regime can be tackled at the mean-field level using the disordered Gross-Pitaevskii equation. The dynamical establishment of thermalization in this context was previously studied in one-dimensional disordered chains [7,8] and in two-dimensional (2D) [9,10] random potentials, in cases where the disorder is typically larger than interactions. In this regime, the short-time dynamics of the gas is governed by multiple scattering, while interactions make the system thermalize on a longer time scale. The opposite limit of a disorder smaller than interactions has, on the other hand, little been explored in an out-of-equilibrium context.

In this letter, we study the quench dynamics of a 2D weakly interacting disordered Bose gas by varying the relative strengths of disorder and interactions over a wide range. We characterize this dynamics by numerically and analytically computing the momentum distribution and the coherence function of the gas, focusing on relatively short times, before the system is fully thermalized. For a disorder larger than interactions, we recover that the gas experiences multiple scattering, while being slowly thermalized by interactions. The spatial coherence is low in this limit. Upon increasing the strength of interactions, we find that scattering diminishes and that the gas behaves more and more like a fluid. Ultimately, when interactions become larger than the disorder, the gas becomes a non-equilibrium disordered superfluid with a long-range, algebraic coherence emerging in a light-cone fashion, which is typical of a 2D pre-thermalization process.

Following the approach of [9–15], we consider a simple out-of-equilibrium protocol where a 2D N-particle Bose gas, initially prepared in the plane-wave state $|\bm{k}_0\rangle$, undergoes at t = 0 a simultaneous interaction and disorder-potential quench. We describe the ensuing dynamics for t > 0 with the time-dependent disordered Gross-Pitaevskii equation

$$\begin{equation} i\partial_t\Psi(\bm{r},t)= \left[-\frac{\bm{\nabla}^2}{2m} + V(\bm{r})+gN|\Psi(\bm{r},t)|^2\right] \Psi(\bm{r},t)\quad \end{equation} \tag{ 1 }$$

for the Bose mean field $\Psi(\bm{r},t)$, normalized according to $\int \textrm{d}^2\bm{r} |\Psi(\bm{r},t)|^2=1$. Here and in the rest of the letter, we set $\hbar=1$. We choose the random potential $V(\bm{r})$ to be a blue-detuned speckle, with mean value $\overline{V(\bm{r})}=0$ and correlation function $\overline{V(\bm{r})V(\bm{r}')}=V_0^2\exp(-|\bm{r}-\bm{r}'|^2/2\sigma^2)$, where V₀ is the amplitude of disorder fluctuations, σ the correlation length, and the overbar refers to disorder averaging. In our numerical simulations, we generate this potential in a standard way, by convoluting a spatially δ-correlated complex Gaussian random field with a Gaussian cut-off function and taking the modulus square of the resulting field [16–17]. The temporal propagation of the initial plane wave is performed using a split-step algorithm that includes a Chebyshev expansion of the linear part of the evolution operator, as explained in [10]. Simulations are performed on a 2D regular grid of size $L\times L$ with periodic boundary conditions along x and y. A cell of surface $(\pi\sigma)^2$ is discretized in typically 5 to 7 steps along both x and y. Throughout the paper, numerical values of lengths, momenta, energies and times will be given in units of σ, $\sigma^{-1}$, $1/(m\sigma^2)$ and $m\sigma^2$, respectively. Finally, all the results are averaged over 3200 to 4800 disorder realizations.

We first compute numerically the disorder-averaged momentum distribution $n_{\bm{k}}(t)\equiv \overline{|\langle{\bm{k}}|\Psi(t)\rangle|^2}$ from eq. (1) using the procedure explained above, for various relative values of the disorder potential fluctuations, V₀, initial kinetic energy $\epsilon_0\equiv k_0^2/2m$ and initial interaction energy $g\rho_0\equiv gN/V$, where $V=L^2$ is the volume of the system. The distribution is normalized according to $\int \textrm{d}^2{\bm{k}}/(2\pi)^2n_{\bm{k}}(t)=1$. Density plots of $n_{\bm{k}}(t)$ are shown in fig. 1, with the columns corresponding to different parameter regimes and the rows to three successive times, chosen relatively short compared to the time where the gas gets completely thermalized by interactions.

• –  
  Regime $\epsilon_0>V_0> g\rho_0$ (figs. 1(a)): Particles of energy $\sim\epsilon_0$ experience elastic multiple scattering on the disorder potential, which randomizes the direction of their momenta. The distribution $n_{\bm{k}}(t)$ thus quickly becomes a ring of radius k₀. As time increases (from fig. 1(a-1) to (a-3)), the ring is slowly smoothed by particle collisions (thermalization).
• –  
  Regime $\epsilon_0> g\rho_0>V_0$ (figs. 1(b)): When the interaction strength is increased, particle collisions scramble the distribution before the scattering ring is fully formed.
• –  
  Regime $g\rho_0>\epsilon_0> V_0$ (figs. 1(c)): When interactions become larger than the initial kinetic energy, the gas starts to behave as a superfluid. There is no longer scattering, rather the disorder-induced field fluctuations are coherently enhanced and manifest themselves as interference rings, which become tighter and tighter as time increases.
• –  
  Regime $g\rho_0> V_0>\epsilon_0$ (figs. 1(d)): Same as in figs. 1(c), but for a Bose gas initially almost motionless, at rest on the plots $(\epsilon_0=0)$.

To better understand the physics at play in the distributions of fig. 1, we now discuss more in details the two extreme regimes of figs. 1(a) and (d).

Zoom In Zoom Out Reset image size

We first address the low-interaction regime of figs. 1(a). Since interactions are typically the smallest energy scale in this limit, the physics at short time is essentially a multiple-scattering process of quasi-particles with energy $\epsilon_k\equiv {\bm{k}}^2/2m$ in the random potential. The relevant time scale is here the scattering mean free time τ, which gives the rate at which these quasi-particles are elastically scattered. For the chosen Gaussian potential, we have

$$\begin{equation} \frac{1}{\tau}=2\pi V_0^2\sigma^2e^{-k_0^2\sigma^2}I_0(k_0^2\sigma^2), \end{equation} \tag{ 2 }$$

in the Born approximation [18], where I₀ is the modified Bessel function of the first kind. When t > τ, the disorder-averaged momentum distribution acquires a ring shape, as seen in figs. 1(a). For weak interactions, it was shown to be given by [10]

$$\begin{equation} n_{\bm{k}}(t)=\int \textrm{d}\epsilon A_\epsilon({\bm{k}})f_\epsilon(t)/\nu, \end{equation} \tag{ 3 }$$

where $f_\epsilon(t)$ is the energy distribution of the quasi-particles (normalized according to $\int \textrm{d}\epsilon f_\epsilon(t)=1$), $A_\epsilon({\bm{k}})$ is their spectral function and $\nu=m/2\pi$ is the 2D density of states per unit volume. In the weak-disorder regime where $k_0\ell\gg1$, we have $A_\epsilon({\bm{k}})\simeq1/(2\pi\tau)\times1/[(\epsilon-\epsilon_k)^2+1/4\tau^2]$. In the absence of interactions, $f_\epsilon(t)=A_\epsilon({\bm{k}}_0)$ is independent of time and coincides with the spectral function at ${\bm{k}}={\bm{k}}_0$. Equation (3) then reduces to [11]

$$\begin{equation} n_{\bm{k}}(g=0)=\frac{1}{\pi\nu\tau}\frac{1}{(\epsilon_k-\epsilon_0)^2+1/\tau^2}, \end{equation} \tag{ 4 }$$

which indeed describes a ring of half-width at half maximum $\Delta k\sim 1/\ell$, where $\ell=k_0\tau/m$ is the mean free path. When $g\ne 0$, the energy distribution $f_\epsilon(t)$ evolves in time due to particle collisions. At weak disorder, its dynamics is governed by a Boltzmann-like kinetic equation [9,10,19]. The particle collisions occur at a mean rate given by the Fermi golden rule

$$\begin{equation} \frac{1}{\tau_\text{col}}\!=\!2\pi\int\!\!\frac{\textrm{d}^2{\bm{k}}}{(2\pi)^2} \overline{|\langle{\bm{k}}_0|gN|\Psi|^2|{\bm{k}}\rangle|^2} \delta(\epsilon_0-\epsilon_{\bm{k}}) \!\sim\!\frac{(g\rho_0)^2}{\epsilon_0}.\quad \end{equation} \tag{ 5 }$$

While $f_\epsilon(t)$ ultimately becomes thermal when $t\gg\tau_\text{col}$ [9], as long as $t<\tau_\text{col}$ the scattering ring remains well resolved (as in figs. 1(a-1) to (a-3)) and we have found that the time evolution of $f_\epsilon(t)$ is approximately captured by a diffusion process in the energy domain:

$$\begin{equation} f_\epsilon(t)\simeq \int \textrm{d}\epsilon' \frac{A_{\epsilon'}({\bm{k}}_0)}{\sqrt{4\pi Dt}} \exp\left[-\frac{(\epsilon-\epsilon')^2}{4Dt}\right], \end{equation} \tag{ 6 }$$

with a diffusion coefficient in energy space such that $2D\tau_\text{col}\sim 1/\tau^2$. When $g\ne 0$, the momentum distribution at short time is thus obtained by convoluting the non-interacting result (4) with the Gaussian distribution in eq. (6). For $\tau\ll t< \tau_\text{col}$ this leads to a broadening of the ring with half-width $\Delta k(t)\sim (1+C t/\tau_\text{col})/\ell$, with C a numerical constant. This broadening is visible in the left panel of fig. 2, which shows radial cuts of the scattering ring, obtained numerically at successive times. As long as $t<\tau_\text{coll}$ we find that the widths extracted from these cuts indeed increase linearly in time, as shown in the right panel of fig. 2.

Zoom In Zoom Out Reset image size

It is also interesting to investigate the spatial coherence of the Bose gas in the multiple-scattering regime. To this aim, we consider the coherence function

$$\begin{equation} g_1(\bm{r},t)\equiv{V}\overline{\Psi^*(0,t)\Psi(\bm{r},t)}= \!\int\!\!\frac{\textrm{d}^2{\bm{k}}}{(2\pi)^2}n_{\bm{k}}(t)e^{i{\bm{k}}\cdot \bm{r}}. \end{equation} \tag{ 7 }$$

When g = 0, $n_{\bm{k}}$ is given by eq. (4), so that $g_1({\bm{r}},t)=J_0(k_0 r)\exp(-{r/\ell})$, with $r\equiv |{\bm{r}}|$ and J₀ the Bessel function of the first kind [18]. This expression is essentially the first-order correlation function of the 2D matter-wave speckle pattern formed by the particles scattered on the random potential [20]. It shows that the coherence length of the Bose gas is rather short in the multiple-scattering regime ¹ , of the order of the de Broglie wavelength $2\pi/k_0$. We show in fig. 3 the g₁ function computed numerically from eq. (1) for g = 0, together with the above theoretical formula. In the presence of interactions and as long as $t<\tau_\text{col}$, the shape of $n_{\bm{k}}$ remains close to the Lorentzian (4), with a broadened width $\Delta k(t)\sim(1+Ct/\tau_\text{col})/\ell$, such that

$$\begin{equation} g_1({\bm{r}},t)\simeq J_0(k_0r)\exp[-\Delta k(t)r]. \end{equation} \tag{ 8 }$$

As shown in fig. 3, eq. (8) describes very well the numerical results for g₁ at $g\ne 0$. Here the main effect of interactions is to smooth the matter-wave speckle, without significantly affecting the coherence of the gas. This picture would of course change at long time $t\gg \tau_\text{col}$, where the momentum distribution starts to significantly deviate from eq. (4) and approaches a thermal law [9].

Zoom In Zoom Out Reset image size

We now address the opposite, superfluid regime where $g\rho_0\gg V_0\gg \epsilon_0$. From now on, we focus for simplicity on the case $\epsilon_0=0$ of a gas initially at rest, corresponding to the numerical distributions in figs. 1(d-1) to (d-3). In such a low-energy regime, the natural approach for describing the dynamics of the 2D Bose gas relies on the density-phase representation of the wave function [21,22]

$$\begin{equation} \Psi({\bm{r}},t)=\sqrt{\rho({\bm{r}},t)}\exp[i\theta({\bm{r}},t)-ig\rho_0 t]{/\sqrt{N}}, \end{equation} \tag{ 9 }$$

with the initial condition $\rho({\bm{r}},0)=\rho_0$ and $\theta({\bm{r}},0)=0$. $-ig\rho_0 t$ is the dynamical phase of the uniform solution in the absence of disorder (${\rho({\bm{r}},t)}\equiv\rho_0$) and the $1/\sqrt{N}\equiv 1/\sqrt{\rho_0 V}$ factor guarantees our choice of normalization at t = 0, $\int \textrm{d}^2{\bm{r}} |\Psi({\bm{r}},t)|^2=1$. As the gas evolves in time, the density $\rho({\bm{r}},t)$ and the phase $\theta({\bm{r}},t)$ fluctuate due to the presence of the disorder potential. We then exploit that the disorder potential is the smallest energy scale to make use of perturbation theory. Precisely, we write $\rho({\bm{r}},t)=\rho_0+\delta \rho({\bm{r}},t)$ and linearize the Gross-Pitaevskii equation (1) with respect to $\delta \rho({\bm{r}},t)$ and $\nabla\theta({\bm{r}},t)$ (without any assumption on the phase itself, which may be strongly fluctuating in this regime) [23]. This leads to the Bogoliubov-de Gennes equations

$$\begin{equation} \frac{\partial \delta\rho}{\partial t} = -\frac{\rho_0}{m}\nabla^2\theta, \end{equation} \tag{ 10 }$$

$$\begin{equation} \frac{\partial \theta}{\partial t} = \frac{1}{4m\rho_0}\nabla^2\delta\rho-V({\bm{r}})-g\delta\rho, \end{equation} \tag{ 11 }$$

which can be readily diagonalized to provide the density,

$$\begin{eqnarray} &&\delta\rho({\bm{r}},t)=-2\rho_0\int\frac{\textrm{d}^2{\bm{k}}}{(2\pi)^2}{V}({\bm{k}})\frac{1-\cos(E_k t)}{\epsilon_k+2g\rho_0}e^{i{\bm{k}}\cdot{\bm{r}}}\qquad \end{eqnarray} \tag{ 12 }$$

and the phase fluctuations,

$$\begin{eqnarray} &&\theta({\bm{r}},t)=-\int\frac{\textrm{d}^2{\bm{k}}}{(2\pi)^2}{V}({\bm{k}})\frac{\sin(E_k t)}{E_k}e^{i{\bm{k}}\cdot{\bm{r}}}, \end{eqnarray} \tag{ 13 }$$

where $V({\bm{k}})$ refers to the Fourier transform of $V({\bm{r}})$. These expressions involve the well-known Bogoliubov dispersion relation, $E_k\equiv\sqrt{\epsilon_k(\epsilon_k+2g\rho_0)}$, where we recall that $\epsilon_k={\bm{k}}^2/2m$. The coherence function follows from the relation

$$\begin{eqnarray} g_1({\bm{r}},t)& = \exp\left\{-\int\frac{\textrm{d}^2{\bm{k}}}{(2\pi)^2}(1-\cos{\bm{k}}\cdot{\bm{r}})\right.\nonumber\\ & \quad\times\left[\overline{|\theta({\bm{k}},t)|^2}+\overline{|\delta\rho({\bm{k}},t)|^2}/(4\rho_0^2)\right] \bigg\},\quad \end{eqnarray} \tag{ 14 }$$

deduced from the Bogoliubov-Popov theory of fluctuations [22–24], here adapted to the classical disorder-induced fluctuations $\theta({\bm{r}},t)$ and $\delta\rho({\bm{r}},t)$. Inserting eqs. (12) and (13) into eq. (14) and performing the disorder averages, we find:

$$\begin{eqnarray} g_1({\bm{r}},t)&= \exp\left\{-\int\frac{\textrm{d}^2{\bm{k}}}{(2\pi)^2}B({\bm{k}})(1-\cos{\bm{k}}\cdot{\bm{r}})\right.\nonumber\\ & \quad\left. \times\left(\left[\frac{1-\cos(E_k t)}{\epsilon_k+2g\rho_0} \right]^2+ \left[\frac{\sin(E_k t)}{E_k} \right]^2\right) \right\},\nonumber\\ \end{eqnarray} \tag{ 15 }$$

where $B({\bm{k}})\!\equiv\!\int \textrm{d}{\bm{r}} e^{i{\bm{k}}\cdot {\bm{r}}}\overline{V(0)V({\bm{r}})}\!=\!2\pi V_0^2\sigma^2e^{-q^2\sigma^2/2}$ is the disorder power spectrum. The physical content of eq. (15) is the dynamical spreading of correlations in a Bose gas initially quenched from a weakly fluctuating state [25], here in two dimensions. Unlike quench configurations more commonly encountered in the literature however [24–29], in our case the fluctuations are neither of quantum or thermal origin, but come from the random potential. In the present regime, scattering does not take place, so that the mean free time τ is no longer a relevant time scale. Instead, the short-time dynamics of the Bose gas is governed by the nonlinear time $\tau_\text{NL}\equiv 1/g\rho_0$. For $t\gg\tau_\text{NL}$ and at the intermediate separations $\xi\ll r\equiv|{\bm{r}}|\ll 2c_s t$, where $\xi\equiv 1/\sqrt{g\rho_0m}$ is the healing length and $c_s=\sqrt{g\rho_0/m}$ is the speed of sound, eq. (15) simplifies to the time-independent algebraic law

$$\begin{equation} g_1({\bm{r}},t)\simeq \biggl[G(\sigma/\xi)\frac{\xi}{r}\biggr]^{\alpha}, \end{equation} \tag{ 16 }$$

where $\alpha\equiv(V_0\sigma/\sqrt{2}g\rho_0\xi)^2$ and $G(x)=\sqrt{2}x\exp\{[(6x^2+1)\exp(2x^2)E_1(2x^2)-3-\gamma]/2\}$, with γ being the Euler-Mascheroni constant and E₁ the first-order exponential integral. A similar 2D algebraic scaling was found in [30], in a configuration where the fluctuations were however encoded in the initial state and not in a disorder potential. Note that eq. (16) is both independent of time and of the precise shape of the disorder spectrum $B({\bm{k}})$. It only depends on the small set of parameters $\{\xi,\sigma,g\rho_0, V_0\}$ and is thus, to a large extent, universal. Time independence and universality are characteristic properties of a pre-thermalization process, where a quenched system quickly converges to a fixed, thermal-like point, from where it only departs very slowly [26–29]. In the present scenario, the existence of a pre-thermalization regime requires the disorder amplitude to be very weak, so that the system is nearly integrable at short time. In the pre-thermal regime, the correlation function (16) looks like the one of a uniform 2D Bose superfluid at equilibrium and finite temperature [31], and the gas exhibits an algebraically decaying coherence [22], in strong contrast with the multiple-scattering regime (compare with eq. (8)). Equation (16) breaks down at $r=2c_s t$, the boundary of a light cone within which correlations can spread. Out of the light cone, $g_1({\bm{r}},t)\simeq[G(\sigma/\xi)\xi/4c_st]^\alpha$ reaches a time-dependent plateau, reminiscent of the perfect coherence of the initial plane-wave state. Let us mention that to observe the algebraic law (16), it is required to use of a correlated potential with $\sigma\leq \xi$ so to select the low-k phonon modes $E_k\simeq c_s k$ in eq. (15). In particular, for an uncorrelated potential we have found no numerical evidence of algebraic decay in the coherence function.

Within the linearization approach that led us to eq. (15), the disorder fluctuations are described in terms of independent Bogoliubov quasi-particles, which reduce to phonons at low energy. Collisions between these phonons make the system slowly depart from the pre-thermal regime described by eq. (16) [32]. While the effect of these collisions remains small at short time, it was found in [30] that even at short time it is necessary to take them into account in order to achieve a quantitative agreement with numerical simulations for g₁. The consequence of phonon collisions is a deviation of the fluid fluctuations $\overline{|\theta({\bm{k}},t)|^2}$ and $\overline{|\delta\rho({\bm{k}},t)|^2}$ in eq. (14) from their quadratic expressions [32]. At short time, this deviation is small and thus reduces to a correction term in the second line of eq. (15). We take it of the form $\beta(t)/k$, i.e., replace the second line of eq. (15) by $[(1\!-\!\cos(E_k t))/(\epsilon_k\!+\!2g\rho_0)]^2+[\sin(E_k t)/E_k]^2+\beta(t)/k$, and use $\beta(t)$ as a fitting parameter. The factor $1/k$ in the correction term is introduced somewhat arbitrarily to reduce the weight of phonon collisions at short scale. We have verified that an alternative fitting option, independent of k, also allows us to reproduce the numerical results [30], albeit less accurately at short scale. A comparison between eq. (15), modified according to this procedure, and the exact simulations for g₁, is shown in fig. 4(parameters are here the same as those in figs. 1(d)). The agreement is excellent, except for the small spatial variations of g₁ in the vicinity of the light-cone boundary, which are not present in the simulation results. One reason might be a smoothing of these variations due to particle collisions.

Zoom In Zoom Out Reset image size

When $V_0/g\rho_0\ll1$, one can expand the exponential in eq. (15), so that the momentum distribution $n_{\bm{k}}(t)$ at ${{\bm{k}}\ne0}$ is in first approximation given by the simple law

$$\begin{eqnarray} &&n_{\bm{k}}(t)\simeq B({\bm{k}}) \left(\left[\frac{1\!-\!\cos(E_k t)}{\epsilon_k+2g\rho_0} \right]^2\!+\! \left[\frac{\sin(E_k t)}{E_k} \right]^2\right).\qquad \end{eqnarray} \tag{ 17 }$$

Such a profile consists of concentric rings whose minima are located at positions $k_n=(\sqrt{2}/\xi)\{-1+[1+(\pi n/g\rho_0 t)^2]^{1/2}\}^{1/2}$, where n is a non-zero positive integer. These rings are well visible in the distributions in figs. 1(d-1) to (d-3). Physically, they originate from both the interference between phonons (interference pattern $\sim\cos(2E_kt)$) and between the phonons and the mean field ($\sim\cos(E_kt)$). The spacing $\sim \pi/(c_s t)$ between the rings decreases in time, signaling that the interfering particles are further and further apart as time grows. Equation (17) is displayed in the inset of fig. 4, together with cuts along k_x of the distributions in figs. 1(d-1) to (d-3) extracted from simulations. The agreement is again very good, in particular for the positions of the minima, indicated by vertical lines.

We have theoretically described the dynamics of a 2D, weakly interacting Bose gas in the presence of a random potential, and have identified two qualitatively different non-equilibrium regimes depending on the relative strengths of disorder and interactions. For interactions weaker than the disorder, the physics is that of multiple scattering, with interactions slowly thermalizing the energy distribution and a short-range coherence at short time. When interactions are stronger than the disorder, on the other hand, we recover a low-energy-physics regime: scattering ceases, the gas becomes a non-equilibrium superfluid, with long-range correlations spreading within a light cone.

NC thanks Giovanni Martone for discussions. NC and DD acknowledge financial support from the Agence Nationale de la Recherche (grants Nos. ANR-19-CE30-0028-01 CONFOCAL and ANR-18-CE30-0017 MANYLOK).