Cloud shape of a molecular Bose–Einstein condensate in a disordered trap: a case study of the dirty boson problem

Disorder is ubiquitous in nature, since there is no perfect order in any realistic physical system. Random disordered media have for long time been studied in terms of percolation and transport phenomena [1–3] or anomalous diffusion [4]. The dirty boson problem first arose experimentally in the context of superfluid helium in porous vycor glass, where it was shown that superfluidity can still persist despite the influence of disorder [5, 6]. It amounts to understanding the emergence of coherence and order for ultracold bosonic atoms in the presence of frozen disorder in the form of a static random potential. Its importance and intriguing aspects stem from the interplay of localization and interactions as well as of disorder and superfluidity. Since the realization of Bose–Einstein condensates (BECs) in 1995 [7, 8] the interest regarding the effects of disorder potentials on the properties of ultracold quantum gases has substantially increased [9]. This was initiated by unavoidable irregularities in the trapping potential induced by wire imperfections [10, 11]. However, ultracold quantum gases became controllable to an unprecedented level of precision. For example, disorder can nowadays be imposed thereon in a highly tunable manner via laser speckles [12, 13], light fields engineered by digital micromirror devices [14, 15] or atomic impurities trapped in an optical lattice [16, 17]. Hence, they represent an extremely promising platform for practically realizing Feynman's quantum simulator [18], which mimics the physics of other quantum many-body systems [19, 20]. This especially applies to many phenomena of condensed matter physics, which are inevitably influenced, or even caused, by various types of uncontrollable disorder. For instance, the phenomenon of Anderson localization, originally used to microscopically describe the absence of diffusion in the presence of disorder [21], was predicted for a flow of a dilute Bose gas through a disordered region [22]. So far it was observed for non-interacting atoms in one dimension [23, 24] as well as presumably in three dimensions [25, 26]. More recently, the so called many-body localization has been observed for interacting particles in disordered lattice systems and has attracted significant interest [27–29]. In this paper, however, we are more interested in analyzing the impact of disorder upon bosonic many-body quantum systems in the bulk. Therefore, we briefly review the literature on the homogeneous or trapped dirty boson problem both for a weak and a strong random potential.

In order to quantitatively study a weakly interacting homogeneous Bose gas in a static random potential, Huang and Meng worked out a Bogoliubov theory [30], where quantum and thermal fluctuations as well as disorder fluctuations were assumed to be small. This approach treated the disorder perturbatively and was subsequently extended by others, either within the original framework of second quantization [31–34] or within the functional integral approach including the replica method [35–37]. In particular, it was also demonstrated that, despite the randomness of the potential, superfluidity persists and that its depletion is larger than the condensate depletion since the bosons scattered by the disorder landscape represent randomly distributed obstacles for the motion of the superfluid. An extension to the situation when the disorder correlation function falls off with a characteristic correlation length, as in the case of a Gaussian [32, 38–40], a Lorentzian [41], or laser speckle disorder [42, 43] revealed that condensate and superfluid depletions decrease with increasing correlation length. Note that, to the best of our knowledge, so far no prediction of the Huang–Meng theory has yet been checked experimentally although it stems from 1992. Finally, some subtleties regarding condensate deformation and its quantum depletion in external, and possibly random, weak potentials were resolved in [44, 45] using an inhomogeneous Bogoliubov theory.

The dirty boson problem was also studied non-perturbatively using various approaches. For a delta-correlated disorder, either in the homogeneous case [46–51] or in a disordered trap [52–54], a paradigmatic result is that an increase of the disorder strength at zero temperature leads to a first-order quantum phase transition from a superfluid to a Bose-glass phase. In the latter, global phase coherence is destroyed and phase coherence only persists locally due to the particle scattering with the disorder. On the other hand, non-perturbative studies of dirty bosons with a finite disorder correlation length are scarce. For instance, the diffusion Monte-Carlo study in [55], where the disorder is implemented within a hard-sphere model, concludes for the homogeneous case that no quantum phase transition from a superfluid to a Bose-glass occurs. This result is confirmed by the theoretical investigations of this paper, where the case of a disordered trap is studied for concrete experimental parameters.

Note that, until now, it has remained a challenge to test all these theoretical predictions within a concrete experimental setup. This motivated us to investigate within an experiment-theory collaboration the static cloud shape of a BEC in a disordered trap and its change with increasing disorder strength. In section 2 we introduce the underlying experimental setup, where the disorder is realized by laser speckles and where the harmonically trapped BEC consists of molecules of fermionic ⁶Li atoms. Subsequent section 3 develops a corresponding theoretical description, which is non-perturbative with respect to disorder and includes quantum fluctuations. Finally, section 4 presents the measured cloud widths as a function of the laser speckle strength and discusses how they can be explained both qualitatively and quantitatively by the developed theory. Although we have achieved a quite significant level of agreement between the measurements and the theory, we also note that some experimental results indicate certain limitations of our theoretical approach, which needs to be refined for future studies.

In the following, relevant components of the experimental setup are presented. First, we give an overview of the sequence for the creation of molecular BECs. Subsequently we focus on statistical properties of the optical speckle potential, its interaction with lithium molecules and a novel method for calibrating the disorder strength. Finally, we describe how to load the quantum gas into a disordered trap, the measurement protocol as well as the data analysis.

2.1. Creation of molecular Bose–Einstein condensates

A general overview of our experimental setup is given in [56]. In short, we prepare quantum degenerate gases of neutral, fermionic ⁶Li atoms in an equal mixture of the two lowest-lying Zeeman substates $| {m}_{J}=-1/2,{m}_{I}=1\rangle$ and $| -1/2,0\rangle$ of the electronic ground state ${}^{2}{{\rm{S}}}_{1/2}$, where m_J and m_I are the magnetic quantum numbers of the electronic and nuclear spin, respectively. A magnetic Feshbach resonance centered at $832\,{\rm{G}}$ [57] is utilized to tune the s-wave scattering length between atoms of opposite spin to $4510\,{{a}}_{0}$ at $763.6\,{\rm{G}}$, where a₀ denotes the Bohr radius. At this magnetic field, the interatomic interaction potential features a molecular state with binding energy ${E}_{{\rm{b}}}/{k}_{{\rm{B}}}=1.5\,\mu {\rm{K}}$ [58], where ${k}_{{\rm{B}}}$ denotes the Boltzmann constant. Such Feshbach molecules form during cooling once the sample temperature approaches ${E}_{{\rm{b}}}/{k}_{{\rm{B}}}$ [59] and they interact with scattering length $2706\,{{a}}_{0}$ [60]. Molecular BECs of typically (4.3 ± 0.2) × 10⁵ molecules are created by forced evaporative cooling in a superposition of optical and magnetic potentials. An optical dipole trap confines the cloud radially, which corresponds to the x- and z-direction in figure 1, with trap frequencies $({\omega }_{x},{\omega }_{z})=2\pi \times (195\,\mathrm{Hz},129\,\mathrm{Hz})$. The weaker axial confinement, along the y-direction in figure 1, is provided by a magnetic saddle potential for which we have ${\omega }_{y}=2\pi \times 22.6\,\mathrm{Hz}$. After evaporation, the atomic sample is held for $500\,\mathrm{ms}$ at a constant trap depth of $500\,\mathrm{nK}\times {k}_{{\rm{B}}}$ to ensure thermal equilibration, before studying its properties in a disordered potential.

Zoom In Zoom Out Reset image size

2.2. Optical speckle potential

The optical speckle potential is realized by passing a Gaussian laser beam through a holographic diffusive plate (Edmund Optics 47-991) and imaging the light onto the position of the atoms. The beam is derived from a frequency-doubled Nd:YVO solid-state laser (Coherent Verdi V18) with wavelength $\lambda =532\,\mathrm{nm}$. We control its power irradiated onto the atomic sample by means of an acousto-optic modulator in a proportional-integral-derivative control loop. After passing through the diffuser, the beam waist is enlarged to $w=29\,\mathrm{mm}$ by a telescope to match the free aperture of the imaging objective with numerical aperture $0.29$ and focal length $f=98.8\,\mathrm{mm}$.

The speckle grain size or correlation length is characterized by the width of the intensity autocorrelation function $C(\delta {\bf{r}})=\overline{I({\bf{r}})I({\bf{r}}+\delta {\bf{r}})}$ [12], where $I({\bf{r}})$ is the intensity at the point ${\bf{r}}$ and the bar denotes averaging over a region that encompasses a large number of speckle grains. For a Gaussian beam with waist w, the autocorrelation function in the image plane takes on the shape of a Gaussian

$$\begin{eqnarray}&&C(\delta {\bf{r}})=\overline{{I}^{2}}\left[1+\exp \left(-\displaystyle \frac{\delta {{\bf{r}}}^{2}}{{\sigma }^{2}}\right)\right],\end{eqnarray} \tag{ 1 }$$

with the diffraction-limited width $\sigma =\lambda f/(\pi w)=576\,\mathrm{nm}$ for $w=29\,\mathrm{mm}$. Along the beam propagation axis the grain size is expected to be larger by a factor of 4f/w = 13.6 [61]. A speckle is called fully developed if the phases of contributing partial waves are distributed uniformly across the interval [0, 2π) [12]. In this case, dark regions are more abundant than speckle grains of any higher intensity [62] and one can show that the probability density of the intensity follows an exponential law [12]

$$\begin{eqnarray}&&P(I)=\displaystyle \frac{1}{\bar{I}}\exp \left(-\displaystyle \frac{I}{\bar{I}}\right).\end{eqnarray} \tag{ 2 }$$

Both the intensity autocorrelation and probability distribution are measured in a test setup featuring an optical beam path for projection of the laser speckle that is identical to the one used in the experiment. It includes a microscope at the intended atomic cloud position which images the speckle onto a CCD camera. The measured intensity probability distribution, shown in figure 2, is exponential. Thus, we conclude that the speckle is, indeed, fully developed. The measured correlation length ${\sigma }_{{\rm{m}}}=750\,{\rm{n}}{\rm{m}}$ is larger than theoretically expected, which is in part due to a truncation of the Gaussian beam by the finite aperture of the objective. Taking this into account, a numerical simulation of the speckle provides a correlation length of $650\,\mathrm{nm}$. The remaining deviation is attributed to potential aberrations of our optical system.

Zoom In Zoom Out Reset image size

Since the atomic transition of lithium, which is affected the most by the laser speckle, is the far detuned D2 line at $671\,\mathrm{nm}$, the interaction with the light field is dominated by the optical dipole force [63]. The corresponding potential reads

$$\begin{eqnarray}&&u({\bf{r}})=-\displaystyle \frac{3\pi {c}^{2}{\rm{\Gamma }}}{2{\omega }_{0}^{3}}\left(\displaystyle \frac{1}{{\omega }_{0}-\omega }+\displaystyle \frac{1}{{\omega }_{0}+\omega }\right)I({\bf{r}}),\end{eqnarray} \tag{ 3 }$$

with the speed of light c, the natural line width ${\rm{\Gamma }}=2\pi \times 5.87\,\mathrm{MHz}$ [64], the speckle laser frequency ω = 2πc/λ, and the atomic transition frequency ω₀. As ω > ω₀, the potential is repulsive and we define the disorder strength as the mean potential $\overline{u}$ induced by the mean intensity $\bar{I}$. Spontaneous scattering is negligible (<10⁻³ s⁻¹) for all intensities used in our experiments. The polarisability of ⁶Li₂ molecules is twice that of unbound atoms [65] and, therefore, the molecules experience double the dipole potential. The maximal disorder strength for molecules we can achieve with the available laser power is $\bar{u}/{k}_{{\rm{B}}}=356\,\mathrm{nK}$.

Another relevant characteristic of optical speckle is the spatial dependence of its mean intensity $\bar{I}$. Due to the finite scattering angle ϑ of the diffusive plate, the light field in the image plane extends over a range far greater than the diffraction-limited spot size and its envelope is described by a Gaussian beam shape with waist ${w}_{{\rm{s}}}=\vartheta f$ [12]. In our setup, we have $\vartheta =0.5^\circ$ and therefore we obtain ${w}_{{\rm{s}}}=862\,\mu {\rm{m}}$, which is larger than the typical size $300\,\mu {\rm{m}}$ of an ultracold atomic cloud. For a precise calibration of the mean intensity it is therefore crucial not only to know the speckle waist but also the cloud's position therein. Other calibration methods are restricted to attractive potentials [66] or rely on measurements of the momentum distribution [67] and are, thus, not applicable to our system.

We have developed a novel method, which enables us to align the speckle center to the atomic cloud and infer the waist simultaneously. It is based on rotating the speckle pattern as a whole and measuring the resulting displacement of an atomic cloud within the speckle, see figure 3. Our diffusive plate is attached to a motorized rotary stage (Thorlabs K10CR1/M). Spinning this stage causes the emerging speckle pattern to revolve around its principal axis in the image plane. For the alignment procedure, we create a molecular BEC and linearly ramp the speckle laser power to its maximum value during $50\,\mathrm{ms}$. The sample is held for $100\,\mathrm{ms}$ in the potential rotating with angular velocity ω = 20° s⁻¹ before being released from the trap and imaged along the speckle principal axis after $11\,\mathrm{ms}$ time of flight. Since the disorder strength around the envelope center at x = y = 0 is large enough to drag the cloud along with its motion, the displacement of the atoms Δs is approximately the distance traversed by the speckle ${\rm{\Delta }}d(x,y)=\omega \sqrt{{x}^{2}+{y}^{2}}{t}_{{\rm{i}}}$ during the illumination time t_i. The value ${t}_{{\rm{i}}}=136\,\mathrm{ms}$ takes into account the weighted ramp duration and the time of flight. For increasing distance from the center, the decrease in mean intensity will eventually lead to a decline of Δs, although Δd continues to grow. We model this situation by assuming Δs to be, up to a constant factor, given by the product of Δd and the envelope function:

$$\begin{eqnarray}&&{\rm{\Delta }}s(x,y)\propto {\rm{\Delta }}d(x,y)\,\displaystyle \frac{\bar{I}(x,y)}{\bar{I}(0,0)}={\rm{\Delta }}d(x,y)\exp \left(-2\displaystyle \frac{{x}^{2}+{y}^{2}}{{w}_{{\rm{s}}}^{2}}\right).\end{eqnarray} \tag{ 4 }$$

In order to infer w_s, we scan the relative position of the quantum gas with respect to the speckle potential across the x- or y-axis. For each position, we measure Δy or Δx for clockwise ($\circlearrowright$) and anticlockwise ($\circlearrowleft$) rotation and evaluate ${\rm{\Delta }}{y}_{\mathrm{diff}}={\rm{\Delta }}{y}_{\circlearrowright }-{\rm{\Delta }}{y}_{\circlearrowleft }$, where ${\rm{\Delta }}{y}_{\circlearrowright }\approx -{\rm{\Delta }}{y}_{\circlearrowleft }$, or ${\rm{\Delta }}{x}_{\mathrm{diff}}$ in order to increase the signal strength. The result of such a measurement is shown in figures 3(b) and (c), where each data point is the average over 10 repetitions for each rotational direction. From fits according to (4) along the respective axes we extract ${w}_{{\rm{s}},x}=808(38)\,\mu {\rm{m}}$ and ${w}_{{\rm{s}},y}=884(63)\,\mu {\rm{m}}$, which agrees with ${w}_{{\rm{s}}}=862\,\mu {\rm{m}}$ within less than 10% deviation. Note that the smaller displacement along y compared to x is due to an additional confining potential that is not turned off during time of flight. As a possible cause for the ellipticity we identify a slight astigmatism induced by the acousto-optic modulator. For further experiments, we shift the peak of the speckle envelope to the atomic cloud position in order to achieve an almost homogeneous random potential with about 6% variation of $\bar{I}$ across the typical cloud size.

Zoom In Zoom Out Reset image size

2.3. Molecular Bose–Einstein condensates in a disordered trap

In order to load the molecular BEC into the disordered potential, the speckle potential is linearly ramped to its final value $\overline{u}$/${k}_{{\rm{B}}}$ within the interval $[0,356]\,{\rm{n}}{\rm{K}}$ during $50\,\mathrm{ms}$. After an additional hold time of $100\,\mathrm{ms}$ we probe the system by performing resonant high-intensity absorption imaging [68] on the transition $| {m}_{J}=-1/2,{m}_{I}=1\rangle \leftrightarrow {| }^{2}{{P}}_{3/2},{m}_{J}=-3/2,{m}_{I}=1\rangle$ along the z-axis and, therefore, measuring the in situ atomic column density in the x–y plane. We repeat this measurement about 140 times for each sampled disorder strength. Before every iteration, the microscopic details of the random potential experienced by the cloud are altered by rotating the speckle plate. This procedure enables us to average observables over many disorder realizations and, therefore, conform with the theoretical model as described in the subsequent section 3. From a semi-classical fit [69] to integrated density profiles without disorder we estimate the temperature to be $70\,\mathrm{nK}$, which amounts to about 25% of the non-interacting critical temperature for the onset of trapped molecular BEC. As this fit does not account for the quantum depletion of the BEC, the determined temperature has to be considered as an upper bound.

Figure 4 depicts exemplary density profiles for different disorder strengths. The impact of increasing disorder manifests itself in the form of emerging density fluctuations and an alteration of the cloud shape, which is characterized by the spatial extension along the principal axes. Since the healing length $\xi =270\,\mathrm{nm}$ of the BEC in the trap center is smaller than all correlation lengths, the wave function resolves the finest structures of the speckle potential [70]. These are, however, not visible in the recorded density distributions due to the limited resolution of our imaging setup, which amounts to $2.2\,\mu {\rm{m}}$. In order to extract the widths and simultaneously reduce the noise introduced by the fluctuations, we integrate the 2D column density profiles from figure 4(a) along the y-direction and fit a Gaussian function from which we get the 1/e half-widths a_x. Analogously, we obtain a_y from Gaussian fits of the density profiles integrated along the x-direction. Although the interaction strength is quite large, this ansatz is well suited to describe the integrated profiles even for the largest explored disorder strength, as illustrated in figures 4(b) and (c).

Zoom In Zoom Out Reset image size

Let us consider a weakly interacting Bose gas at T = 0 in the ground state in an external trapping potential $U({\bf{r}})$ with the addition of a disorder potential $u({\bf{r}})$. We assume that the disorder is homogeneous after the ensemble average, i.e. $\langle u({\bf{r}})\rangle =\langle u\rangle$. Otherwise, any non-homogeneity of the average can be absorbed into $U({\bf{r}})$. Since the correlation length is much smaller than the spatial extension of the cloud in x- and y-direction, we conclude that $\langle u\rangle$ matches the spatial average $\overline{u}$. The maximal value of the gas parameter in our experimental setup is $n({\bf{0}}){a}_{s}^{3}=0.017$, where $n({\bf{0}})$ denotes the clean-case particle density in the trap center, see (6) below. Comparing with the typical gas parameters for weakly interacting ⁸⁷Rb gas of about 0.0007 [71] and for a strongly interacting ⁴He of about 0.24 [72], we conclude that our experiment is in an intermediate regime, where quantum fluctuations have to be accounted for. Hence, the ground state wave function $\psi ({\bf{r}})$ satisfies the extended Gross–Pitaevskii equation that includes the quantum-fluctuation corrections of Lee et al [73] on the mean-field level [74]

$$\begin{eqnarray}&&\hspace{-22mm}-\displaystyle \frac{{{\hslash }}^{2}}{2m}{{\rm{\nabla }}}^{2}\psi ({\bf{r}})+U({\bf{r}})\psi ({\bf{r}})+u({\bf{r}})\psi ({\bf{r}})-\mu \psi ({\bf{r}})+g\left[\psi {\left({\bf{r}}\right)}^{3}+\displaystyle \frac{40}{3\sqrt{\pi }}{a}_{s}^{3/2}\psi {\left({\bf{r}}\right)}^{4}\right]=0,\end{eqnarray} \tag{ 5 }$$

where a_s is the s-wave scattering length and $g=4\pi {{\hslash }}^{2}{a}_{s}^{}/m$ denotes the interaction strength. The ground state chemical potential μ is taken to be independent of the disorder realization. This leads to a distribution of the total number of particles for different realizations of the disorder. The value of μ is fixed by the requirement that on average one has

$$\begin{eqnarray}&&N={\int }_{{{\mathbb{R}}}^{3}}{{\rm{d}}}^{3}{\bf{r}}\,\langle n({\bf{r}})\rangle ,\quad \mathrm{with}\quad n({\bf{r}})=\psi {\left({\bf{r}}\right)}^{2}+\displaystyle \frac{8}{3\sqrt{\pi }}{a}_{s}^{3/2}\psi {\left({\bf{r}}\right)}^{3},\end{eqnarray} \tag{ 6 }$$

where the second term of the particle density $n({\bf{r}})$ accounts for the quantum depletion of the condensate [74].

Averaging the condensate wave function over many disorder realizations gives $\langle \psi ({\bf{r}})\rangle$ that satisfies the average of (5)

$$\begin{eqnarray}&&\begin{array}{l}-\displaystyle \frac{{{\hslash }}^{2}}{2m}{{\rm{\nabla }}}^{2}\langle \psi ({\bf{r}})\rangle +U({\bf{r}})\langle \psi ({\bf{r}})\rangle +\langle u({\bf{r}})\psi ({\bf{r}})\rangle -\mu \langle \psi ({\bf{r}})\rangle \\ \qquad \qquad \quad +g\left[\langle \psi {\left({\bf{r}}\right)}^{3}\rangle +\displaystyle \frac{40}{3\sqrt{\pi }}{a}_{s}^{3/2}\langle \psi {\left({\bf{r}}\right)}^{4}\rangle \right]=0.\end{array}\end{eqnarray} \tag{ 7 }$$

Let us define $\varphi ({\bf{r}})=\psi ({\bf{r}})-\langle \psi ({\bf{r}})\rangle$ to be the part of $\psi ({\bf{r}})$ that fluctuates due to disorder. It obeys $\langle \varphi ({\bf{r}})\rangle =0$, $\langle \varphi {\left({\bf{r}}\right)}^{2}\rangle =\langle \psi {\left({\bf{r}}\right)}^{2}\rangle -\langle \psi ({\bf{r}}){\rangle }^{2}$ and satisfies

$$\begin{eqnarray}&&\begin{array}{l}-\displaystyle \frac{{{\hslash }}^{2}}{2m}{{\rm{\nabla }}}^{2}\varphi ({\bf{r}})+U({\bf{r}})\varphi ({\bf{r}})+(u({\bf{r}})\psi ({\bf{r}})-\langle u({\bf{r}})\psi ({\bf{r}})\rangle )-\mu \varphi ({\bf{r}})\\ \qquad \qquad \quad +\,g\left[\psi {\left({\bf{r}}\right)}^{3}-\langle \psi {\left({\bf{r}}\right)}^{3}\rangle +\displaystyle \frac{40}{3\sqrt{\pi }}{a}_{s}^{3/2}\left(\psi {\left({\bf{r}}\right)}^{4}-\langle \psi {\left({\bf{r}}\right)}^{4}\rangle \right)\right]=0.\end{array}\end{eqnarray} \tag{ 8 }$$

From the previous two equations one observes that the averages and the fluctuations are mutually coupled. Moreover, $\langle \psi ({\bf{r}})\rangle$ depends on $\langle u({\bf{r}})\psi ({\bf{r}})\rangle$ and higher-order averages like $\langle \psi {\left({\bf{r}}\right)}^{3}\rangle$ and $\langle \psi {\left({\bf{r}}\right)}^{4}\rangle$. The latter signals the necessity to deal with an infinite hierarchy of averages of increasing order. In order to use analytical tools, we will introduce several simplifying approximations.

First, we are going to use the cumulant expansion method [75–77], where in the Gaussian approximation all the cumulants of higher than second order are neglected. The general expansion method and the approximation procedure are outlined in appendix A. In this manner we have, for instance, $\langle \psi {\left({\bf{r}}\right)}^{3}\rangle \approx 3\langle \psi ({\bf{r}})\rangle \langle \psi {\left({\bf{r}}\right)}^{2}\rangle \,-2\langle \psi ({\bf{r}}){\rangle }^{3}$ and $\langle \psi {\left({\bf{r}}\right)}^{4}\rangle \approx 3\langle \psi {\left({\bf{r}}\right)}^{2}{\rangle }^{2}-2\langle \psi ({\bf{r}}){\rangle }^{4}$. Such a procedure truncates the hierarchy of the disorder averages of increasing order and enables an approximate but essentially non-perturbative treatment of the disorder up to the second cumulant order. Thus, applying the cumulant expansion to (7) leads to

$$\begin{eqnarray}&&\begin{array}{l}-\displaystyle \frac{{{\hslash }}^{2}}{2m}{{\rm{\nabla }}}^{2}\langle \psi ({\bf{r}})\rangle +\left(U({\bf{r}})-\mu +3g\langle \psi {\left({\bf{r}}\right)}^{2}\rangle -2g\langle \psi ({\bf{r}}){\rangle }^{2}\right)\langle \psi ({\bf{r}})\rangle \\ \qquad \qquad \quad +\,\langle u({\bf{r}})\psi ({\bf{r}})\rangle +\displaystyle \frac{40g}{3\sqrt{\pi }}{a}_{s}^{3/2}\left(3\langle \psi {\left({\bf{r}}\right)}^{2}{\rangle }^{2}-2\langle \psi ({\bf{r}}){\rangle }^{4}\right)=0,\end{array}\end{eqnarray} \tag{ 9 }$$

while the normalization condition (6) becomes

$$\begin{eqnarray}&&N={\int }_{{{\mathbb{R}}}^{3}}{{\rm{d}}}^{3}{\bf{r}}\left[\langle \psi {\left({\bf{r}}\right)}^{2}\rangle +\displaystyle \frac{8}{3\sqrt{\pi }}{a}_{s}^{3/2}\left(3\langle \psi ({\bf{r}})\rangle \langle \psi {\left({\bf{r}}\right)}^{2}\rangle -2\langle \psi ({\bf{r}}){\rangle }^{3}\right)\right].\end{eqnarray} \tag{ 10 }$$

Secondly, we assume that all the disorder-averaged quantities vary in space much more slowly than their fluctuations do. The smaller the correlation length with respect to the other characteristic lengths of the system, the better this assumption is justified. Therefore, in order to determine $\langle u({\bf{r}})\psi ({\bf{r}})\rangle$ and $\langle \psi {\left({\bf{r}}\right)}^{2}\rangle$ from (8) we employ the local density approximation (LDA) by treating $\langle \psi ({\bf{r}})\rangle$, $\langle \psi {\left({\bf{r}}\right)}^{2}\rangle$, etc, alongside with $U({\bf{r}})$ as constant background quantities $\langle \psi \rangle$, $\langle {\psi }^{2}\rangle$, and U, respectively, that will at the end of calculation be restored to their original values at the spatial point ${\bf{r}}$. Of course, $\langle \psi ({\bf{r}})\rangle$ and $\langle \psi {\left({\bf{r}}\right)}^{2}\rangle$ have to be determined in conjunction with (9) and (10). In this way, (8) under LDA describes the fluctuations of a spatially infinite condensate in a fixed self-consistent locally determined background. The final resulting equation, after applying the cumulant expansion as well, is given in appendix B.

Furthermore, according to the experimental data, the correlation length of the disorder potential in the propagation direction is about one order of magnitude larger than that in the transverse speckle plane. Hence, we will consider the scenario of a quasi two-dimensional disorder that is constant along the z-axis and has a Gaussian correlation function in the transverse x–y plane. After a lengthy but straightforward calculation, which is sketched in appendix B, we obtain

$$\begin{eqnarray}&&\langle u({\bf{r}})\psi ({\bf{r}})\rangle =\langle u\rangle \langle \psi ({\bf{r}})\rangle -\langle \psi ({\bf{r}})\rangle {\int }_{{{\mathbb{R}}}^{2}}\displaystyle \frac{{{\rm{d}}}^{2}{{\bf{k}}}_{\perp }}{{\left(2\pi \right)}^{2}}\displaystyle \frac{\langle u,u{\rangle }_{c}({{\bf{k}}}_{\perp })}{{\varepsilon }_{{{\bf{k}}}_{\perp }}+V({\bf{r}})+{V}_{{\rm{Q}}}({\bf{r}})},\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}&&\langle \psi {\left({\bf{r}}\right)}^{2}\rangle =\langle \psi ({\bf{r}}){\rangle }^{2}+\langle \psi ({\bf{r}}){\rangle }^{2}{\int }_{{{\mathbb{R}}}^{2}}\displaystyle \frac{{{\rm{d}}}^{2}{{\bf{k}}}_{\perp }}{{\left(2\pi \right)}^{2}}\displaystyle \frac{\langle u,u{\rangle }_{c}({{\bf{k}}}_{\perp })}{{\left({\varepsilon }_{{{\bf{k}}}_{\perp }}+V\left({\bf{r}}\right)+{V}_{{\rm{Q}}}\left({\bf{r}}\right)\right)}^{2}},\end{eqnarray} \tag{ 12 }$$

where ${\varepsilon }_{{{\bf{k}}}_{\perp }}\,=\,{{\hslash }}^{2}{{\bf{k}}}_{\perp }^{2}/(2m)={{\hslash }}^{2}({k}_{x}^{2}+{k}_{y}^{2})/(2m)$ denotes the transverse free-particle dispersion, $V({\bf{r}})=U({\bf{r}})\,+\langle u\rangle -\mu +3g\langle \psi {\left({\bf{r}}\right)}^{2}\rangle$ has the role of an effective potential and ${V}_{{\rm{Q}}}({\bf{r}})=\tfrac{160g}{3\sqrt{\pi }}{a}_{s}^{3/2}\langle \psi ({\bf{r}})\rangle \left(3\langle \psi {\left({\bf{r}}\right)}^{2}\rangle -2\langle \psi ({\bf{r}}){\rangle }^{2}\right)$ is a correction due to quantum fluctuations. The quantity $\langle u,u{\rangle }_{c}({{\bf{k}}}_{\perp })$ is the Fourier transform of the second-order cumulant of the disorder potential $\langle u,u{\rangle }_{c}({{\bf{r}}}_{\perp }^{}-{{\bf{r}}}_{\perp }^{{\prime} })\equiv \langle u({{\bf{r}}}_{\perp }^{}),u({{\bf{r}}}_{\perp }^{{\prime} }){\rangle }_{c}=\langle u({{\bf{r}}}_{\perp }^{})u({{\bf{r}}}_{\perp }^{{\prime} })\rangle -\langle u{\rangle }^{2}$. Had we considered higher-order cumulants, they would have as well contributed to (11) and (12). The set of equations (9)–(12) has to be solved upon specifying the total number of particles as well as the external and the disorder potential.

In accordance with the experimental situation, we assume that the Bose gas is trapped in an external harmonic potential

$$\begin{eqnarray}&&U({\bf{r}})=\displaystyle \frac{m}{2}\left({\omega }_{x}^{2}{x}^{2}+{\omega }_{y}^{2}{y}^{2}+{\omega }_{z}^{2}{z}^{2}\right),\end{eqnarray} \tag{ 13 }$$

where ω_x, ω_y, and ω_z are the corresponding trapping frequencies. In addition, we consider the isotropic 2D speckle disorder that has a Gaussian correlation function of the underlying electric field in the transverse plane

$$\begin{eqnarray}&&\gamma ({{\bf{r}}}_{\perp }^{})=\displaystyle \frac{\langle {{ \mathcal E }}^{* }({{\bf{r}}}_{\perp }^{}+{{\bf{r}}}_{\perp }^{{\prime} }){ \mathcal E }({{\bf{r}}}_{\perp }^{{\prime} })\rangle }{\langle | { \mathcal E }({{\bf{r}}}_{\perp }^{{\prime} }){| }^{2}\rangle }=\exp \left(-\displaystyle \frac{{{\bf{r}}}_{\perp }^{2}}{2{\sigma }^{2}}\right),\end{eqnarray} \tag{ 14 }$$

where σ is the correlation length. Since the resulting speckle potential is proportional to the square of the field, i.e. $u({\bf{r}})\propto | { \mathcal E }({\bf{r}}){| }^{2}$, it does not have a Gaussian distribution. However, due to the Gaussian character of the field ${ \mathcal E }({\bf{r}})$, all potential correlation functions can be expressed via sums of products of the field correlation function (14) [78]. Hence, we find

$$\begin{eqnarray}&&\langle u({{\bf{r}}}_{\perp }^{}),u({{\bf{r}}}_{\perp }^{{\prime} }){\rangle }_{c}=\langle u{\rangle }^{2}| \gamma ({{\bf{r}}}_{\perp }^{}-{{\bf{r}}}_{\perp }^{{\prime} }){| }^{2}=\langle u{\rangle }^{2}\exp \left[-\displaystyle \frac{{\left({{\bf{r}}}_{\perp }^{}-{{\bf{r}}}_{\perp }^{{\prime} }\right)}^{2}}{{\sigma }^{2}}\right],\end{eqnarray} \tag{ 15 }$$

where $\langle u\rangle =\langle u({{\bf{r}}}_{\perp }^{})\rangle$ is the average speckle potential value. In the transverse ${{\bf{k}}}_{\perp }$-plane, we get

$$\begin{eqnarray}&&\langle u,u{\rangle }_{c}({{\bf{k}}}_{\perp })=\pi {\sigma }^{2}\langle u{\rangle }^{2}\exp \left(-\displaystyle \frac{{\sigma }^{2}{{\bf{k}}}_{\perp }^{2}}{4}\right).\end{eqnarray} \tag{ 16 }$$

Now we aim for an approximate semi-analytic variational solution of the coupled set of equations (9)–(12). As no obvious energy functional exists, we cannot rely upon the standard variational procedure of [79, 80], which has been so successful in the realm of ultracold quantum gases. Instead, we follow [81, 82], where a similar variational solution is worked out directly on the basis of the underlying equations of motion within a so-called cumulant approach. To this end, we use a Gaussian ansatz for the disorder-ensemble averaged condensate wave function

$$\begin{eqnarray}&&\langle \psi ({\bf{r}})\rangle =\displaystyle \frac{A}{{\pi }^{3/4}\sqrt{{w}_{x}{w}_{y}{w}_{z}}}\exp \left(-\displaystyle \frac{{x}^{2}}{2{w}_{x}^{2}}-\displaystyle \frac{{y}^{2}}{2{w}_{y}^{2}}-\displaystyle \frac{{z}^{2}}{2{w}_{z}^{2}}\right),\end{eqnarray} \tag{ 17 }$$

as this reduces the complexity of the problem (9)–(12) to determining only five unknown parameters: the condensate widths w_x, w_y, and w_z, the normalization constant A, whose square quantifies the number of bosons in the condensate, and the chemical potential μ. Next, we multiply (9) by $\langle \psi ({\bf{r}})\rangle$, ${x}^{2}\langle \psi ({\bf{r}})\rangle$, ${y}^{2}\langle \psi ({\bf{r}})\rangle$, and ${z}^{2}\langle \psi ({\bf{r}})\rangle$, respectively, and integrate over the whole space. This provides four algebraic equations for the unknowns. The fifth equation is given by the normalization condition (10).

As a first result, we show in figure 5 how the condensate fraction A²/N monotonously decreases with increasing disorder strength $\overline{u}$, which equals $\langle u\rangle$, for the typical experimental parameters mentioned in section 2. For vanishing disorder we obtain a quantum depletion of 9% for the molecular BEC, which is much larger than for the usual atomic counterparts. Furthermore, in the range of experimentally accessible disorder strengths, i.e. up to $\overline{u}/{k}_{{\rm{B}}}=356\,\mathrm{nK}$, the condensate fraction only drops down to 62%. Even for much larger values of the disorder strength, our theory shows that the condensate fraction never decreases below 45%. With this we conclude that a significant part of the molecular Bose gas remains globally phase coherent despite the mutual impact of both quantum depletion and disorder fluctuations. Thus, as the condensate fraction is far from vanishing, a quantum phase transition from a superfluid to a Bose-glass does not occur for any disorder strength.

Zoom In Zoom Out Reset image size

In the following we focus our analysis to the cloud widths, which are directly accessible observables. To this end, we first have to calculate the cloud density profile, that is the integrand of (10). Then we find the 1/e half-widths along the x-, y- and z-axis of the density profile. In the subsequent section we compare the obtained theoretical cloud widths with the corresponding experimental ones.

A comparison of our measurements of the cloud widths along the x- and y-axis for different disorder strengths with theoretical predictions is presented in figure 6. According to section 3, our theory provides the disorder-ensemble averaged density profiles and their corresponding widths non-perturbatively with respect to disorder, including the quantum fluctuations. On the other hand, in the experiment we determine for each disorder realization the projected density profile and its widths. From the measured width distributions we obtain their average values and standard deviations, respectively. Thus, the experimental and theoretical procedures to obtain the cloud widths are not exactly the same. However, we have carefully checked that an additional analysis of the experimental data in the same manner as in the theoretical model would lead, for all practical purposes, to the same quantitative results for the average widths. Namely, according to the upper plot of figure 6(b), the cloud widths extracted from the averaged density profiles in the range of the considered disorder strengths turn out to deviate only up to 2% from the averages of the individual cloud widths. This deviation is smaller than the statistical error of 5% for extracting the average widths.

Zoom In Zoom Out Reset image size

In figure 6(a) we find that the geometric mean of the experimentally measured widths $\sqrt{{a}_{x}{a}_{y}}$, i.e. the cross-section area of the cloud in the imaging plane, is well reproduced by our theory. Additionally, our calculations reveal that the width a_z in the propagation direction also displays qualitatively the same growth. Hence, the cloud steadily expands while the disorder strength is increased. The corresponding particle number of 4.5 × 10⁵ is determined from the widths in the clean case, i.e. for $\overline{u}=0$, and is used in the subsequent calculations for all disorder strengths. Its agreement with the experimentally measured value (4.3 ± 0.2) × 10⁵ a posteriori justifies the necessity to include quantum fluctuations in our model. Had we neglected quantum fluctuation effects, the theoretically inferred particle number would have been 6.5 × 10⁵, which is about 45% larger than in the experiment.

In addition, for small disorder strengths the final equations for the unknown parameters w_x, w_y, w_z, A, and μ are solved perturbatively up to the second order in $\overline{u}$. Such an approach could be considered as an extension of the homogeneous dirty boson theory of Huang and Meng [30] to the trapped scenario. We find that the geometric mean of the measured transverse widths is reproduced quite well for $\overline{u}$/${k}_{{\rm{B}}}$ up to about 30 nK, see the green dashed curves in figure 6(a). To the best of our knowledge, this represents the first experimental verification of the Huang–Meng theory from 1992.

Furthermore, we observe in figure 6(b) that the averages of both transverse widths monotonously increase with the disorder strength in the experiment, as well as in the theoretical description. For small enough disorder we find again an agreement with the extended Huang–Meng-like theory. However, we have also to note that the present theory does not reproduce quantitatively the dependence of the individual widths for larger disorder strengths. The smaller width in x-direction is systematically underestimated theoretically, while the larger one in y-direction is overestimated. For the largest disorder strength the relative differences with respect to the measured values is about 12% for both widths. This noticeable discrepancy is attributed to the fact that the cloud aspect ratio turns out to decrease with increasing disorder strength, while our LDA-based theory predicts a constant aspect ratio, the same as in the clean-case, see the inset of figure 6(a). A physical explanation of this finding is that the disorder potential tends, on average, to make the cloud isotropic as the correlation lengths in both x- and y-direction coincide. Thus, the transverse isotropic disorder competes with the anisotropic trap and, by increasing its strength, the cloud becomes more and more circular. At the same time, the disorder counteracts the trap confinement and tends to spread out the cloud, which also explains the above mentioned disorder-induced cloud expansion. In our theoretical model, the cloud inevitably inherits the trap anisotropy due to the applied LDA. Hence, in order to capture the measured behavior, one must go beyond the present LDA by, for instance, including systematically the gradients of the confining potential $U({\bf{r}})$ and of the disorder-ensemble averages $\langle \psi ({\bf{r}})\rangle$ and $\langle \psi {\left({\bf{r}}\right)}^{2}\rangle$. However, that would make the calculations considerably more numerically demanding.

Beyond the so far presented comparison of the data common to both experiment and theory, further experimental measurements (theoretical considerations) can provide additional insights into the disordered Bose gas, that are currently not accessible to theory (experiment). Let us substantiate this conclusion by mentioning exemplarily the following two points.

First, we can easily extract theoretical values for disorder-ensemble averaged condensate widths, which elude any experimental observation as the condensate density is not measurable with the current setup. For weak disorder the condensate widths are closely following the cloud widths. In addition, for the largest experimentally accessible disorder strength $\overline{u}/{k}_{{\rm{B}}}=356\,\mathrm{nK}$ the condensate widths turn out to be still 91% of the cloud widths in both x- and y-direction. This finding is in line with the result from figure 5 that the condensate fraction remains to be above 60% in the performed experiment, even despite the presence of quite strong disorder.

Secondly, we examine the fluctuations of the cloud widths and their dependence on the disorder strength, since they are naturally accessible in the experiment. For all the strengths, the distributions of the widths can be well approximated by a Gaussian, see figure 7(a). The obtained distributions are compared based on the ratio of their standard deviation to the respective average value, normalized to the relative fluctuations in the clean case. With increasing disorder strength, the relative fluctuations increase approximately linearly for both transverse directions, as shown in figure 7(b). Notably, the smaller width in x-direction turns out to be more sensitive to the disorder variations than the larger one in y-direction. This is not unexpected, since the cloud is much narrower in x-direction and, hence, more disturbed by the laser speckle potential landscape. Note that the cloud size may vary not only due to different disorder realizations, but also due to the fluctuating particle number. Similarly as in the clean case Thomas–Fermi limit, we find that in the presence of disorder the widths are closely proportional to N^1/5. Therefore, we also present the relative fluctuations of N^1/5, given by the triangles in figure 7(b). The respective distribution is nearly constant and clearly excludes the possibility that the cloud widths fluctuate due to a change of the number of molecules.

Zoom In Zoom Out Reset image size

Within a theory-experiment collaboration we have performed a detailed case study of the dirty boson problem by analyzing the cloud shape of a molecular BEC in a disordered trap. The geometric mean of the experimentally measured transverse cloud widths, which corresponds to the cross-section area of the cloud in the imaging plane, turned out to perfectly agree within the error bars with a theory, which is non-perturbative with respect to disorder and includes quantum fluctuations. However, the experimentally measured and theoretically calculated individual cloud widths display deviations, which we attributed to the theoretical limitations of treating the harmonic confinement within LDA. Furthermore, we also investigated the fluctuations of the measured cloud widths and showed that they do not stem from particle number fluctuations, but only depend on the disorder strength. This experimental observation goes beyond the scope of the present theory, which only considers disorder-ensemble averaged quantities, without their corresponding fluctuations. Thus, we conclude that the performed experiment on a molecular BEC in a disordered trap paves the way and asks for a more refined theory in future studies.

We thank Benjamin Gänger and Jan Phieler for their help in initially constructing the experiment and Tobias Lausch for developing the software that runs the lab. Furthermore, we thank Antun Balaž, Piet Brouwer, and Aristeu Lima for insightful discussions. This work was supported by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) via the Collaborative Research Center SFB/TR185 (Project No. 277625399). BN receives support from a DFG Fellowship through the Excellence Initiative by the Graduate School Materials Science in Mainz (GSC 266).

Here, we introduce the basics of the cumulant expansion method and provide some examples that are used in our calculations. Let ${X}_{1}^{}$, ${X}_{2}^{}$, ..., ${X}_{n}^{}$ denote random variables. Their respective moments (averages) are defined by the moment generating function

$$\begin{eqnarray}\begin{array}{rcl}M({z}_{1}^{},{z}_{2}^{},\ldots ,{z}_{n}^{}) & = & \langle \exp ({z}_{1}^{}{X}_{1}^{}+{z}_{2}^{}{X}_{2}^{}\,+\ldots +\,{z}_{n}^{}{X}_{n}^{})\rangle \\ & = & \sum _{{\nu }_{1}=0}^{\infty }\sum _{{\nu }_{2}=0}^{\infty }\ldots \sum _{{\nu }_{n}=0}^{\infty }\displaystyle \frac{{z}_{1}^{{\nu }_{1}}}{{\nu }_{1}!}\displaystyle \frac{{z}_{2}^{{\nu }_{2}}}{{\nu }_{2}!}\ldots \displaystyle \frac{{z}_{n}^{{\nu }_{n}}}{{\nu }_{n}!}\langle {X}_{1}^{{\nu }_{1}}{X}_{2}^{{\nu }_{2}}\ldots {X}_{n}^{{\nu }_{n}}\rangle ,\end{array}\end{eqnarray} \tag{ A.1 }$$

while the cumulants are defined by the cumulant generating function

$$\begin{eqnarray}\begin{array}{rcl}K({z}_{1}^{},{z}_{2}^{},\ldots ,{z}_{n}^{}) & = & \mathrm{log}M({z}_{1}^{},{z}_{2}^{},\ldots ,{z}_{n}^{})\\ & = & \sum _{{\nu }_{1}=0}^{\infty }\sum _{{\nu }_{2}=0}^{\infty }\ldots \sum _{{\nu }_{n}=0}^{\infty\prime} \,\displaystyle \frac{{z}_{1}^{{\nu }_{1}}}{{\nu }_{1}!}\displaystyle \frac{{z}_{2}^{{\nu }_{2}}}{{\nu }_{2}!}\ldots \displaystyle \frac{{z}_{n}^{{\nu }_{n}}}{{\nu }_{n}!}\langle {X}_{1}^{{\nu }_{1}},{X}_{2}^{{\nu }_{2}},\ldots ,{X}_{n}^{{\nu }_{n}}{\rangle }_{c},\end{array}\end{eqnarray} \tag{ A.2 }$$

where the prime sign denotes that the term with ${\nu }_{1}={\nu }_{2}$ = ... = ${\nu }_{n}=0$ is excluded from the summation [75]. Moments can be expressed by cumulants of equal and lower orders and vice versa. For convenience, we note the following lowest-order examples:

$$\begin{eqnarray}&&\langle {X}_{1}^{}\rangle =\langle {X}_{1}^{}{\rangle }_{c},\end{eqnarray} \tag{ A.3 }$$

$$\begin{eqnarray}&&\langle {X}_{1}^{}{X}_{2}^{}\rangle =\langle {X}_{1}^{}\rangle \langle {X}_{2}^{}\rangle +\langle {X}_{1}^{},{X}_{2}^{}{\rangle }_{c},\end{eqnarray} \tag{ A.4 }$$

$$\begin{eqnarray}\begin{array}{rcl}\langle {X}_{1}^{}{X}_{2}^{}{X}_{3}^{}\rangle & = & \langle {X}_{1}^{}\rangle \langle {X}_{2}^{}\rangle \langle {X}_{3}^{}\rangle +\langle {X}_{3}^{}\rangle \langle {X}_{1}^{},{X}_{2}^{}{\rangle }_{c}+\langle {X}_{2}^{}\rangle \langle {X}_{1}^{},{X}_{3}^{}{\rangle }_{c}\\ & & +\,\langle {X}_{1}^{}\rangle \langle {X}_{2}^{},{X}_{3}^{}{\rangle }_{c}+\langle {X}_{1}^{},{X}_{2}^{},{X}_{3}^{}{\rangle }_{c},\end{array}\end{eqnarray} \tag{ A.5 }$$

$$\begin{eqnarray}\begin{array}{rcl}\langle {X}_{1}^{}{X}_{2}^{}{X}_{3}^{}{X}_{4}^{}\rangle & = & \langle {X}_{1}^{}\rangle \langle {X}_{2}^{}\rangle \langle {X}_{3}^{}\rangle \langle {X}_{4}^{}\rangle +\langle {X}_{3}^{}\rangle \langle {X}_{4}^{}\rangle \langle {X}_{1}^{},{X}_{2}^{}{\rangle }_{c}\\ & & +\,\langle {X}_{2}^{}\rangle \langle {X}_{4}^{}\rangle \langle {X}_{1}^{},{X}_{3}^{}{\rangle }_{c}+\langle {X}_{2}^{}\rangle \langle {X}_{3}^{}\rangle \langle {X}_{1}^{},{X}_{4}^{}{\rangle }_{c}\\ & & +\,\langle {X}_{1}^{}\rangle \langle {X}_{4}^{}\rangle \langle {X}_{2}^{},{X}_{3}^{}{\rangle }_{c}+\langle {X}_{1}^{},{X}_{4}^{}{\rangle }_{c}\langle {X}_{2}^{},{X}_{3}^{}{\rangle }_{c}\\ & & +\,\langle {X}_{1}^{}\rangle \langle {X}_{3}^{}\rangle \langle {X}_{2}^{},{X}_{4}^{}{\rangle }_{c}+\langle {X}_{1}^{},{X}_{3}^{}{\rangle }_{c}\langle {X}_{2}^{},{X}_{4}^{}{\rangle }_{c}\\ & & +\,\langle {X}_{1}^{}\rangle \langle {X}_{2}^{}\rangle \langle {X}_{3}^{},{X}_{4}^{}{\rangle }_{c}+\langle {X}_{1}^{},{X}_{2}^{}{\rangle }_{c}\langle {X}_{3}^{},{X}_{4}^{}{\rangle }_{c}\\ & & +\,\langle {X}_{4}^{}\rangle \langle {X}_{1}^{},{X}_{2}^{},{X}_{3}^{}{\rangle }_{c}+\langle {X}_{3}^{}\rangle \langle {X}_{1}^{},{X}_{2}^{},{X}_{4}^{}{\rangle }_{c}\\ & & +\,\langle {X}_{2}^{}\rangle \langle {X}_{1}^{},{X}_{3}^{},{X}_{4}^{}{\rangle }_{c}+\langle {X}_{1}^{}\rangle \langle {X}_{2}^{},{X}_{3}^{},{X}_{4}^{}{\rangle }_{c}\\ & & +\,\langle {X}_{1}^{},{X}_{2}^{},{X}_{3}^{},{X}_{4}^{}{\rangle }_{c}.\end{array}\end{eqnarray} \tag{ A.6 }$$

In the Gaussian approximation the cumulants of higher than second order are neglected, so that (A.5) and (A.6) reduce to

$$\begin{eqnarray}\begin{array}{rcl}\langle {X}_{1}^{}{X}_{2}^{}{X}_{3}^{}\rangle & \approx & \langle {X}_{1}^{}\rangle \langle {X}_{2}^{}\rangle \langle {X}_{3}^{}\rangle +\langle {X}_{3}^{}\rangle \langle {X}_{1}^{},{X}_{2}^{}{\rangle }_{c}+\langle {X}_{2}^{}\rangle \langle {X}_{1}^{},{X}_{3}^{}{\rangle }_{c}+\langle {X}_{1}^{}\rangle \langle {X}_{2}^{},{X}_{3}^{}{\rangle }_{c}\\ & = & \langle {X}_{3}^{}\rangle \langle {X}_{1}^{}{X}_{2}^{}\rangle +\langle {X}_{2}^{}\rangle \langle {X}_{1}^{}{X}_{3}^{}\rangle +\langle {X}_{1}^{}\rangle \langle {X}_{2}^{}{X}_{3}^{}\rangle -2\langle {X}_{1}^{}\rangle \langle {X}_{2}^{}\rangle \langle {X}_{3}^{}\rangle ,\end{array}\end{eqnarray} \tag{ A.7 }$$

$$\begin{eqnarray}\begin{array}{rcl}\langle {X}_{1}^{}{X}_{2}^{}{X}_{3}^{}{X}_{4}^{}\rangle & \approx & \langle {X}_{1}^{}\rangle \langle {X}_{2}^{}\rangle \langle {X}_{3}^{}\rangle \langle {X}_{4}^{}\rangle +\langle {X}_{3}^{}\rangle \langle {X}_{4}^{}\rangle \langle {X}_{1}^{},{X}_{2}^{}{\rangle }_{c}+\langle {X}_{2}^{}\rangle \langle {X}_{4}^{}\rangle \langle {X}_{1}^{},{X}_{3}^{}{\rangle }_{c}\\ & & +\,\langle {X}_{2}^{}\rangle \langle {X}_{3}^{}\rangle \langle {X}_{1}^{},{X}_{4}^{}{\rangle }_{c}+\langle {X}_{1}^{}\rangle \langle {X}_{4}^{}\rangle \langle {X}_{2}^{},{X}_{3}^{}{\rangle }_{c}+\langle {X}_{1}^{},{X}_{4}^{}{\rangle }_{c}\langle {X}_{2}^{},{X}_{3}^{}{\rangle }_{c}\\ & & +\,\langle {X}_{1}^{}\rangle \langle {X}_{3}^{}\rangle \langle {X}_{2}^{},{X}_{4}^{}{\rangle }_{c}+\langle {X}_{1}^{},{X}_{3}^{}{\rangle }_{c}\langle {X}_{2}^{},{X}_{4}^{}{\rangle }_{c}+\langle {X}_{1}^{}\rangle \langle {X}_{2}^{}\rangle \langle {X}_{3}^{},{X}_{4}^{}{\rangle }_{c}\\ & & +\,\langle {X}_{1}^{},{X}_{2}^{}{\rangle }_{c}\langle {X}_{3}^{},{X}_{4}^{}{\rangle }_{c}=\langle {X}_{1}^{}{X}_{3}^{}\rangle \langle {X}_{2}^{}{X}_{4}^{}\rangle +\langle {X}_{1}^{}{X}_{2}^{}\rangle \langle {X}_{3}^{}{X}_{4}^{}\rangle \\ & & +\,\langle {X}_{2}^{}{X}_{3}^{}\rangle \langle {X}_{1}^{}{X}_{4}^{}\rangle -2\langle {X}_{1}^{}\rangle \langle {X}_{2}^{}\rangle \langle {X}_{3}^{}\rangle \langle {X}_{4}^{}\rangle ,\end{array}\end{eqnarray} \tag{ A.8 }$$

where (A.3) and (A.4) were invoked to get the final results. As a special case, we obtain the approximations of the cubic and quartic moments of a random variable

$$\begin{eqnarray}&&\langle {X}_{1}^{3}\rangle \approx 3\langle {X}_{1}^{}\rangle \langle {X}_{1}^{2}\rangle -2\langle {X}_{1}^{}{\rangle }^{3},\end{eqnarray} \tag{ A.9 }$$

$$\begin{eqnarray}&&\langle {X}_{1}^{4}\rangle \approx 3\langle {X}_{1}^{2}{\rangle }^{2}-2\langle {X}_{1}^{}{\rangle }^{4}.\end{eqnarray} \tag{ A.10 }$$

For the sake of generality, we note that within the Gaussian approximation the following recursive relation holds

$$\begin{eqnarray}&&\langle {X}_{1}^{k}\rangle =\langle {X}_{1}^{}\rangle \langle {X}_{1}^{k-1}\rangle +(k-1)\left(\langle {X}_{1}^{2}\rangle -\langle {X}_{1}^{}{\rangle }^{2}\right)\langle {X}_{1}^{k-2}\rangle ,\quad k\geqslant 2,\end{eqnarray} \tag{ A.11 }$$

where $\langle {X}_{1}^{0}\rangle =1$.

Here, we briefly go through the main steps of our non-perturbative account of the disorder potential in a theory that also contains the quantum fluctuations corrections. Of course, the same procedure is applicable to a simple Gross–Pitaevskii equation [83]. Under the approximations mentioned in the main text, (8) becomes

$$\begin{eqnarray}&&\begin{array}{l}-\displaystyle \frac{{{\hslash }}^{2}}{2m}{{\rm{\nabla }}}_{\perp }^{2}\varphi ({{\bf{r}}}_{\perp }^{})+U\varphi ({{\bf{r}}}_{\perp }^{})+(u({{\bf{r}}}_{\perp }^{})-\langle u\rangle )\langle \psi \rangle +(u({{\bf{r}}}_{\perp }^{})\varphi ({{\bf{r}}}_{\perp }^{})-\langle u({{\bf{r}}}_{\perp }^{})\varphi ({{\bf{r}}}_{\perp }^{})\rangle )-\mu \varphi ({{\bf{r}}}_{\perp }^{})\\ \quad +\,g\,\left\{3\langle \psi {\rangle }^{2}\varphi ({{\bf{r}}}_{\perp }^{})+3\langle \psi \rangle \left(\varphi {\left({{\bf{r}}}_{\perp }^{}\right)}^{2}-\langle {\varphi }^{2}\rangle \right)+\varphi {\left({{\bf{r}}}_{\perp }^{}\right)}^{3}+\displaystyle \frac{40}{3\sqrt{\pi }}{a}_{s}^{3/2}\left[4\langle \psi {\rangle }^{3}\varphi ({{\bf{r}}}_{\perp }^{})\right.\right.\\ \quad \left.\left.+\,6\langle \psi {\rangle }^{2}\left(\varphi {\left({{\bf{r}}}_{\perp }^{}\right)}^{2}-\langle {\varphi }^{2}\rangle \right)+4\langle \psi \rangle \varphi {\left({{\bf{r}}}_{\perp }^{}\right)}^{3}+\left(\varphi {\left({{\bf{r}}}_{\perp }^{}\right)}^{4}-3\langle {\varphi }^{2}{\rangle }^{2}\right)\right]\right\}=0,\end{array}\end{eqnarray} \tag{ B.1 }$$

where now there is only a transverse spatial dependence due to the assumed two-dimensional character of the disorder and due to the LDA. Hence, we are here effectively solving the problem of a homogeneous two-dimensional disordered Bose gas. An extension to three dimensions is straightforward. Note that the cumulant expansion has also been used. Multiplying the previous equation with $u({{\bf{r}}}_{\perp }^{{\prime} })$ and averaging produces

$$\begin{eqnarray}&&\begin{array}{l}-\displaystyle \frac{{{\hslash }}^{2}}{2m}{{\rm{\nabla }}}_{\perp }^{2}\langle u({{\bf{r}}}_{\perp }^{{\prime} })\varphi ({{\bf{r}}}_{\perp }^{})\rangle +U\langle u({{\bf{r}}}_{\perp }^{{\prime} })\varphi ({{\bf{r}}}_{\perp }^{})\rangle +(\langle u({{\bf{r}}}_{\perp }^{{\prime} })u({{\bf{r}}}_{\perp }^{})\rangle -\langle u{\rangle }^{2})\langle \psi \rangle \\ \quad +\,(\langle u({{\bf{r}}}_{\perp }^{{\prime} })u({{\bf{r}}}_{\perp }^{})\varphi ({{\bf{r}}}_{\perp }^{})\rangle -\langle u\rangle \langle u({{\bf{r}}}_{\perp }^{})\varphi ({{\bf{r}}}_{\perp }^{})\rangle )-\mu \langle u({{\bf{r}}}_{\perp }^{{\prime} })\varphi ({{\bf{r}}}_{\perp }^{})\rangle \\ \quad +\,g\,\left\{3\langle \psi {\rangle }^{2}\langle u({{\bf{r}}}_{\perp }^{{\prime} })\varphi ({{\bf{r}}}_{\perp }^{})\rangle +3\langle \psi \rangle (\langle u({{\bf{r}}}_{\perp }^{{\prime} })\varphi {\left({{\bf{r}}}_{\perp }^{}\right)}^{2}\rangle -\langle u\rangle \langle {\varphi }^{2}\rangle )+\langle u({{\bf{r}}}_{\perp }^{{\prime} })\varphi {\left({{\bf{r}}}_{\perp }^{}\right)}^{3}\rangle \right.\\ \quad +\,\displaystyle \frac{40}{3\sqrt{\pi }}{a}_{s}^{3/2}\left[4\langle \psi {\rangle }^{3}\langle u({{\bf{r}}}_{\perp }^{{\prime} })\varphi ({{\bf{r}}}_{\perp }^{})\rangle +6\langle \psi {\rangle }^{2}\left(\langle u({{\bf{r}}}_{\perp }^{{\prime} })\varphi {\left({{\bf{r}}}_{\perp }^{}\right)}^{2}\rangle -\langle u\rangle \langle {\varphi }^{2}\rangle \right)\right.\\ \quad +\,\left.\left.4\langle \psi \rangle \langle u({{\bf{r}}}_{\perp }^{{\prime} })\varphi {\left({{\bf{r}}}_{\perp }^{}\right)}^{3}\rangle +\left(\langle u({{\bf{r}}}_{\perp }^{{\prime} })\varphi {\left({{\bf{r}}}_{\perp }^{}\right)}^{4}\rangle -3\langle u\rangle \langle {\varphi }^{2}{\rangle }^{2}\right)\right]\right\}=0,\end{array}\end{eqnarray} \tag{ B.2 }$$

which reduces after the cumulant expansion to

$$\begin{eqnarray}&&-\displaystyle \frac{{{\hslash }}^{2}}{2m}{{\rm{\nabla }}}_{\perp }^{2}\langle u({{\bf{r}}}_{\perp }^{{\prime} }),\varphi ({{\bf{r}}}_{\perp }^{}){\rangle }_{c}+(V+{V}_{{\rm{Q}}})\langle u({{\bf{r}}}_{\perp }^{{\prime} }),\varphi ({{\bf{r}}}_{\perp }^{}){\rangle }_{c}+\langle u({{\bf{r}}}_{\perp }^{{\prime} }),u({{\bf{r}}}_{\perp }^{}){\rangle }_{c}\langle \psi \rangle =0,\end{eqnarray} \tag{ B.3 }$$

where $V=U+\langle u\rangle -\mu +3g\langle {\psi }^{2}\rangle$ and ${V}_{{\rm{Q}}}=\tfrac{160g}{3\sqrt{\pi }}{a}_{s}^{3/2}\langle \psi \rangle \left(3\langle {\psi }^{2}\rangle -2\langle \psi {\rangle }^{2}\right)$ have been introduced in the main text. Similarly, multiplying (B.1) with $\varphi ({{\bf{r}}}_{\perp }^{{\prime} })$ and averaging gives

$$\begin{eqnarray}&&\begin{array}{l}-\displaystyle \frac{{{\hslash }}^{2}}{2m}{{\rm{\nabla }}}_{\perp }^{2}\langle \varphi ({{\bf{r}}}_{\perp }^{{\prime} })\varphi ({{\bf{r}}}_{\perp }^{})\rangle +U\langle \varphi ({{\bf{r}}}_{\perp }^{{\prime} })\varphi ({{\bf{r}}}_{\perp }^{})\rangle +\langle \varphi ({{\bf{r}}}_{\perp }^{{\prime} })u({{\bf{r}}}_{\perp }^{})\rangle \langle \psi \rangle +\langle \varphi ({{\bf{r}}}_{\perp }^{{\prime} })u({{\bf{r}}}_{\perp }^{})\varphi ({{\bf{r}}}_{\perp }^{})\rangle \\ \quad -\,\mu \langle \varphi ({{\bf{r}}}_{\perp }^{{\prime} })\varphi ({{\bf{r}}}_{\perp }^{})\rangle +g\,\left\{3\langle \psi {\rangle }^{2}\langle \varphi ({{\bf{r}}}_{\perp }^{{\prime} })\varphi ({{\bf{r}}}_{\perp }^{})\rangle +3\langle \psi \rangle \langle \varphi ({{\bf{r}}}_{\perp }^{{\prime} })\varphi {\left({{\bf{r}}}_{\perp }^{}\right)}^{2}\rangle \right.\\ \quad +\langle \varphi ({{\bf{r}}}_{\perp }^{{\prime} })\varphi {\left({{\bf{r}}}_{\perp }^{}\right)}^{3}\rangle +\displaystyle \frac{40}{3\sqrt{\pi }}{a}_{s}^{3/2}\left[4\langle \psi {\rangle }^{3}\langle \varphi ({{\bf{r}}}_{\perp }^{{\prime} })\varphi ({{\bf{r}}}_{\perp }^{})\rangle +6\langle \psi {\rangle }^{2}\langle \varphi ({{\bf{r}}}_{\perp }^{{\prime} })\varphi {\left({{\bf{r}}}_{\perp }^{}\right)}^{2}\rangle \right.\\ \quad +\,\left.\left.4\langle \psi \rangle \langle \varphi ({{\bf{r}}}_{\perp }^{{\prime} })\varphi {\left({{\bf{r}}}_{\perp }^{}\right)}^{3}\rangle +\langle \varphi ({{\bf{r}}}_{\perp }^{{\prime} })\varphi {\left({{\bf{r}}}_{\perp }^{}\right)}^{4}\rangle \right]\right\}=0,\end{array}\end{eqnarray} \tag{ B.4 }$$

which simplifies after the cumulant expansion to

$$\begin{eqnarray}&&-\displaystyle \frac{{{\hslash }}^{2}}{2m}{{\rm{\nabla }}}_{\perp }^{2}\langle \varphi ({{\bf{r}}}_{\perp }^{{\prime} }),\varphi ({{\bf{r}}}_{\perp }^{}){\rangle }_{c}+(V+{V}_{{\rm{Q}}})\langle \varphi ({{\bf{r}}}_{\perp }^{{\prime} }),\varphi ({{\bf{r}}}_{\perp }^{}){\rangle }_{c}+\langle u({{\bf{r}}}_{\perp }^{{\prime} }),\varphi ({{\bf{r}}}_{\perp }^{}){\rangle }_{c}\langle \psi \rangle =0.\end{eqnarray} \tag{ B.5 }$$

Note that in the last term we replaced $\langle u({{\bf{r}}}_{\perp }^{}),\varphi ({{\bf{r}}}_{\perp }^{{\prime} }){\rangle }_{c}$ with $\langle u({{\bf{r}}}_{\perp }^{{\prime} }),\varphi ({{\bf{r}}}_{\perp }^{}){\rangle }_{c}$, which should be exact in the homogeneous case. The obtained two coupled equations (B.3) and (B.5) can then be rewritten as

$$\begin{eqnarray}&&-\displaystyle \frac{{{\hslash }}^{2}}{2m}{{\rm{\nabla }}}_{\perp }^{2}\langle u,\varphi {\rangle }_{c}({{\bf{r}}}_{\perp }^{{\prime} }-{{\bf{r}}}_{\perp }^{})+(V+{V}_{{\rm{Q}}})\langle u,\varphi {\rangle }_{c}({{\bf{r}}}_{\perp }^{{\prime} }-{{\bf{r}}}_{\perp }^{})+\langle u,u{\rangle }_{c}({{\bf{r}}}_{\perp }^{{\prime} }-{{\bf{r}}}_{\perp }^{})\langle \psi \rangle =0,\end{eqnarray} \tag{ B.6 }$$

$$\begin{eqnarray}&&-\displaystyle \frac{{{\hslash }}^{2}}{2m}{{\rm{\nabla }}}_{\perp }^{2}\langle \varphi ,\varphi {\rangle }_{c}({{\bf{r}}}_{\perp }^{{\prime} }-{{\bf{r}}}_{\perp }^{})+(V+{V}_{{\rm{Q}}})\langle \varphi ,\varphi {\rangle }_{c}({{\bf{r}}}_{\perp }^{{\prime} }-{{\bf{r}}}_{\perp }^{})+\langle u,\varphi {\rangle }_{c}({{\bf{r}}}_{\perp }^{{\prime} }-{{\bf{r}}}_{\perp }^{})\langle \psi \rangle =0,\end{eqnarray} \tag{ B.7 }$$

due to spatial homogeneity after the disorder averaging. In the transverse ${{\bf{k}}}_{\perp }$-plane one finds

$$\begin{eqnarray}&&({\varepsilon }_{{{\bf{k}}}_{\perp }}+V+{V}_{{\rm{Q}}})\langle u,\varphi {\rangle }_{c}({{\bf{k}}}_{\perp })+\langle u,u{\rangle }_{c}({{\bf{k}}}_{\perp })\langle \psi \rangle =0,\end{eqnarray} \tag{ B.8 }$$

$$\begin{eqnarray}&&({\varepsilon }_{{{\bf{k}}}_{\perp }}+V+{V}_{{\rm{Q}}})\langle \varphi ,\varphi {\rangle }_{c}({{\bf{k}}}_{\perp })+\langle u,\varphi {\rangle }_{c}({{\bf{k}}}_{\perp })\langle \psi \rangle =0,\end{eqnarray} \tag{ B.9 }$$

where ${\varepsilon }_{{{\bf{k}}}_{\perp }}={{\hslash }}^{2}{{\bf{k}}}_{\perp }^{2}/(2m)={{\hslash }}^{2}({k}_{x}^{2}+{k}_{y}^{2})/(2m)$ denotes the transverse dispersion, so that

$$\begin{eqnarray}&&\langle u,\varphi {\rangle }_{c}({{\bf{k}}}_{\perp })=-\displaystyle \frac{\langle \psi \rangle \langle u,u{\rangle }_{c}({{\bf{k}}}_{\perp })}{{\varepsilon }_{{{\bf{k}}}_{\perp }}+V+{V}_{{\rm{Q}}}},\end{eqnarray} \tag{ B.10 }$$

$$\begin{eqnarray}&&\langle \varphi ,\varphi {\rangle }_{c}({{\bf{k}}}_{\perp })=\displaystyle \frac{\langle \psi {\rangle }^{2}\langle u,u{\rangle }_{c}({{\bf{k}}}_{\perp })}{{\left({\varepsilon }_{{{\bf{k}}}_{\perp }}+V+{V}_{{\rm{Q}}}\right)}^{2}}.\end{eqnarray} \tag{ B.11 }$$

Integration over the entire ${{\bf{k}}}_{\perp }$-plane of the previous equations finally yields

$$\begin{eqnarray}&&\langle u({{\bf{r}}}_{\perp }^{}),\varphi ({{\bf{r}}}_{\perp }^{}){\rangle }_{c}=-\langle \psi \rangle {\int }_{{{\mathbb{R}}}^{2}}\displaystyle \frac{{{\rm{d}}}^{2}{{\bf{k}}}_{\perp }}{{\left(2\pi \right)}^{2}}\displaystyle \frac{\langle u,u{\rangle }_{c}({{\bf{k}}}_{\perp })}{{\varepsilon }_{{{\bf{k}}}_{\perp }}+V+{V}_{{\rm{Q}}}},\end{eqnarray} \tag{ B.12 }$$

$$\begin{eqnarray}&&\langle \varphi ({{\bf{r}}}_{\perp }^{}),\varphi ({{\bf{r}}}_{\perp }^{}){\rangle }_{c}=\langle \psi {\rangle }^{2}{\int }_{{{\mathbb{R}}}^{2}}\displaystyle \frac{{{\rm{d}}}^{2}{{\bf{k}}}_{\perp }}{{\left(2\pi \right)}^{2}}\displaystyle \frac{\langle u,u{\rangle }_{c}({{\bf{k}}}_{\perp })}{{\left({\varepsilon }_{{{\bf{k}}}_{\perp }}+V+{V}_{{\rm{Q}}}\right)}^{2}}.\end{eqnarray} \tag{ B.13 }$$

The obtained cumulants are homogeneous across the entire ${{\bf{r}}}_{\perp }^{}$-space and depend parametrically on $\langle \psi \rangle ,\langle {\psi }^{2}\rangle$ and U. In a trapped system according to the LDA, we get

$$\begin{eqnarray}&&\langle u({\bf{r}}),\varphi ({\bf{r}}){\rangle }_{c}=-\langle \psi ({\bf{r}})\rangle {\int }_{{{\mathbb{R}}}^{2}}\displaystyle \frac{{{\rm{d}}}^{2}{{\bf{k}}}_{\perp }}{{\left(2\pi \right)}^{2}}\displaystyle \frac{\langle u,u{\rangle }_{c}({{\bf{k}}}_{\perp })}{{\varepsilon }_{{{\bf{k}}}_{\perp }}+V({\bf{r}})+{V}_{{\rm{Q}}}({\bf{r}})},\end{eqnarray} \tag{ B.14 }$$

$$\begin{eqnarray}&&\langle \varphi ({\bf{r}}),\varphi ({\bf{r}}){\rangle }_{c}=\langle \psi ({\bf{r}}){\rangle }^{2}{\int }_{{{\mathbb{R}}}^{2}}\displaystyle \frac{{{\rm{d}}}^{2}{{\bf{k}}}_{\perp }}{{\left(2\pi \right)}^{2}}\displaystyle \frac{\langle u,u{\rangle }_{c}({{\bf{k}}}_{\perp })}{{\left({\varepsilon }_{{{\bf{k}}}_{\perp }}+V\left({\bf{r}}\right)+{V}_{{\rm{Q}}}\left({\bf{r}}\right)\right)}^{2}},\end{eqnarray} \tag{ B.15 }$$

by restoring the explicit spatial dependence of $\langle \psi ({\bf{r}})\rangle ,\langle {\psi }^{2}({\bf{r}})\rangle$ and $U({\bf{r}})$. The last two equations correspond to (11) and (12) of the main text.