Analytical and numerical study of dirty bosons in a quasi-one-dimensional harmonic trap

The dirty boson problem is defined as a system of interacting bosons in a random potential [1]. The combined effect of disorder and two-particle interaction represents one of the most challenging problems in condensed matter physics due to the intriguing interplay between localisation and superfluidity. Cold atoms provide a controlled experimental setup in which that fundamental question of interacting bosons in a random environment can be addressed in both a quantitative and a tunable way.

The earliest relevant experiments, which were central for motivating the research of the dirty boson problem, dealt with superfluidity of thin films of ⁴He adsorbed in porous Vycor glass in the low-density limit [2]. There it was proven that, despite the presence of disorder, superfluidity can still persist. For ultracold Bose gases disorder appears either naturally as, e.g., in magnetic wire traps [3–7], where imperfections of the wire itself can induce local disorder, or it may be created artificially and controllably as, e.g., by using laser speckle fields [8–12]. A setup more in the spirit of condensed matter physics relies on a Bose gas with impurity atoms of another species trapped in a deep optical lattice, so the latter represent randomly distributed scatterers [13, 14]. Furthermore, an incommensurate optical lattice can provide a pseudo-random potential for an ultracold Bose gas [15–17].

Non-interacting particles in a random environment can be localised provided that the disorder is sufficiently strong. This phenomenon of Anderson localisation occurs as the particles are repeatedly reflected back in the random potential, so interferences yield exponentially localised one-body wave functions [18]. In one dimension Anderson localisation was experimentally found in an ultracold Bose gas in [11, 17]. Within a Bose–Einstein condensation (BEC), which is a many-particle interacting system, the presence of disorder causes the emergence of a new phase besides the superfluid phase (SF), which is called a Bose-glass phase due to the localisation of bosons in the respective minima of the random potential landscape. This Bose-glass phase contains no superfluid fraction and is characterised by a finite compressibility, by the absence of a gap, and by an infinite superfluid susceptibility [1]. Indications for the existence of the Bose-glass phase were found, for instance, in the experiments of [6, 7, 9, 19]. There it was shown within the superfluid phase that an increasing disorder strength first yields a fragmentation of the condensate due to the formation of tiny BEC droplets in the minima of the random environment. For a sufficiently strong disorder the condensate then turns out to be completely destroyed as all bosons are localised in the minima of the random potential, which represents the Bose-glass phase. But a more quantitative investigation of that elusive phase is still lacking both from an experimental and a theoretical point of view.

One of the first important theoretical results of the dirty boson problem was obtained by Huang and Meng in 1992 [20]. Within a Bogoliubov theory for a weakly interacting Bose–Einstein condensate it was found that a weak random disorder potential leads to a depletion of both the global condensate density and the superfluid density due to the localisation of bosons in the respective minima of the random potential [20–31]. Beyond the weak disorder, a perturbative approach was worked out in [32, 33], where the impact of the external random potential upon the quantum fluctuations was studied in detail. In order to analyse the BEC in the strong disorder regime, there are, in principle, two complementary non-perturbative approaches. The first starts from the superfluid phase and ends up in the Bose-glass phase for increasing disorder strength. To this end [34] applies the random phase approximation and yields a self-consistent integral equation for the disorder-averaged particle density, whereas [35, 36] work out a stochastic self-consistent mean-field approach using two chemical potentials, one for the condensed and one for the exited particles. The second approach starts conversely from the Bose-glass phase and proceeds towards the superfluid phase for decreasing disorder strength. For instance, [37, 38] work this out on the basis of a careful energetic analysis, where the disorder turns out to strongly influence the size, shape, and structure of the mini-condensates in the minima of the random potential. Furthermore, an order parameter was introduced and applied in the context of a Hartree–Fock theory in [39] in order to describe the possible emergence of the Bose-glass phase.

The localisation due the presence of the disorder in one-dimensional BEC systems was tackled in the literature in different directions. For instance, it was shown analytically that the one-dimensional gas of short-range-interacting bosons in the presence of disorder can undergo a finite-temperature phase transition between superfluid and insulator [40]. Furthermore, solving numerically and variationally the Gross–Pitaevskii equation of the weakly interacting BEC in a weakly disordered lattice and a speckle potential, the localised BECs are found to have an exponential tail [41]. Using quantum Monte Carlo simulations it turned out that, surprisingly, disorder-induced phase coherence could occur [42]. Furthermore, in the context of optical lattices, the quantum phase diagram of a dirty BEC in one dimension was also investigated via different methods [1]. [43] proved that approximative description of all quantum phases can be obtained via the site-dependent decoupling mean-field approach. By means of the density matrix renormalisation group technique, the existence of a critical value of the disorder strength for the Bose-glass phase was proven in [44]. In addition, the exact Bose-Fermi mapping demonstrated that the superfluid Bose-glass transition and the general phase diagram of trapped incommensurate optical lattices can be uniquely determined from finite-temperature density distributions of the trapped disordered system [45]. Despite all those previous investigations there is still a lack of knowledge concerning the emergence of the Bose-glass phase and its elusive properties.

In the present paper we treat a quasi-one-dimensional trapped BEC in a disorder potential both analytically and numerically. In particular, we focus on the question how the bosons, which are localised in the minima of the random potential, are distributed within the harmonic confinement. For sufficiently large disorder we even expect to find a Bose-glass region in the trap, where the global condensate vanishes and only localised bosons exist. Note that the corresponding three-dimensional trapped case is treated separately in [46]. We start first by describing the underlying BEC model and by developing a Hartree–Fock mean-field theory for the weak disorder regime and apply it within the Thomas–Fermi approximation to the dirty BEC system in section 2. However, since we study a system that is not fully amenable to the Thomas–Fermi approximation, we also employ numerical and variational treatment, described in section 3. We solve the corresponding Gross–Pitaevskii (GP) equation of the BEC model, and then apply a variational ansatz for the intermediate disorder regime. The results of those three different methods are discussed and compared in section 4. For instance, we find that the density of fragmented mini-condensates is redistributed for increasing disorder strength. Whereas for small disorder bosons tend to localise at the border of the trap, for intermediate disorder strength they concentrate in the trap centre.

It has been suggested in [22] that a Gaussian-correlated disorder potential constitutes an appropriate model to describe realistic random landscapes. Therein the final correlation length corresponds to the average width of the mountains or valleys in the disorder potential. Such a Gaussian correlation has been explored in more detail both for a BEC with contact as well as dipole-dipole interaction in any geometry [25, 29, 30]. Qualitatively similar results are obtained for other disorder potentials with finite correlation lengths as, for instance, laser speckles [26, 29] and Lorentz correlation [27]. All these studies have in common that disorder effects typically decrease with increasing correlation length and are, thus, most pronounced for δ-correlation. Therefore, we restrict ourselves in the following to the case of δ-correlation, to focus on the study of disorder effects.

The model of a three-dimensional weakly interacting homogeneous Bose gas in a δ-correlated disorder potential was studied within the Hartree–Fock mean-field theory in [39] by applying the Parisi replica method [47–49]. As a result, the corresponding phase diagram for the occurrence of the superfluid, the Bose-glass, and the normal phase was determined in the control parameter plane spanned by disorder strength and temperature. This Hartree–Fock theory is extended in [50] to a harmonic confinement, and is applied in the following to one-dimensional systems.

To this end, we consider a model of one-dimensional harmonically trapped BEC in a δ-correlated disorder potential with contact interactions between the particles. The corresponding free energy is calculated in appendix A and is given by equation (A5). It depends on the superfluid order parameter ${n}_{0}(x)$, which represents the condensate density, the Bose-glass order parameter q(x), which stands for the density of atoms being localised in the local minima of disorder potential, and an auxiliary function ${Q}_{0}(x)$. The self-consistency equations are obtained by extremizing the free energy with respect to these functions, i.e., $\tfrac{\delta { \mathcal F }}{\delta {n}_{0}({x}^{\prime })}=0$, $\tfrac{\delta { \mathcal F }}{\delta {Q}_{0}({x}^{\prime })}=0$, and $\tfrac{\delta { \mathcal F }}{\delta q({x}^{\prime })}=0$. This yields, together with the particle number equation (A6), four coupled equations: an algebraic equation for q(x),

$$\begin{eqnarray}&&q(x)=\displaystyle \frac{D}{{\hslash }M}{Q}_{0}{(x)}^{3}\;\displaystyle \frac{{n}_{0}(x)}{1-\tfrac{D}{{\hslash }M}{Q}_{0}{(x)}^{3}},\end{eqnarray} \tag{ 1 }$$

a nonlinear differential equation for ${n}_{0}(x)$,

$$\begin{eqnarray}&&\left[-\mu +2{gn}(x)+V(x)-{{gn}}_{0}(x)-\displaystyle \frac{D}{{\hslash }}{Q}_{0}(x)-\displaystyle \frac{{{\hslash }}^{2}}{2M}\displaystyle \frac{{\partial }^{2}}{\partial {x}^{2}}\right]\;\sqrt{{n}_{0}(x)}=0,\end{eqnarray} \tag{ 2 }$$

the equation for the total density n(x), which is sum of the previous two densities,

$$\begin{eqnarray}&&n(x)=q(x)+{n}_{0}(x),\end{eqnarray} \tag{ 3 }$$

and the auxiliary function ${Q}_{0}(x)$, which is a solution of the cubic equation

$$\begin{eqnarray}&&-\displaystyle \frac{D}{{\hslash }}{Q}_{0}{(x)}^{3}+[-\mu +2{gn}(x)+V(x)]{Q}_{0}{(x)}^{2}-\displaystyle \frac{M}{2}=0.\end{eqnarray} \tag{ 4 }$$

Here D denotes the strength of disorder, as defined in appendix A. Note that in the clean case (D = 0) equations (1)–(4) reduce to the standard Gross–Pitaevskii theory. Furthermore, for finite disorder strength D the homogeneous case, where $V(x)=0$, is treated semi-analytically in appendix B. There it is shown that increasing the disorder strength D yields a first-order quantum phase transition from the superfluid to the Bose-glass phase.

For the harmonically trapped case, however, no analytic approach is known which gives the exact solution of the differential equation (2) even in the absence of disorder. Therefore, we approximate its solution via the Thomas–Fermi (TF) approximation method, which is based on neglecting the kinetic energy. To this end, it turns out that we have to distinguish between two different spatial regions: the superfluid region, where the bosons are distributed in the condensate as well as in the minima of the disorder potential, and the Bose-glass region, where there are no bosons in the global condensate so that all bosons contribute only to the local Bose–Einstein condensates. In the following the radius of the superfluid region, i.e., the condensate radius, is denoted by ${R}_{{\rm{TF1}}}$, while the radius of the whole bosonic cloud ${R}_{{\rm{TF2}}}$ is called the cloud radius.

Within the TF approximation the algebraic equations (1), (3), and (4) remain the same, but the differential equation (2) reduces to an algebraic relation in the superfluid region:

$$\begin{eqnarray}&&-\mu +2{gn}(x)+V(x)-{{gn}}_{0}(x)-\displaystyle \frac{D}{{\hslash }}{Q}_{0}(x)=0.\end{eqnarray} \tag{ 5 }$$

Outside the superfluid region, i.e., in the Bose-glass region, equation (2) reduces simply to ${n}_{0}(x)=0$. The advantage of the TF approximation is that now we have only four coupled algebraic equations.

At first we consider the superfluid region. In the TF approximation the dependency on the auxiliary function ${Q}_{0}(x)$ in equation (4) can be eliminated and equations (1), (3), and (5) reduce in the superfluid region to:

$$\begin{eqnarray}&&-\tilde{\mu }+2\tilde{n}(\tilde{x})+{\tilde{x}}^{2}-{\tilde{n}}_{0}(\tilde{x})-2\displaystyle \frac{\tilde{D}}{\sqrt{{\tilde{n}}_{0}(\tilde{x})}}=0,\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&\tilde{q}(\tilde{x})=\tilde{D}\;\displaystyle \frac{{\tilde{n}}_{0}(\tilde{x})}{{\tilde{n}}_{0}{(\tilde{x})}^{3/2}-\tilde{D}},\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&\tilde{n}(\tilde{x})=\tilde{q}(\tilde{x})+{\tilde{n}}_{0}(\tilde{x}),\end{eqnarray} \tag{ 8 }$$

where $\tilde{n}(\tilde{x})=n(x)/\bar{n}$ denotes the dimensionless total density, ${\tilde{n}}_{0}(\tilde{x})={n}_{0}(x)/\bar{n}$ the dimensionless condensate density, $\tilde{q}(\tilde{x})=q(x)/\bar{n}$ the dimensionless Bose-glass order parameter, $\tilde{x}=x/{R}_{\mathrm{TF}}$ the dimensionless coordinate, $\bar{n}=\bar{\mu }/g$ the maximal total density in the clean case, $\tilde{\mu }=\mu /\bar{\mu }$ the dimensionless chemical potential, $\tilde{D}=\tfrac{{\xi }^{3}}{{{ \mathcal L }}^{3}}$ the dimensionless disorder strength, $\xi =\tfrac{{l}^{2}}{{R}_{\mathrm{TF}}}$ the coherence length in the centre of the trap, $l=\sqrt{\tfrac{{\hslash }}{M{\rm{\Omega }}}}$ the oscillator length, and ${R}_{\mathrm{TF}}=l\sqrt{\tfrac{2\bar{\mu }}{{\hslash }{\rm{\Omega }}}}$ the TF cloud radius in the clean case, i.e., when D = 0. The chemical potential in the absence of the disorder $\bar{\mu }={\hslash }{\omega }_{r}{\left(\tfrac{3}{2\sqrt{2}}N\tfrac{a}{l}\sqrt{\tfrac{{\rm{\Omega }}}{{\omega }_{r}}}\right)}^{2/3}$, which provides the energy scale, is deduced from the normalisation condition (A6) in the clean case.

Inserting equations (7) and (8) into equation (6) gives us one self-consistency equation for the condensate density in the superfluid region:

$$\begin{eqnarray}&&{\tilde{n}}_{0}{(\tilde{x})}^{3}+(-\tilde{\mu }+{\tilde{x}}^{2}){\tilde{n}}_{0}{(\tilde{x})}^{2}-\tilde{D}{\tilde{n}}_{0}{(\tilde{x})}^{3/2}-\tilde{D}(-\tilde{\mu }+{\tilde{x}}^{2})\sqrt{{\tilde{n}}_{0}(\tilde{x})}+2{\tilde{D}}^{2}=0.\end{eqnarray} \tag{ 9 }$$

This equation is of sixth order with respect to $\sqrt{{\tilde{n}}_{0}(\tilde{x})}$, which makes it impossible to solve analytically. Therefore, we solve it numerically and insert the result into equations (7) and (8) in order to determine the Bose-glass order parameter $\tilde{q}(\tilde{x})$ and the total density $\tilde{n}(\tilde{x})$, respectively.

Now we come to the Bose-glass region, where we have from equation (2) in TF approximation ${\tilde{n}}_{0}(\tilde{x})=0$, and equation (3) reduces to $\tilde{n}(\tilde{x})=\tilde{q}(\tilde{x})$. Inserting this into equation (1), we get ${Q}_{0}(x)$ = ${\left(\tfrac{{\hslash }M}{D}\right)}^{1/3}$, which reduces equation (4) to:

$$\begin{eqnarray}&&\tilde{q}(\tilde{x})=\displaystyle \frac{1}{2}(3{\tilde{D}}^{2/3}+\tilde{\mu }-{\tilde{x}}^{2}).\end{eqnarray} \tag{ 10 }$$

We also need to write down the dimensionless equivalent of the normalisation condition (A6), which reads:

$$\begin{eqnarray}&&{\int }_{-{\tilde{R}}_{{\rm{TF2}}}}^{{\tilde{R}}_{{\rm{TF2}}}}\tilde{n}(\tilde{x}){\rm{d}}\tilde{x}=\displaystyle \frac{4}{3},\end{eqnarray} \tag{ 11 }$$

where ${\tilde{R}}_{{\rm{TF2}}}={R}_{{\rm{TF2}}}/{R}_{{\rm{TF}}}$ denotes the dimensionless cloud radius, and the total density $\tilde{n}(\tilde{x})$ in equation (11) is the combination of the total densities from both the superfluid region and the Bose-glass region.

Before considering any particular parameter value for our BEC system, we have first to justify using the TF approximation and determine its range of validity. To this end we rewrite equation (2) in the clean case, i.e., for D = 0, and divide it with $\bar{\mu }$. This yields:

$$\begin{eqnarray}&&\left[-1+\tilde{n}(\tilde{x})+{\tilde{x}}^{2}-{\left(\displaystyle \frac{\xi }{{R}_{{\rm{TF}}}}\right)}^{2}\displaystyle \frac{{\partial }^{2}}{\partial {\tilde{x}}^{2}}\right]\;\sqrt{\tilde{n}(\tilde{x})}=0.\end{eqnarray} \tag{ 12 }$$

Note that in the clean case the total density coincides with the condensate one. The TF approximation is only justified when the prefactor of the kinetic term ${\left(\tfrac{\xi }{{R}_{{\rm{TF}}}}\right)}^{2}$ is small enough, so the kinetic term can be neglected, which yields

$$\begin{eqnarray}&&\xi \ll {R}_{{\rm{TF}}}.\end{eqnarray} \tag{ 13 }$$

The corresponding TF results are presented and discussed in section 4. In order to asses their validity, in particular for large disorder strengths, however, we turn first to two other complementary approaches to the dirty BEC problem.

In this section we start with working out a numerical method that relies on solving the Gross–Pitaevskii equation for an ensemble of realisations of disorder landscape. Furthermore, we also introduce a variational approach, which is tailored to describe the numerical results analytically.

3.1. Numerical method

Now we perform a numerical study for the Bose-condensed gas in one dimension at zero temperature in a harmonic trapping potential $V(x)=\tfrac{1}{2}M{{\rm{\Omega }}}^{2}{x}^{2}$. Furthermore, we assume a Gaussian-distributed disorder potential U(x), which satisfies the conditions

$$\begin{eqnarray}&&\overline{U(x)}=0,\end{eqnarray} \tag{ 14 }$$

and

$$\begin{eqnarray}&&\bar{U(x)U({x}^{\prime })}=D(x-{x}^{\prime }),\end{eqnarray} \tag{ 15 }$$

where $D(x-{x}^{\prime })$ denotes the correlation function.

A one-dimensional BEC in the mean-field Hartree approximation is given by a generalised time-independent Gross–Pitaevskii equation for the condensate wave function $\psi (x)$:

$$\begin{eqnarray}&&\left[-\displaystyle \frac{{{\hslash }}^{2}}{2M}\displaystyle \frac{{\partial }^{2}}{\partial {x}^{2}}-\mu +U(x)+V(x)+g| \psi (x){| }^{2}\right]\psi (x)=0.\end{eqnarray} \tag{ 16 }$$

Equation (16) represents a stochastic nonlinear differential equation which can not be solved exactly, and, therefore, we apply a numerical approach. To this end we have first to generate the random potential U(x) before inserting it into equation (16), and then calculate the disorder average over many realisations of U(x).

Motivated by Fourier series, a simple ansatz for generating a random Gaussian function U(x) is performed as follows. The potential is written as a finite superposition of $\mathrm{sin}{kx}$ and $\mathrm{cos}{kx}$ terms with properly selected amplitudes A_n, B_n, and wave numbers k_n [51, 52]:

$$\begin{eqnarray}&&U(x)=\displaystyle \frac{1}{\sqrt{{\mathtt{N}}}}\displaystyle \sum _{n=0}^{{\mathtt{N}}-1}({A}_{n}\mathrm{cos}{k}_{n}x+{B}_{n}\mathrm{sin}{k}_{n}x),\end{eqnarray} \tag{ 17 }$$

where N denotes the number of terms, which should be large enough in order to obtain a good approximation for the random potential. Furthermore, we assume$\ {A}_{n}$ and$\ {B}_{n}\$to be mutually independent Gaussian random variables with zero mean, and variance equal to D(0):

$$\begin{eqnarray}&&\bar{{A}_{n}{B}_{n}}=0,\qquad \bar{{A}_{n}{A}_{m}}=\bar{{B}_{n}{B}_{m}}=D(0){\delta }_{{nm}}.\end{eqnarray} \tag{ 18 }$$

The wave numbers k_n are independent random variables, as well, selected from the probability distribution:

$$\begin{eqnarray}&&p({k}_{n})=\displaystyle \frac{S({k}_{n})}{{\int }_{-\infty }^{\infty }S({k}^{\prime }){{\rm{d}}{k}}^{\prime }},\end{eqnarray} \tag{ 19 }$$

where S(k) defines the spectral density as the Fourier transform of the correlation function:

$$\begin{eqnarray}&&S(k)={\int }_{-\infty }^{\infty }{\rm{d}}{x}\;{{\rm{e}}}^{-{\rm{i}}{k}{x}}D(x).\end{eqnarray} \tag{ 20 }$$

In the special case of the Gaussian-correlated disorder we have

$$\begin{eqnarray}&&D(x-{x}^{\prime })=\displaystyle \frac{D}{\sqrt{2\pi }\lambda }\;{{\rm{e}}}^{-\tfrac{{(x-{x}^{\prime })}^{2}}{2{\lambda }^{2}}},\end{eqnarray} \tag{ 21 }$$

where λ denotes the correlation length and D the disorder strength. The probability distribution (19) reads in this case:

$$\begin{eqnarray}&&p({k}_{n})=\displaystyle \frac{\lambda }{\sqrt{2\pi }}\;{{\rm{e}}}^{-\tfrac{{\lambda }^{2}{k}_{n}^{2}}{2}}.\end{eqnarray} \tag{ 22 }$$

Note that the analytical study in section 2 is done for δ-correlated disorder, but since it is impossible to treat the δ-correlated disorder numerically, we use the Gaussian-distributed disorder (21), which specialises to a δ-distributed one in the limit $\lambda \to 0$, i.e., ${\mathrm{lim}}_{\lambda \to 0}D(x)=D\delta (x).$

In order to numerically generate the correlation function (21) with sufficient accuracy, two numbers have to be appropriately large enough. The first one is the number N of terms in equation (17), the second one is the number M of realisations of the disorder potential, which are used to evaluate the disorder ensemble average (21). It can be shown analytically that the error in reproducing the correlation function (21) in the case ${\mathtt{N}}\to \infty$ is of the order of ${\mathtt{1}}/{\mathtt{N}}$ [52]. All Gaussian-correlated quantities are generated using the Box–Müller algorithm [53]. We insert the generated disorder potential (17) into the Gross–Pitaevskii equation (16), and then use a C computer program that solves the time-independent Gross–Pitaevskii equation in one space dimension in a harmonic trap using the imaginary-time propagation [54–58]. In this way we obtain the numerical solution of the ground-state wave function $\psi (x)$ of equation (16) for ${\mathtt{M}}=1000$ realisations of the disorder potential and ${\mathtt{N}}=\mathrm{100\; 00}$ terms in equation (17). To this end we use different values of the disorder strength D in order to cover the range from the weak to the intermediate disorder regime. We have chosen the disorder correlation length to be $\lambda =0.01\;l$, which is small enough in order to approach the case of δ-correlated disorder.

Performing disorder ensemble averages, we have access to the particle density $n(x)=\bar{\psi {(x)}^{2}}$, to the condensate density ${n}_{0}(x)={\overline{\psi (x)}}^{2}$, and to the Bose-glass order parameter $q(x)=n(x)-{n}_{0}(x)$. In order to compare the numerical results with the analytical ones obtained in section 2, we use the same rescaling parameters for all densities, coordinates, chemical potential, and disorder strength, as already explained below equation (8).

Before discussing the numerical results in detail, we show first one typical example in two graphs in figure 1, where the total density $\tilde{n}(\tilde{x})$ is plotted for two different values of the disorder strength (solid, green line), showing the original data for ${\mathtt{M}}=1000$ and ${\mathtt{N}}=10000$ terms in equation (17). We remark that the resulting density is fluctuating around a Gaussian-like curve. Comparing figure 1(a) with figure 1(b) we conclude that the fluctuations are increasing with the disorder strength. The origin of those fluctuations is that the ${\mathtt{M}}=1000$ realisations of the disorder potential for performing the disorder ensemble average are not sufficient to produce a smooth curve. One solution of this problem would be to increase the number M of the realisations of the disorder potential, which would need longer execution time, especially for the intermediate disorder regime, where the numerics has to be run for a larger spatial range. Another solution is to extract a continuous smooth curve that fits best to our data, as it is done in figure 1 (dotted-dashed, red line). This method is applied to all numerical densities in this paper. Furthermore, from the Gaussian fit in figure 1 (dotted, blue line), we remark that the original data of the total density approach a Gaussian form in the intermediate disorder regime much better than in the weak disorder regime. This can be explained with the argument that increasing the disorder reduces effectively the repulsive interaction between the particles [50] and, thus, approaches the case of non-interacting bosons, where the total density is given by a Gaussian.

Zoom In Zoom Out Reset image size

3.2. Variational method

Since the four self-consistency equations (1)–(4) are obtained by extremizing the free energy (A5), we can apply the variational method in the spirit of [59–66] to obtain approximate results. In order to be able to compare the variational results with the analytical and the numerical ones from section 2 and the previous subsection, respectively, we use the same rescaling parameters already introduced below equation (8) for all functions and parameters. To this end, we have to multiply (A5) with the factor $1/(\bar{\mu }\bar{n}{R}_{{\rm{TF}}})$ to obtain:

$$\begin{eqnarray}\tilde{{ \mathcal F }} & = & \displaystyle \int {\rm{d}}\tilde{x}\left[-{[\tilde{q}(\tilde{x})+{\tilde{n}}_{0}(\tilde{x})]}^{2}-\sqrt{{\tilde{n}}_{0}(\tilde{x})}\left\{\tilde{\mu }+{\left(\displaystyle \frac{\xi }{{R}_{{\rm{TF}}}}\right)}^{2}\displaystyle \frac{{\partial }^{2}}{\partial {\tilde{x}}^{2}}-2[\tilde{q}(\tilde{x})+{\tilde{n}}_{0}(\tilde{x})]-{\tilde{x}}^{2}+2\tilde{D}{\tilde{Q}}_{0}(\tilde{x})\right\}\sqrt{{\tilde{n}}_{0}(\tilde{x})}\right.\\ & & -\left.\displaystyle \frac{1}{2}{\tilde{n}}_{0}{(\tilde{x})}^{2}+2\tilde{D}{\tilde{Q}}_{0}(\tilde{x})[\tilde{q}(\tilde{x})+{\tilde{n}}_{0}(\tilde{x})]-2\tilde{D}\displaystyle \frac{\tilde{q}(\tilde{x})+{\tilde{n}}_{0}(\tilde{x})}{\sqrt{-\tilde{\mu }+2[\tilde{q}(\tilde{x})+{\tilde{n}}_{0}(\tilde{x})]+{\tilde{x}}^{2}-2\tilde{D}{\tilde{Q}}_{0}(\tilde{x})}}\right],\end{eqnarray} \tag{ 23 }$$

where $\tilde{{ \mathcal F }}={ \mathcal F }/(\bar{\mu }\bar{n}{R}_{{\rm{TF}}})$ denotes the dimensionless free energy and ${\tilde{Q}}_{0}(\tilde{x})=\sqrt{\tfrac{2\bar{\mu }}{M}}{Q}_{0}(x)$.

Motivated by the numerical results presented in figure 1, we suggest the three following ansätze for the condensate density ${\tilde{n}}_{0}(\tilde{x})$, the Bose-glass order parameter $\tilde{q}(\tilde{x})$, and the auxiliary function ${\tilde{Q}}_{0}(\tilde{x})$:

$$\begin{eqnarray}&&{\tilde{n}}_{0}(\tilde{x})=\alpha {{\rm{e}}}^{-\sigma {\tilde{x}}^{2}},\end{eqnarray} \tag{ 24 }$$

$$\begin{eqnarray}&&\tilde{q}(\tilde{x})+{\tilde{n}}_{0}(\tilde{x})=\gamma {{\rm{e}}}^{-\theta {\tilde{x}}^{2}},\end{eqnarray} \tag{ 25 }$$

$$\begin{eqnarray}&&{\tilde{Q}}_{0}(\tilde{x})=\displaystyle \frac{\tilde{q}(\tilde{x})+{\tilde{n}}_{0}(\tilde{x})}{\tilde{D}}-(\zeta +\eta {\tilde{x}}^{2}),\end{eqnarray} \tag{ 26 }$$

where α, σ, γ, θ, ζ, and η denote variational parameters. The parameters α and γ are proportional to the number of particles in the condensate and the total number of particles, while parameters σ and θ represent the width of the condensate density and the total density, respectively.

Inserting the ansätze (24)–(26) into the free energy (23) and performing the integral yields:

$$\begin{eqnarray}\tilde{{ \mathcal F }} & = & \sqrt{\pi }\left\{\displaystyle \frac{{\gamma }^{2}}{\sqrt{2\theta }}+\displaystyle \frac{\alpha }{2{\sigma }^{3/2}}-\displaystyle \frac{\alpha }{4\sqrt{\sigma }}(4\tilde{\mu }+\sqrt{2}\alpha )+{\left(\displaystyle \frac{\xi }{{R}_{{\rm{TF}}}}\right)}^{2}\displaystyle \frac{\alpha \sqrt{\sigma }}{2}+\tilde{D}\left(\displaystyle \frac{2\alpha \zeta }{\sqrt{\sigma }}+\displaystyle \frac{\alpha \eta }{{\sigma }^{3/2}}-\displaystyle \frac{\gamma (\eta +2\zeta \theta )}{{\theta }^{3/2}}\right)\right.\\ & & -\left.\displaystyle \frac{2\tilde{D}\gamma }{\sqrt{1+2\tilde{D}\eta }}{{\rm{e}}}^{\tfrac{2\tilde{D}\zeta -\tilde{\mu }}{2+4\tilde{D}\eta }\theta }{K}_{0}\left(\displaystyle \frac{2\tilde{D}\zeta -\tilde{\mu }}{2+4\tilde{D}\eta }\theta \right)\right\},\end{eqnarray} \tag{ 27 }$$

where ${K}_{0}(s)$ represents the modified Bessel function of the second kind.

The free energy (27) has now to be extremized with respect to the variational parameters α, σ, γ, θ, ζ, and η. Together with the thermodynamic condition $-\tfrac{\partial \tilde{{ \mathcal F }}}{\partial \tilde{\mu }}=\tfrac{4}{3}$, we have seven coupled equations for seven variables α, σ, γ, θ, ζ, η, and $\tilde{\mu }$ that we solve numerically. From all physical solutions we select the one with the smallest free energy (27), then we insert the resulting variational parameters α, σ, γ, and θ into the ansätze (24) and (25) in order to get the total density $\tilde{n}(\tilde{x})$, the condensate density ${\tilde{n}}_{0}(\tilde{x})$, and the Bose-glass order parameter $\tilde{q}(\tilde{x})$.

The results presented in this section correspond to a dirty BEC with $N={10}^{6}$ atoms of ⁸⁷Rb, with the s-wave scattering length $a=100\;{a}_{0}=5.29\;\mathrm{nm}$, where a₀ represents the Bohr radius. For the trap frequencies we use experimentally realistic parameters: the longitudinal frequency is chosen to be ${\rm{\Omega }}=2\pi \times 50\;\mathrm{Hz}$, and the radial one ${\omega }_{r}=2\pi \times 179\;\mathrm{Hz}$. For those parameters the longitudinal and the transversal oscillator lengths read $l=\sqrt{\tfrac{{\hslash }}{m{\rm{\Omega }}}}=1.52\;\mu {\rm{m}}$ and ${l}_{r}=\sqrt{\tfrac{{\hslash }}{m{\omega }_{r}}}=806.04\;\mathrm{nm}$, respectively, while the coherence length in the trap centre in the clean case turns out to be $\xi =45.6\;\mathrm{nm}$, and the Thomas–Fermi radius reads ${R}_{\mathrm{TF}}=50.9\;\mu {\rm{m}}$. Regarding the geometry of the system, we see that the transversal oscillator length is much larger than the scattering length, $a\ll {l}_{r}$, but still smaller than the longitudinal oscillator length, ${l}_{r}\lt l$, so we are indeed in the quasi one-dimensional regime [67, 68]. If we estimate the value of the dimensionless quantity ${\gamma }_{1{\rm{D}}}={E}_{\mathrm{int}}/{E}_{\mathrm{kin}}={Mg}/{{\hslash }}^{2}\bar{n}$ [69], which compares the interaction and the kinetic energy of the system, where the effective 1D interaction strength is given by $g=2{{\hslash }}^{2}{Ma}/{l}_{r}^{2}$ and $\bar{n}$ is defined in section 2, we get ${\gamma }_{1{\rm{D}}}=2\times {10}^{-7}\ll 1$. This clearly shows that, for the chosen parameters, our system is in the weakly interacting regime, and that we can describe it using the Hartree–Fock mean-field theory [69]. Furthermore, if we calculate the value of the dimensionless quantity ${\alpha }_{1{\rm{D}}}={Mgl}/{{\hslash }}^{2}=2{al}/{l}_{r}^{2}$, as defined in [69], which relates the effective 1D interaction strength g and the longitudinal trap frequency Ω, we obtain ${\alpha }_{1{\rm{D}}}=0.024\ll 1$. Since the condition $N{\alpha }_{1{\rm{D}}}\gg 1$ [69] is satisfied, we see that the system is close to the TF regime. This justifies our approach that, while using the TF approximation to obtain some analytical results and understand behaviour of the system in general, full numerical treatment is still necessary in order to describe the properties of 1D dirty bosons. In particular this will become obvious when we observe unphysical features of the TF approximation.

We first present results of the TF approximation and afterwards give in parallel numerical and variational results and compare them.

4.1. Thomas–Fermi results

In section 2 we have presented an analytical theory for the dirty boson problem in 1D. Using the above specified parameter values, we now solve equation (9) numerically. To this end we select from all its real solutions for ${\tilde{n}}_{0}(\tilde{x})$ the physical one, i.e., the one with the smallest energy, and denote the region where it is non-trivial as a superfluid region. The Bose-glass order parameter is here determined by equation (7). Then we combine the superfluid region solution with equation (10), describing the pure Bose-glass region, in which ${\tilde{n}}_{0}(\tilde{x})=0$ and $\tilde{n}(\tilde{x})=\tilde{q}(\tilde{x})$. After that we fix the chemical potential $\tilde{\mu }$ using the normalisation condition (11). The resulting densities are combined and plotted in figure 2.

Zoom In Zoom Out Reset image size

The cloud radius ${\tilde{R}}_{\mathrm{TF}2}$ for the system is determined by the condition that the total density vanishes, $\tilde{n}({\tilde{R}}_{\mathrm{TF}2}=0$, while the condensate radius ${\tilde{R}}_{\mathrm{TF}1}$ is maximal value of the coordinate $\tilde{x}$ for which equation (9) still has a solution, and separates the superfluid region ($| \tilde{x}| \leqslant {\tilde{R}}_{\mathrm{TF}1}$) from the Bose-glass one (${\tilde{R}}_{\mathrm{TF}1}\lt | \tilde{x}| \leqslant {\tilde{R}}_{\mathrm{TF}2}$). This is illustrated in figure 2 for the dimensionless disorder strength $\tilde{D}=0.016$, where we can see that the total density has a small jump at the condensate radius, while the condensate density exhibits a jump to zero. The Bose-glass order parameter $\tilde{q}(\tilde{x})$ has a double-bump structure, exhibits a jump at the condensate radius, and also vanishes at the cloud radius, by definition. In the Bose-glass region, the Bose-glass order parameter and the total density coincide. The jump exhibited by all three densities at the condensate radius is not a physical one, and is an artefact of the applied TF approximation.

To study the influence of the disorder on the BEC properties, we plot the resulting TF radii in figure 3(a) as a function of the disorder strength $\tilde{D}$. We see that both cloud and condensate radius coincide in the clean case, as expected. The condensate radius decreases with increasing disorder strength and vanishes at the critical value ${\tilde{D}}_{{\rm{c}}}=0.143$, which marks a quantum phase transition from the superfluid to the Bose-glass phase. This corresponds to the value ${\tilde{D}}_{{\rm{c}}}=0.333$, which was found in the non-perturbative approach of [37, 38] to be the critical disorder strength, where the Bose-glass phase becomes energetically unstable and goes over into the superfluid phase. On the other side, the cloud radius ${\tilde{R}}_{\mathrm{TF}2}$ increases with the disorder in the superfluid phase, but remains constant in the Bose-glass phase at the value ${\tilde{R}}_{\mathrm{TF}2}=1.256$. This means that beyond the critical disorder strength ${\tilde{D}}_{{\rm{c}}}$ the bosonic cloud is not extending anymore and has a maximal size. The same conclusion can be deduced from figure 3(b), where we depict the fractional number of condensed particles ${N}_{0}/N=\tfrac{3}{4}{\int }_{-{\tilde{R}}_{\mathrm{TF}1}}^{{\tilde{R}}_{\mathrm{TF}1}}{\tilde{n}}_{0}(\tilde{x}){\rm{d}}\tilde{x}$. Here ${N}_{0}/N$ equals to one in the clean case, i.e., all particles are in the condensate. Afterwards it decreases as the disorder strength $\tilde{D}$ increases, until it vanishes at ${\tilde{D}}_{{\rm{c}}}$, marking the end of the superfluid phase and the beginning of the Bose-glass phase. The fraction of the atoms in the disconnected local mini-condensates, $Q/N=\tfrac{3}{4}{\int }_{-{R}_{\mathrm{TF}2}}^{{R}_{\mathrm{TF}2}}\tilde{q}(\tilde{x}){\rm{d}}\tilde{x}$, behaves conversely. It increases with the increasing disorder until reaching the maximal value of one at ${\tilde{D}}_{{\rm{c}}}$, then it remains equal to one in the Bose-glass phase since all particles are stuck in the local minima of the disorder potential.

Zoom In Zoom Out Reset image size

4.2. Numerical and variational results

Now we turn to numerical and variational results obtained using the methods presented in section 3. Figure 4 presents in parallel numerically and variationally obtained densities $\tilde{n}(\tilde{x})$, ${\tilde{n}}_{0}(\tilde{x})$, and $\tilde{q}(\tilde{x})$ for various values of the disorder strength $\tilde{D}$. The first notable difference compared to TF results is that the condensate and the cloud radius are not clearly defined, since the densities do not vanish at a well-defined point, but gradually converge to zero at the respective borders. Therefore, we define the corresponding radii, which, for simplicity, we again denote by ${\tilde{R}}_{\mathrm{TF}1}$ and ${\tilde{R}}_{\mathrm{TF}2}$, by the conditions ${\tilde{n}}_{0}({\tilde{R}}_{\mathrm{TF}1})=\varepsilon$ and $\tilde{n}({\tilde{R}}_{\mathrm{TF}2})=\varepsilon$, where $\varepsilon ={10}^{-4}$ represents a conveniently chosen small number. From figures 4(a) and (d) we see that the cloud radius increases with increasing disorder strength, while the maximal density at the trap centre decreases. Figures 4(b) and (e) show the same type of behaviour for the condensate, and generally reveal good agreement between the full numerical and variational results, as for the total density.

Zoom In Zoom Out Reset image size

The numerically and variationally obtained values of the Bose-glass order parameter are also plotted for different values of $\tilde{D}$ in figures 4(c) and (f), respectively. While the variational results have similar form as the total particle density and the condensate density, due to the assumed functional dependence in equations (24) and (25), the numerical results show completely different behaviour. In the weak disorder regime, the numerically calculated order parameter $\tilde{q}(\tilde{x})$ reveals a double bump structure and is maximal at the border of the condensate, while in the intermediate disorder regime it resembles a Gaussian-like form. This redistribution takes place, according to figure 4(c), at a disorder strength value between $\tilde{D}=0.151$ and $\tilde{D}=0.268$. Thus, the main difference between the weak and the intermediate disorder regime is that the local condensates concentrate at the border of the condensate in the former case, but sit in the trap centre in the latter case. Despite this marked difference, using either method yields that the width as well as the maximum of the Bose-glass order parameter increase with the disorder strength.

In order to obtain further information on the behaviour of the system, we plot in figures 5(a) and (c) the numerical and variational fractional number of condensed particles ${N}_{0}/N$ and fractional number of particles in the disconnected mini-condensates Q/N, as functions of the disorder strength, respectively. The condensed fraction ${N}_{0}/N$ decreases with the disorder strength and, conversely, Q/N increases, meaning that more and more particles are leaving the condensate with increasing $\tilde{D}$. Unfortunately, the employed numerical algorithm breaks down for larger values of the disorder strength, and one would have to use other approaches in this case. Starting from the disorder strength value $\tilde{D}=0.393$, the variational equations turn out to have only complex solutions, so we cannot extract further information about our system for higher disorder strengths using this approach. Therefore, we focus on the regimes of weak to moderate disorder.

Zoom In Zoom Out Reset image size

Numerically and variationally calculated cloud radius ${\tilde{R}}_{\mathrm{TF}2}$ and the condensate radius ${\tilde{R}}_{\mathrm{TF}1}$, as defined in this subsection, are plotted in figures 5(b) and (d) as functions of the disorder strength, respectively. Both radii are almost identical in the weak disorder regime, and afterwards both increase linearly with the disorder strength in the moderate disorder regime. Since the applied method breaks down for larger disorder strengths, we cannot determine if a quantum phase transition exists, and this question remains still open.

4.3. Comparison

Now we compare the physical quantities obtained via the three different methods: the TF approximation, the numerical method, and the variational method. For the small disorder strength $\tilde{D}=0.016$, the three total densities $\tilde{n}(\tilde{x})$ in figure 6(a) agree qualitatively well, but quantitatively the TF-approximated function $\tilde{n}(\tilde{x})$ is a better approximation for the numerical total density, especially in the centre of the bosonic cloud, where the variational result does not agree well with the numerical one. The same remark can be made for the condensate density ${\tilde{n}}_{0}(\tilde{x})$ in figure 6(b). For the Bose-glass order parameter $\tilde{q}(\tilde{x})$ in figure 7(c) the double-bump structure, which exists in both numerical and TF-approximated results, is missing in the variational result, which has just a bell form, as assumed by the variational ansatz. This makes again the TF-approximated Bose-glass order parameter $\tilde{q}(\tilde{x})$ a better approximation for the numerical one.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

For the moderate disorder strength $\tilde{D}=0.386$, the TF-approximated total density $\tilde{n}(\tilde{x})$ in figure 6(d) is also a better approximation for the numerical one in the centre of the bosonic cloud, while at the trap borders the variational approximation wins. According to the TF results shown in figure 3, at the disorder strength value $\tilde{D}=0.386$ we are already in the Bose-glass phase, thus, the TF-approximated condensate density ${\tilde{n}}_{0}(\tilde{x})$ in figure 6(e) vanishes. This is not the case for both the numerical and the variational condensate densities, which are compatible and match quite well at trap borders. The variational Bose-glass order parameter $\tilde{q}(\tilde{x})$ in figure 6(f) also agrees well with the numerical one and both have the same bell shape, while the TF-approximated Bose-glass order parameter has a significant deviation, which is expected since the TF approximation breaks down in the moderate disorder regime.

In figures 7(a) and (b) we see that the variational and the numerical condensate radius ${\tilde{R}}_{\mathrm{TF}1}$ and cloud radius ${\tilde{R}}_{\mathrm{TF}2}$ have the same behaviour, namely both increase with the disorder strength. This is in stark contrast to the TF result, where the condensate radius is found to decrease with $\tilde{D}$, evenatually leading to a quantum phase transition at ${\tilde{D}}_{{\rm{c}}}=0.143$. Such a quantum phase transition is predicted only in TF approximation, which is known to fail for moderate disorder strengths.

The question that still remains to be answered concerns the possible existence of the quantum phase transition from the superfluid to the Bose-glass phase for large disorder strengths. According to [38], the disorder has to energetically overcome the interaction in order to yield such a quantum phase transition. However, numerical and variational results displayed in figure 7(a) suggest that this is not the case neither in the weak nor in the moderate disorder regime. One would have to investigate much stronger disorder regime, employing a different set of methods, in order to be able to detect a possible quantum phase transition, which is beyond the scope of this paper.

From the discussion in the previous section, we conclude that the TF approximation yields good results for the quasi-one-dimensional dirty bosons in the weak disorder regime, which agree well with the numerical ones, especially in the centre of the bosonic cloud, where the kinetic energy can be neglected. However, this approximation breaks down in the moderate disorder regime, and is unable to describe the dirty BEC system properly. The origin of this failure is the fact that the condition (13) is not fulfilled in the moderate disorder regime. The coherence length becomes significantly larger as we increase the disorder strength, especially at the border of the bosonic cloud, and when it becomes of the order of the Thomas–Fermi radius, the TF approximation breaks down. Furthermore, quantum fluctuations are more prominent in lower dimensions, which also restricts the validity range of the TF approximation. On the other side, the variational method with the ansätze (24)–(26) turns out to be a good approximation to describe the dirty BEC system in the moderate disorder regime and works there better than in the weak disorder regime, especially at the cloud border, where the Bose-glass region is located. This is due to the fact that a stronger disorder reduces significantly the repulsive interaction between the particles and approaches the case of non-interacting bosons, where the densities are Gaussian-like, as in our variational ansätze. Although the variational method breaks down for larger disorder strengths, it still provides results in an important range of disorder strengths. The combination of applying the TF approximation for the weak disorder together with the variational method for the moderate disorder covers a significant range of disorder strengths. With this we can analytically describe the redistribution of the local disconnected mini-condensates from the edge of the atomic cloud to the trap centre for increasing disorder strengths, as obtained from detailed numerical simulations. We expect that all these results are useful for a quantitative analysis of experiments for dirty bosons in quasi-one-dimensional harmonic traps. The problem of the large disorder strengths still persists with the current approach and remains to be addressed by other methods.

The authors gratefully acknowledge financial support from German Academic and Exchange Service (DAAD), the Ministry of Education, Science, and Technological Development of the Republic of Serbia under projects ON171017, NAI-DBEC, and IBEC, and German Research Foundation (DFG) via the Collaborative Research Center SFB/TR49 'Condensed Matter Systems with Variable Many-Body Interactions'. Numerical simulations were run on the PARADOX supercomputing facility at the Scientific Computing Laboratory of the Institute of Physics Belgrade.

Here we briefly summarise the main result of [50], which relies on deriving a Hartree–Fock approximation for the free energy of one-dimensional harmonically trapped dirty bosons. The starting point is the functional integral for the grand-canonical partition function

$$\begin{eqnarray}&&{ \mathcal Z }=\oint { \mathcal D }{\psi }^{\ast }\oint { \mathcal D }\psi {{\rm{e}}}^{-{ \mathcal A }[{\psi }^{\ast },\psi ]/{\hslash }},\end{eqnarray} \tag{ A1 }$$

where the integration is performed over all Bose fields ${\psi }^{\ast }(x,\tau ),\psi (x,\tau )$, which are periodic in imaginary time τ, i.e., $\psi (x,\tau )=\psi (x,\tau +{\hslash }\beta )$. The Euclidean action is given in standard notation by

$$\begin{eqnarray}{ \mathcal A }[{\psi }^{\ast },\psi ] & = & {\displaystyle \int }_{0}^{{\hslash }\beta }{\rm{d}}\tau \displaystyle \int {\rm{d}}{x}\left[{\psi }^{\ast }(x,\tau )\left\{{\hslash }\displaystyle \frac{\partial }{\partial \tau }-\displaystyle \frac{{{\hslash }}^{2}}{2M}{\rm{\Delta }}+V(x)+U(x)-\mu \right\}\psi (x,\tau )\right.\\ & & +\left.\displaystyle \frac{1}{2}\displaystyle \int {{\rm{d}}{x}}^{\prime }{\psi }^{\ast }(x,\tau )\psi (x,\tau ){V}^{({\rm{int}})}(x-{x}^{\prime }){\psi }^{\ast }({x}^{\prime },\tau )\psi ({x}^{\prime },\tau )\right],\end{eqnarray} \tag{ A2 }$$

where $V(x)=M{{\rm{\Omega }}}^{2}{x}^{2}/2$ denotes the harmonic trap with the trap frequency Ω, M the particle mass, μ the chemical potential, and ${V}^{({\rm{int}})}(x-{x}^{\prime })=g\delta (x-{x}^{\prime })$ the contact interaction potential. The interaction coupling strength $g=2a{\hslash }{\omega }_{r}$ depends on the s-wave scattering length a, which has to be positive in order to obtain a stable BEC, and the transversal trap frequency ${\omega }_{r}$. Note that the latter has to be large enough, i.e., ${\omega }_{r}\gg {\rm{\Omega }}$, in order to ensure a quasi-one-dimensional setup [67, 68].

We assume for the disorder potential U(x) that it is homogeneous after performing the disorder ensemble average (denoted by $\bar{\cdot }$) over all possible realisations. Thus, the expectation value of the disorder potential can be set to vanish without loss of generality, as defined by equation (14). The disorder correlation function is given by equation (15), where we assume $D(x-{x}^{\prime })=D\;\delta (x-{x}^{\prime })$, and D denotes the disorder strength.

Within the Hartree–Fock mean-field approximation with the replica method, [50] obtains self-consistency equations, which determine the particle density n(x) as well as the order parameter of the superfluid ${n}_{0}(x)$, which represents the condensate density, and the order parameter of the Bose-glass phase q(x), which stands for the density of the particles being condensed in the respective minima of the disorder potential.

More precisely, the two order parameters ${n}_{0}(x)$ and q(x) of our mean-field theory at T = 0 are defined by following the notion of spin-glass theory [70–72]. On the one hand, the off-diagonal long-range limit of the two-point correlation function defines the condensate density [73],

$$\begin{eqnarray}&&\underset{| x-{x}^{\prime }| \to \infty }{\mathrm{lim}}\bar{\langle \psi (x,\tau ){\psi }^{\ast }({x}^{\prime },\tau )\rangle }=\sqrt{{n}_{0}(x){n}_{0}({x}^{\prime })}.\end{eqnarray} \tag{ A3 }$$

On the other hand, the Bose-glass order parameter q(x) was introduced in [39] in close analogy to the Edward-Anderson order parameter of spin-glasses [74] by the off-diagonal long-range limit of the four-point correlation function,

$$\begin{eqnarray}&&\underset{| x-{x}^{\prime }| \to \infty }{\mathrm{lim}}\bar{{| \langle \psi (x,\tau ){\psi }^{\ast }({x}^{\prime },\tau )\rangle | }^{2}}=[{n}_{0}(x)+q(x)][{n}_{0}({x}^{\prime })+q({x}^{\prime })].\end{eqnarray} \tag{ A4 }$$

Note that the validity of equations (A3) and (A4) within the Hartree–Fock approximation was analysed in detail in [50]. There it is also shown that two- and four-point correlations decay exponentially on a length scale which can be physically interpreted as the localisation length named after Larkin [37, 75]. A similar exponential decay of the one-body density matrix in the presence of disorder was found in [1, 31, 76], although in the quasi-condensed phase a power-lay decay is expected as in the spatially uniform gas [77, 78].

In the one-dimensional case and at T = 0, the mean-field Hartree–Fock theory with the help of the replica method leads to the free energy [50]:

$$\begin{eqnarray}&&{ \mathcal F } \displaystyle \int {\rm{d}}{x}\quad \left[{\displaystyle \frac{q(x)+{n}_{0}(x)}{\sqrt{-\mu +2g[q(x)+{n}_{0}(x)]+\tfrac{1}{2}M{{\rm{\Omega }}}^{2}{x}^{2}-\tfrac{D}{{\hslash }}{Q}_{0}(x)}}\phantom{\rule{-20.5em}{0ex}}}\right.-g{[q(x)+{n}_{0}(x)]}^{2}-\sqrt{{n}_{0}(x)}\left\{\mu +\displaystyle \frac{{{\hslash }}^{2}}{2M}\displaystyle \frac{{\partial }^{2}}{\partial {x}^{2}}-2g[q(x)+{n}_{0}(x)]-\displaystyle \frac{1}{2}M{{\rm{\Omega }}}^{2}{x}^{2}\right.\\ & & +\left.\displaystyle \frac{D}{{\hslash }}{Q}_{0}(x)\right\}\sqrt{{n}_{0}(x)}-\displaystyle \frac{g}{2}{n}_{0}{(x)}^{2}+\displaystyle \frac{D}{{\hslash }}{Q}_{0}(x)[q(x)+{n}_{0}(x)]\\ & & -\left.\displaystyle \frac{D}{{\hslash }}\sqrt{\displaystyle \frac{M}{2}}\displaystyle \frac{q(x)+{n}_{0}(x)}{\sqrt{-\mu +2g[q(x)+{n}_{0}(x)]+\tfrac{1}{2}M{{\rm{\Omega }}}^{2}{x}^{2}-\tfrac{D}{{\hslash }}{Q}_{0}(x)}}\right].\end{eqnarray} \tag{ A5 }$$

Here ${Q}_{0}(x)$ represents an auxiliary function within the Hartree–Fock theory. From the thermodynamic relation $-\tfrac{\partial { \mathcal F }}{\partial \mu }=N$ we obtain

$$\begin{eqnarray}&&{\int }_{-\infty }^{\infty }n(x){\rm{d}}{x}=N,\end{eqnarray} \tag{ A6 }$$

which defines the particle density n(x). Extremizing the free energy (A5) with respect to ${n}_{0}(x)$, q(x), and ${Q}_{0}(x)$ yields, together with equation (A6), the self-consistency equations (1)–(4).

The simplest case to discuss for dirty bosons is the homogeneous one, where $V(x)=0$. Since all densities are spatially constant in the homogeneous case, we drop in this section the x dependency of all densities. With this, equations (1)–(4) reduce, after eliminating Q₀, to:

$$\begin{eqnarray}&&q=\displaystyle \frac{D}{{\hslash }M}\;\displaystyle \frac{{n}_{0}}{{\left(\tfrac{2{{gn}}_{0}}{M}\right)}^{\tfrac{3}{2}}-\tfrac{D}{{\hslash }M}},\end{eqnarray} \tag{ B1 }$$

$$\begin{eqnarray}&&{{gn}}_{0}=-\mu +2{gn}-\displaystyle \frac{D}{{\hslash }}\sqrt{\displaystyle \frac{M}{2{{gn}}_{0}}},\end{eqnarray} \tag{ B2 }$$

$$\begin{eqnarray}&&n=q+{n}_{0}.\end{eqnarray} \tag{ B3 }$$

From equations (B1) and (B3) we get an algebraic fifth-order equation for the condensate fraction ${n}_{0}/n$:

$$\begin{eqnarray}&&{\left(\displaystyle \frac{{n}_{0}}{n}\right)}^{5/2}-{\left(\displaystyle \frac{{n}_{0}}{n}\right)}^{3/2}+\bar{D}=0,\end{eqnarray} \tag{ B4 }$$

where $\bar{D}=\tfrac{{\xi }^{3}}{{{ \mathcal L }}^{3}}$ denotes the dimensionless disorder strength, $\xi =\tfrac{{\hslash }}{\sqrt{2{Mgn}}}$ the coherence length, and ${ \mathcal L }={\left(\tfrac{{{\hslash }}^{4}}{{M}^{2}D}\right)}^{1/3}$ the Larkin length [37, 75].

Figure 8 shows how condensate fraction ${n}_{0}/n$ decreases with increasing the disorder strength $\bar{D}$ according to equation (B4). Thus, the Hartree–Fock mean-field theory predicts a first-order quantum phase transition from the superfluid phase to the Bose-glass phase at the critical value ${\bar{D}}_{{\rm{c}}}=\tfrac{6}{25}\sqrt{\tfrac{3}{5}}\simeq 0.185$. This corresponds to the value ${\bar{D}}_{{\rm{c}}}=1$, that was found in the non-perturbative approach of [37, 38], which investigate at which disorder strength the Bose-glass phase becomes energetically unstable and goes over into the superfluid phase. Therefore, we expect that a quantum phase transition will also appear in the trapped case within the Thomas–Fermi approximation.

Zoom In Zoom Out Reset image size

Now we check whether our results are compatible with the Huang–Meng theory [20–31], where the Bose-glass order parameter of a homogeneous dilute Bose gas at zero temperature in case of weak disorder regime is deduced within the seminal Bogoliubov theory. The Bose-glass order parameter in one dimension is according to the Huang–Meng theory proportional to the disorder strength, which yields in dimensionless form

$$\begin{eqnarray}&&\displaystyle \frac{{q}_{{\rm{HM}}}}{n}=\displaystyle \frac{\bar{D}}{{2}^{3/2}}.\end{eqnarray} \tag{ B5 }$$

In our Hartree–Fock mean-field theory the Bose-glass order parameter in case of weak disorder strength turns out to be:

$$\begin{eqnarray}&&\displaystyle \frac{{q}_{w}}{n}=\bar{D}.\end{eqnarray} \tag{ B6 }$$

Thus, from equation (B6) we conclude that our result agrees qualitatively with the Huang–Meng theory. But quantitatively the comparison of equations (B5) and (B6) reveals that a factor of ${2}^{3/2}$ is missing in our result (B6). This is due to the fact that the Hartree–Fock theory does not contain the Bogoliubov channel, which is included in the Huang–Meng theory.