Many-body localization in continuum systems: two-dimensional bosons

TL;DR

This paper shows that two-dimensional weakly interacting bosons can remain localized in disordered systems at low temperatures, with stability depending on energy distribution and disorder strength, relevant for cold atom experiments.

Contribution

It demonstrates the stability of many-body localization in 2D bosonic systems and identifies conditions under which localization persists at finite temperatures.

Findings

Localization stable at low temperatures in 2D bosons

High-energy particles can destroy localization above a critical temperature

Truncating high-energy states can maintain localization at any temperature

Abstract

We demonstrate that many-body localization of two-dimensional weakly interacting bosons in disorder remains stable in the thermodynamic limit at sufficiently low temperatures. Highly energetic particles destroy the localized state only above a critical temperature, which increases with the strength of the disorder. If the particle distribution is truncated at high energies, as it does for cold atom systems, the localization can be stable at any temperature.

Figures

Equations

$$P_{\alpha }={\sum _{{\alpha }^{\prime }\beta {\beta }^{\prime }}}^{\prime }\frac{⟨{\alpha }^{\prime },{\beta }^{\prime }|H_{\text{int}}|\alpha ,\beta ⟩}{{\Delta }_{\alpha \beta }^{{\alpha }^{\prime }{\beta }^{\prime }}}\sim C,$$

$${\eta }_{hc}=\frac{V_{hc}{N}_{\text{eff}}}{{\Delta }_h}.$$

$$\hat{H}={\hat{H}}_0+{\hat{H}}_{int},$$

$${\hat{H}}_0=\int d^{2}r\left(-{\hat{\Psi }}^{†}(\text{r})\frac{{\hslash }^2}{2m}{\nabla }^2\hat{\Psi }(\text{r})+{\hat{\Psi }}^{†}(\text{r})U(\text{r})\hat{\Psi }(\text{r})\right),$$

$${\hat{H}}_{int}=g\int d^{2}r{\hat{\Psi }}^{†}(\text{r}){\hat{\Psi }}^{†}(\text{r})\hat{\Psi }(\text{r})\hat{\Psi }(\text{r}).$$

$${ϵ}_{\ast }=\frac{m{U}_0^2{\sigma }^2}{\pi {\hslash }^2};{\zeta }_{\ast }=\sqrt{\frac{2e^2}{\pi }}\frac{{\hslash }^2}{m{U}_0\sigma }.$$

$$\zeta (ϵ)=\frac{{\zeta }_{\ast }}{e}\sqrt{\frac{ϵ}{{ϵ}_{\ast }}}\exp \left(\frac{ϵ}{{ϵ}_{\ast }}\right),$$

$$\frac{ng}{T_d}=\frac{mg}{2\pi {\hslash }^2}\ll 1.$$

$${ϵ}_{\ast }\ll T_d.$$

$${\eta }_{hc}\sim \frac{n^{2}g{\zeta }_{\ast }^2}{e^2{T}_d}\frac{ϵ}{{ϵ}_{\ast }}\exp \left[-ϵ\left(\frac{1}{T}-\frac{2}{{ϵ}_{\ast }}\right)\right].$$

$$P_{\alpha }={\sum _{{\alpha }^{\prime }\beta {\beta }^{\prime }}}^{\prime }\frac{⟨{\alpha }^{\prime },{\beta }^{\prime }|H_{int}|\alpha ,\beta ⟩}{{\Delta }_{\alpha \beta }^{{\alpha }^{\prime }{\beta }^{\prime }}}={\sum _{{\alpha }^{\prime }\beta {\beta }^{\prime }}}^{\prime }\frac{U_{\alpha \beta }^{{\alpha }^{\prime }{\beta }^{\prime }}{N}_{\alpha \beta }^{{\alpha }^{\prime }{\beta }^{\prime }}}{{\Delta }_{\alpha \beta }^{{\alpha }^{\prime }{\beta }^{\prime }}}\sim C,$$

$$N_{\alpha \beta }^{{\alpha }^{\prime }{\beta }^{\prime }}=\overline{\sqrt{|N_{\beta }(1+N_{{\alpha }^{\prime }})(1+N_{{\beta }^{\prime }})-N_{{\alpha }^{\prime }}{N}_{{\beta }^{\prime }}(1+N_{\beta })|}}$$

$$U_{\alpha \beta }^{{\alpha }^{\prime }{\beta }^{\prime }}=g\int {\psi }_{\alpha }(\text{r}){\psi }_{{\alpha }^{\prime }}(\text{r}){\psi }_{\beta }(\text{r}){\psi }_{{\beta }^{\prime }}(\text{r})d^{2}r$$

$${\Delta }_{\alpha \beta }^{{\alpha }^{\prime }{\beta }^{\prime }}=|{ϵ}_{\alpha }+{ϵ}_{\beta }-{ϵ}_{{\alpha }^{\prime }}-{ϵ}_{{\beta }^{\prime }}|.$$

$${\psi }_{\alpha }(\text{r})=\{\begin{array}{ll}{\zeta }_{\alpha }^{-1}& |\text{r}-{\text{r}}_{\alpha }|<{\zeta }_{\alpha }/2.\\ 0& \text{otherwise.}\end{array}$$

$$U_{\alpha \beta }^{{\alpha }^{\prime }{\beta }^{\prime }}\approx g\frac{\min ({\zeta }_{\alpha }^2,{\zeta }_{{\alpha }^{\prime }}^2,{\zeta }_{\beta }^2,{\zeta }_{{\beta }^{\prime }}^2)}{{\zeta }_{\alpha }{\zeta }_{{\alpha }^{\prime }}{\zeta }_{\beta }{\zeta }_{{\beta }^{\prime }}}$$

$$g_c({ϵ}_{\alpha })=C(\sum _{\beta ,{\alpha }^{\prime },{\beta }^{\prime }>\alpha }\frac{N_{\alpha \beta }^{{\alpha }^{\prime }{\beta }^{\prime }}}{{\Delta }_{\alpha \beta }^{{\alpha }^{\prime }{\beta }^{\prime }}}\frac{{\zeta }_{\alpha }}{{\zeta }_{{\alpha }^{\prime }}{\zeta }_{\beta }{\zeta }_{{\beta }^{\prime }}}+\sum _{{\alpha }^{\prime }<\alpha ,\beta ,{\beta }^{\prime }}\frac{N_{\alpha \beta }^{{\alpha }^{\prime }{\beta }^{\prime }}}{{\Delta }_{\alpha \beta }^{{\alpha }^{\prime }{\beta }^{\prime }}}\frac{{\zeta }_{{\alpha }^{\prime }}}{{\zeta }_{\alpha }{\zeta }_{\beta }{\zeta }_{{\beta }^{\prime }}}$$

$$+\sum _{\beta <\alpha ,{\alpha }^{\prime },{\beta }^{\prime }}\frac{N_{\alpha \beta }^{{\alpha }^{\prime }{\beta }^{\prime }}}{{\Delta }_{\alpha \beta }^{{\alpha }^{\prime }{\beta }^{\prime }}}\frac{{\zeta }_{\beta }}{{\zeta }_{{\alpha }^{\prime }}{\zeta }_{\alpha }{\zeta }_{{\beta }^{\prime }}}+{\sum _{{\beta }^{\prime }<\alpha ,{\alpha }^{\prime },\beta }\frac{N_{\alpha \beta }^{{\alpha }^{\prime }{\beta }^{\prime }}}{{\Delta }_{\alpha \beta }^{{\alpha }^{\prime }{\beta }^{\prime }}}\frac{{\zeta }_{{\beta }^{\prime }}}{{\zeta }_{{\alpha }^{\prime }}{\zeta }_{\beta }{\zeta }_{\alpha }})}^{-1}.$$

$$g<g_c=\min \{g({ϵ}_{\alpha })\}.$$

$$\min \{{\Delta }_{\alpha \beta \min }^{{\alpha }^{\prime }{\beta }^{\prime }}\}=\min \{{\delta }_{{\alpha }^{\prime }},{\delta }_{\beta },{\delta }_{{\beta }^{\prime }}\},$$

$$n={\int }_{-{ϵ}_{\ast }}^{\infty }{\rho }_0{\overline{N}}_{ϵ}dϵ.$$

$$E(\{N_{\alpha }\})=\sum _{\alpha }({ϵ}_{\alpha }{N}_{\alpha }+g{N}_{\alpha }(N_{\alpha }-1)/2{\zeta }_{\alpha }^2).$$

$$Z=\prod _{\alpha }{Z}_{\alpha },$$

$$Z_{\alpha }=\sum _{n=0}^{\infty }\exp (-(({ϵ}_{\alpha }-\mu )n+gn(n-1)/2{\zeta }_{\alpha }^2)/T).$$

$$Z_{\alpha }\approx \frac{1}{1-\exp (-({ϵ}_{\alpha }-\mu +g{\overline{N}}_{\alpha }/{\zeta }_{\alpha }^2)/T)}.$$

$${\overline{N}}_{ϵ}=T\frac{\partial \ln Z}{\partial \mu }\approx {\left[\exp \left(\frac{ϵ-\mu +{\overline{N}}_{ϵ}g/{\zeta }^2(ϵ)}{T}\right)-1\right]}^{-1}.$$

$${\overline{N}}_{ϵ}=\frac{{\zeta }^2(ϵ)}{2g}\left(\mu -ϵ+\sqrt{(\mu -ϵ{)}^2+\frac{4Tg}{{\zeta }^2(ϵ)}}\right).$$

$${\overline{N}}_{ϵ}=\frac{{\zeta }^2(ϵ)(\mu -ϵ)}{g}\theta (\mu -ϵ),$$

$${ϵ}_{\ast }(0{)}^{MBL}=0.87ng.$$

$${\overline{N}}_{ϵ}=\{\begin{array}{ll}\frac{{\zeta }^2(ϵ)}{2g}\left(\mu -ϵ+\sqrt{(ϵ-\mu {)}^2+\frac{4Tg}{{\zeta }^2(ϵ)}}\right)& ;ϵ<\mu \\ e^{-(ϵ-\mu )/T}& ;\mu <ϵ.\end{array}$$

$$N_{\alpha \beta }^{{\alpha }^{\prime }{\beta }^{\prime }}=\overline{\sqrt{|N_{\beta }(1+N_{{\alpha }^{\prime }})(1+N_{{\beta }^{\prime }})-N_{{\alpha }^{\prime }}{N}_{{\beta }^{\prime }}(1+N_{\beta })|}}.$$

$$N_{\alpha \beta }^{{\alpha }^{\prime }{\beta }^{\prime }}=\sqrt{|{\overline{N}}_{\beta }(1+{\overline{N}}_{{\alpha }^{\prime }})(1+{\overline{N}}_{{\beta }^{\prime }})-{\overline{N}}_{{\alpha }^{\prime }}{\overline{N}}_{{\beta }^{\prime }}(1+{\overline{N}}_{\beta })|}=\sqrt{|{\overline{N}}_{\beta }{\overline{N}}_{{\alpha }^{\prime }}+{\overline{N}}_{\beta }{\overline{N}}_{{\beta }^{\prime }}+{\overline{N}}_{\beta }-{\overline{N}}_{{\alpha }^{\prime }}{\overline{N}}_{{\beta }^{\prime }}|};{\overline{N}}_{\beta },{\overline{N}}_{{\beta }^{\prime }},{\overline{N}}_{{\alpha }^{\prime }}\gg 1,$$

Full text

After several decades since its original formulation, Anderson localization of quantum particles in disorder Anderson continues to attract a lot of interest 50years . One of the key questions is the effect of interactions on the localization Fleishman ; Altshuler . It was established that localization may persist when interactions are present Basko ; Gornyi1 , and the study of this phenomenon, known as many-body localization (MBL), is rapidly progressing ReviewHuse ; ReviewAbanin . A substantial amount of work, with a major role played by numerics, has been performed on lattice systems. The existence of the MBL phase has been proven in a rather general case of one-dimensional (1D) spin chains Imbrie , and a growing number of experimental observations provided further relevant evidence Tanzi ; Schreiber ; Smith ; Bordia ; Choi ; Bordia2 ; Rubio .

In continuum systems, there is no bound on the bandwidth as in lattice systems. The seminal work on the problem of MBL in continuum systems demonstrated that interacting particles can undergo many-body localization-delocalization transition (MBLDT), i.e. the transition from insulator to fluid state, and used an energy-independent single-particle localization length Basko ; Gornyi1 . This condition was relaxed in a number of subsequent works, which took into account the growth of the localization length with energy Aleiner ; Michal ; Michal2 ; Nandkishore ; Gornyi ; Bertoli . The question is whether the growth of the localization length eventually destabilizes the MBL phase, or whether the contribution of highly energetic states is hindered by the decrease in the thermal distribution function.

In the recent work Bertoli , it was shown that two-dimensional (2D) disordered bosons may display a finite-temperature insulator due to truncation of the energy distribution function at high energies, a generic phenomenon in cold atom systems walraven ; ketterle . The truncation ensures the survival of the insulating phase, because the energy $\epsilon$ of a localized single-particle state does not grow unbounded: hot particles quickly escape from the system. In the thermodynamic limit, however, the exponential increase of the localization length leads to the disappearance of the insulating phase above a critical temperature $T_{c}$. A similar conclusion was also obtained in Ref. Nandkishore through another approach. An opposite situation for the thermodynamic limit is found in Ref. Gornyi , where many-body localization is claimed to be unstable in any continuum system, irrespective of dimensionality, with a notable exception of one-dimensional Gaussian (white noise) disorder.

In this paper, we discuss the stability of many-body localization in continuum two-dimensional (2D) systems. The 2D case shows a stronger (exponential) growth of the localization length with energy compared to one dimension. The nature of quantum statistics is not crucial for this problem, and we consider disordered bosons for convenience. First, we review the arguments of Refs. Nandkishore ; Gornyi ; Bertoli , in which different criteria characterizing the MBLDT in continuum systems are employed. Then, building and improving on the results of Ref. Bertoli , we rule out the arguments of Gornyi on the absence of MBL in continuum systems of ultracold particles.

In general, the stability of the MBL phase is controlled by a parameter accounting for an increase of the phase space available for the transition when raising the temperature $T$. The key point is to compare the matrix element of interaction to the accessible level spacing Basko . When this ratio exceeds a model-dependent value of order unity, then delocalization takes place. The parameter controlling many-body delocalization was derived in Ref. Bertoli on the basis of methods developed in Refs. Basko ; Aleiner . It is given by the probability $P_{\alpha}$ that for a given one-particle localized state $\ket{\alpha}$ there exist three other states $\ket{\alpha^{\prime}},\ket{\beta},\ket{\beta^{\prime}}$ for which the matrix element of interaction for the transition from the two-particle state $\ket{\alpha,\beta}$ to $\ket{\alpha^{\prime},\beta^{\prime}}$ exceeds the energy mismatch $\Delta_{\alpha\beta}^{\alpha^{\prime}\beta^{\prime}}=|\epsilon_{\alpha}+\epsilon_{\beta}-\epsilon_{\alpha^{\prime}}-\epsilon_{\beta^{\prime}}|$. For a short-range interaction $H_{\text{int}}$, one gets the probability:

[TABLE]

where $C$ is a parameter of order unity.

In Ref. Bertoli , the localized phase is protected by the high-energy truncation of the energy distribution function. On the contrary, if the particle energy grows unbounded, then delocalization takes place above a critical temperature $T_{c}$ that is interaction-independent, i.e. the insulating phase disappears even without interaction between particles. This conclusion, which apparently contradicts the commonly accepted Anderson localization of all single-particle eigenstates in 2D, was made and interpreted in Ref. Nandkishore . The interpretation is based on the exponential increase of the localization length with energy. In order to estimate the “conductivity” one has to integrate the Bose distribution function multiplied by $\exp(-L/\zeta(\epsilon))$, where $L$ is the linear dimension of the system. Evaluating the integral by using the saddle point approximation, one obtains a power-like rather than exponential decrease of the “conductivity” with increasing $L$. The decrease that is slower than exponential can be interpreted as a disappearance of the insulating phase. The peculiarity of this rather academic problem follows from the fact that single-particle energies, which dominate the energy integral, increase logarithmically with $L$ and thus become infinite in the thermodynamic limit. For realistic systems, the exponential growth of the localization length is limited by e.g. a finite size of the system or, as in Ref. Bertoli , by a truncation of the energy distribution function, and the 2D localization is restored.

As noted, in Ref. Nandkishore the author obtained the same result of an interaction-independent critical temperature, albeit with a different method. Namely, the MBLDT criterion in Ref. Nandkishore contained the occupation number of the initial single-particle state. This statement is erroneous as noted in Refs. Gornyi ; Bertoli . However, the authors of Ref. Gornyi concluded that MBL is unstable in continuum systems at any non-zero temperature, after rephrasing the MBLDT criterion as an energy exchange between “hot” and “cold” particles, i.e. particles with high and low energy, with intermediate energies playing no role.

As it is clear from Eq. (1) (see, e.g. Bertoli ), the initial state single-particle occupation number does not enter the criterion for the MBLDT. Nevertheless, it was found Bertoli that for a fixed interaction there exists a range of temperatures $T<T_{c}$ and disorder strengths, where the competition between the exponentials from the localization length (given by Eq. (7) below) and the distribution function $f(\epsilon)=(e^{(\epsilon-\mu)/T)}-1)^{-1}$ is “won” by the latter when increasing energy. An assumption adopted in Ref. Bertoli is that the corresponding initial and final single-particle states are nearest neighbors in energy space, i.e. $\epsilon_{\alpha}\approx\epsilon_{\alpha^{\prime}}$ and $\epsilon_{\beta}\approx\epsilon_{\beta^{\prime}}$. The energies $\epsilon_{\alpha}$ and $\epsilon_{\beta}$ may differ at will. Below, we relax this approximation in the criterion of delocalization. Before describing our results, let us briefly summarize the MBLDT criterion of Ref. Gornyi . The key point of Ref. Gornyi is that one should think of the whole system as containing two subsystems: “hot” and “cold” particles. Initially the cold particles act as a bath for the hot ones, creating delocalized excitations in the hot system, which in turn act as a bath for the cold system. This extends many-body delocalization over the whole spectrum, including intermediate states. From equations (1) and (9) of Ref. Gornyi , we can write the delocalization parameter $\eta_{hc}$, which plays the same role as $C$ in our equation (1), as:

[TABLE]

Here $V_{hc}$ is the matrix element coupling two hot states with two cold states, $N_{\text{eff}}$ is the characteristic number of pairs in the cold system with which a hot particle can hybridize, and $\Delta_{h}$ is just the level spacing for the hot particles. Note that the structure of Eq. (2) is similar to Eq. (1), but involves only coupling between hot and cold particles. The authors argue in favour of an effective $T$-dependent single-particle mobility edge so that states above such energy are always delocalized. This provides transport at any finite temperature, and no insulator is found other than the zero-temperature Bose glass Falco .

Let us now address the question of the stability of the finite-temperature insulator in the context of the model introduced in Ref. Bertoli for 2D interacting disordered bosons. The Hamiltonian of the system reads:

[TABLE]

where

[TABLE]

[TABLE]

Here $\hat{\Psi}(\textbf{r})$ are bosonic field operators. The first term in Eq. (4) is the kinetic energy of particles with mass $m$, while the second term accounts for the random potential $U(\textbf{r})$. The interparticle interaction, Eq. (5), describes a contact interaction between particles, with coupling constant $g>0$. We will consider the case of a disordered Gaussian short-range potential $U(\textbf{r})$ with zero mean, amplitude $U_{0}$, and correlation length $\sigma$. The energy and length scales of the disorder are Lifshitz ; Zittartz :

[TABLE]

In two dimensions, single-particle states are localized with a localization length depending exponentially on the energy. The single-particle localization length in two dimensions at $\epsilon>\epsilon_{\ast}$ can be written in the form Lee :

[TABLE]

so that $\zeta(\epsilon_{\ast})=\zeta_{\ast}$. At energies $|\epsilon|\lesssim\epsilon_{\ast}$, the energy dependence of $\zeta$ is weak, and one can approximate the localization length as $\zeta(\epsilon)\approx\zeta_{\ast}$. The density of states (DoS) for clean 2D bosons in the continuum is energy independent, $\rho_{0}=m/2\pi\hbar^{2}$. The presence of the disordered potential creates the so-called Lifshitz tails, negative-energy states with a DoS decaying exponentially with increasing the absolute value of the energy Lifshitz ; Zittartz . We approximate $\rho(\epsilon)\simeq\rho_{0}$ for energies $\epsilon\geq-\epsilon_{*}$, omitting exponentially small values of the DoS for $\epsilon<-\epsilon_{*}$. Besides that, we consider the weakly interacting regime, with a small parameter given by:

[TABLE]

Here $T_{d}=2\pi\hbar^{2}n/m$ is the degeneracy temperature, with $n$ being the mean density. Finally, we assume weak disorder, so that

[TABLE]

What should happen when a hot localized particle with energy $\epsilon$ and localization length $\zeta(\epsilon)$ is present in the system, and it interacts with a cold cloud? A rough estimate similar to Eq. (2) suggests that the matrix element $V_{hc}$ in this case is proportional to $\sim g/\zeta^{2}(\epsilon)$. The level spacing is just $\Delta_{h}\sim 1/\rho_{0}\zeta^{2}(\epsilon)$. For particles of the cloud which have energies approaching $\epsilon$, the effective number of channels is $N_{\text{eff}}\sim n\zeta^{2}(\epsilon)\exp(-\epsilon/T)$. Putting this altogether we obtain:

[TABLE]

Accordingly, it is possible that $\eta_{hc}\lesssim 1$ provided that $T<\epsilon_{\ast}/2$, which means that there is an insulating phase at these temperatures in the thermodynamic limit, in agreement with Ref. Bertoli . However, the finite-temperature insulator might be merely a consequence of our approximations. Indeed, in Ref. Gornyi the question is addressed by inserting the effective mobility edge, which provides the system with an effective conduction band populated by many delocalized excitations.

We now go back to the MBLDT criterion given by Eq. (1), which takes the form:

[TABLE]

where:

[TABLE]

The prime in the summation in Eq. (11) means that we are summing over a length scale of the localization length (the lowest one among the four states). The factor $N_{\alpha\beta}^{\alpha^{\prime}\beta^{\prime}}$ accounts for the number of possible (direct and inverse) processes $\alpha,\beta\leftrightarrow\alpha^{\prime},\beta^{\prime}$ involving a given state $\alpha$, and $N_{\alpha^{\prime}},N_{\beta},N_{\beta^{\prime}}$ are the occupation numbers (not the averages, and the average is only taken for the square root expression in Eq. (12)). Following Ref. Basko , we replace the difference of square roots with the square root of the difference. The quantity $\psi_{\alpha}(\textbf{r})$ is the one-body wavefunction of a localized state and the bar in the right hand side of (12) means average value. While the true form of the wavefunction is given by an exponentially decaying wavepacket, we make the following approximation:

[TABLE]

Our main goal is to demonstrate the stability of the MBL state. In order to do this we neglect the well known (see, e.g. Mirlin ) dependence of the matrix element of the two-body $\alpha,\beta\rightarrow\alpha^{\prime},\beta^{\prime}$ scattering, $U_{\alpha\beta}^{\alpha^{\prime}\beta^{\prime}}$, on the energy differencies and estimate it as

[TABLE]

Let us now look at equation (12). The quantities $N_{\alpha^{\prime}},N_{\beta},N_{\beta^{\prime}}$ are actually integers representing the occupation of the state in Fock space. If the corresponding average values $\overline{N}_{\alpha^{\prime}},\overline{N}_{\beta},\overline{N}_{\beta^{\prime}}$ are large, then fluctuations are small and we may substitute the average values of the occupation numbers in the r.h.s. of Eq. (12). Assuming that energies $\epsilon$ are almost equal to each other pairwise, i.e. $\epsilon_{\alpha}\approx\epsilon_{\alpha^{\prime}}$ and $\epsilon_{\beta}\approx\epsilon_{\beta^{\prime}}$, one then recovers the expression $N_{\alpha\beta}^{\alpha^{\prime}\beta^{\prime}}\approx\overline{N}_{\beta}$ as in Ref. Bertoli . Going beyond this approximation involves the calculation for each distinct case when some of the average occupation numbers are small. Detailed calculations are given in the Appendix. It turns out that the MBL state is satable if $g<g_{c}(\epsilon_{\alpha})$, where

[TABLE]

For the MBL state to remain stable the condition $g<g_{c}(\epsilon_{\alpha})$ should be satisfied for all $\alpha$, i.e. the condition of the stability is

[TABLE]

According to Eq. (17), the critical coupling $g_{c}(\epsilon_{\alpha})$ is determined by the typical smallest value of the energy mismatch. In order to estimate ${\rm min}\{\Delta_{\alpha\beta}^{\alpha^{\prime}\beta^{\prime}}\}$ relying on Eq. (14) note that all four states $\alpha,\beta,\alpha^{\prime},\beta^{\prime}$ should be localized nearby. The nearest neighbour spacing between such states can be estimated as $\delta=(\rho_{0}\zeta^{2})^{-1}$, where $\rho_{0}$ is the density of states. Since $\zeta$ depends on the energy, one has $\delta=\delta(\epsilon)$. Therefore, we have

[TABLE]

where $\delta_{\beta}\equiv\delta(\epsilon_{\beta})$ and the same for $\alpha^{\prime}$, $\beta^{\prime}$. It should be noted that the assumption of the energy-independence of the matrix element $U_{\alpha\beta}^{\alpha^{\prime}\beta^{\prime}}$, which we adopted, substantially reduces the estimate for ${\rm min}\{\Delta_{\alpha\beta\;\min}^{\alpha^{\prime}\beta^{\prime}}\}$ compared with Ref. Bertoli : the requirement of Bertoli that $\alpha$ and $\alpha^{\prime}$, as well as $\beta$ and $\beta^{\prime}$ should be nearest neighbours in energy leads to ${\rm min}\{\Delta_{\alpha\beta\;\min}^{\alpha^{\prime}\beta^{\prime}}\}={\rm max}(\delta_{\alpha},\delta_{\beta})$. Including the possibility of a smaller denominator in Eq. (17) can favour delocalization.

Below we solve Eq. (17) numerically and first check whether the MBL phase exists in the low-temperature limit. As $N_{\alpha\beta}^{\alpha^{\prime}\beta^{\prime}}$ depends on the chemical potential $\mu$, we establish a relation between $\mu$, $n$, and $T$ from the normalization condition

[TABLE]

In order to calculate the critical coupling $g_{c}$ and its temperature dependence, one has to evaluate the occupation numbers of single-particle states. Following Ref. Michal , we write an expression for the energy corresponding to the configuration $\{N_{\alpha}\}$ of the occupation numbers as:

[TABLE]

The grand canonical partition function becomes:

[TABLE]

where:

[TABLE]

For a large average occupation number $\overline{N}_{\alpha}\gg 1$, fluctuations are small. One linearizes the exponent around $\overline{N}_{\alpha}$ and the partition function reads:

[TABLE]

Dropping the index $\alpha$, we have the following expression for the average occupation numbers of single-particle states on the insulator side:

[TABLE]

For $\overline{N}_{\epsilon}\gg 1$ we expand the exponent in Eq. (25) and obtain:

[TABLE]

For small average occupation numbers $\overline{N}_{\epsilon}\ll 1$ we neglect the interaction term in the exponent of (23). This gives the Boltzmann distribution $\overline{N}_{\epsilon}\approx e^{-(\epsilon-\mu)/T}$.

Below, we set the value of the constant $C=1$. At $T=0$, Eq. (26) gives

[TABLE]

where $\theta(\mu-\epsilon)$ is the Heaviside theta function. Combining equations (17), (19), and (27) we find a critical disorder:

[TABLE]

The corresponding chemical potential is $\mu=1.64ng$. This result is in good agreement with the one from Ref. Bertoli , as well as with the microscopic analysis of tunneling between bosonic lakes Falco . The critical disorder is higher than the one in Ref. Bertoli . This is expected, as we include smaller level spacing than before and more highly energetic processes.

At temperatures $T\ll\epsilon_{\ast}^{2}/T_{d}=m\epsilon_{\ast}^{2}/2\pi\hbar^{2}n$, the critical disorder is practically temperature independent. In the temperature interval, $\epsilon_{\ast}^{2}/T_{d}\ll T\ll\epsilon_{\ast}$, in the thermodynamic limit the average occupation number $\overline{N}_{\epsilon}$ is large when $\epsilon<\mu$. This is because the chemical potential decreases with increasing $T$, and becomes negative when $T$ is a fraction of $\epsilon_{*}$ Bertoli . One integrates Eq. (17) with Eq. (20) using average occupation numbers given by (see also Refs. Michal ; Bertoli ):

[TABLE]

Note that at this stage we do not consider a truncation in the energy distribution.

Fig. 1 shows the obtained results. Remarkably, the physical picture resulting from Ref. Bertoli survives in the presence of a number of processes that were not taken into account there. For very low temperatures $T\ll\epsilon_{*}/2$, the insulator is stable and the critical coupling is only slightly reduced by an increase in temperature. Most importantly, delocalization is driven in this regime by the low-energy states, as we see from Fig. 2, where we plot the value of $\epsilon_{\alpha}$ as a function of $T$. Indeed, we may identify the value of $\epsilon_{\alpha}$ as the energy of the typical states that cause delocalization through the interaction. Finding a low value of $\epsilon_{\alpha}$ for low temperatures implies that the resonant subnetwork in Fock space typically involves states at low energies.

Hot particles start to dominate only when $T\approx 0.47\epsilon_{*}$. Then, delocalization takes place as a result of hybridization of high-energy particles, decaying into three-body excitations. This is compatible with the picture of a hot-cold mixture proposed in Ref. Gornyi and was already noted in Ref. Bertoli when looking at the behavior of $\epsilon_{\alpha}$ in the thermodynamic limit.

Let us remark that in Eq. (17) we take into account all processes that are resonant irrespective of their energies. The only assumption is that two highly energetic states may be taken as energy neighbors when they interact with two cold states, which is consistent with the analysis in the hot-cold mixture. Nevertheless, the obtained results show that localization is present in the low-temperature regime. This is in direct contrast with the proposal of considering the hot and cold subsystems as two separated entities that act as a bath for one another. However, it may well be that this is the situation when $T\to\epsilon_{*}/2$ and beyond.

To complete our analysis, we introduce a truncation in the energy distribution function. The value of the truncation energy $\epsilon_{b}$ for cooling to temperatures $T\lesssim ng$ is typically equal to $ng+\eta T$, with $\eta$ varying between 5 and 8 walraven ; ketterle . We match therefore the zero-temperature result by setting $\epsilon_{b}=1.64ng+\eta T$, and below we use $\eta=5$. Increasing the value of $\eta$ does not change significantly the physical picture. The result is in good agreement with the one in Ref. Bertoli , as shown in Fig. 3. Even though it is slightly reduced with respect to the phase diagram presented in Ref. Bertoli , the insulating phase survives at all temperatures.

It is actually not surprising that the behavior of the system is almost the same, even after relaxing some of the approximations made in Ref. Bertoli . As we already remarked, very hot particles are not present and cannot drive delocalization.

In conclusion, we see that for short-range interacting bosons the MBL phase is stable at temperatures below a critical disorder-dependent value in spite of the presence of processes involving highly energetic particles.

The research leading to these results has received funding from the European Research Council under European Community’s Seventh Framework Programme (FP7/2007-2013 Grant Agreement no.341197).

Appendix

We need to calculate the average value of the square root

[TABLE]

When the average occupation numbers are large, we substitute them directly into equation (30) because fluctuations are small. We get:

[TABLE]

where the average occupation number is given by Eq, (29). Similarly, we keep the largest terms when one of the average occupation numbers is small:

[TABLE]

When two of the average occupation numbers are small, we omit the interparticle interaction for these states and calculate the probability of having $j$ particles in the state with energy $\epsilon$ as:

[TABLE]

Let us look first at the case where $\overline{N}_{\alpha^{\prime}}\gg 1$ and $\overline{N}_{\beta^{\prime}},\overline{N}_{\beta}\ll 1$. We then have:

[TABLE]

The main contribution comes from the terms with

[TABLE]

Taking into account that $e^{-(\epsilon_{\beta(\beta^{\prime})}-\mu)/T}\ll 1$, this yields:

[TABLE]

where $\overline{N}_{\beta^{\prime}}=e^{-(\epsilon_{\beta^{\prime}}-\mu)/T}$ and $\overline{N}_{\beta}=e^{-(\epsilon_{\beta}-\mu)/T}$. For small average occupation numbers $\overline{N}_{\beta}$ and $\overline{N}_{\beta^{\prime}}$, the states $\beta$ and $\beta^{\prime}$ have large energies. We take now the approximation $\epsilon_{\beta}\approx\epsilon_{\beta^{\prime}}$ to find

[TABLE]

Similar calculations give:

[TABLE]

When all three of the average occupation numbers are small, using the same method we have:

[TABLE]

The value of $\epsilon_{\alpha}$ must be chosen in such a way that off-resonant processes are avoided. As in Eq. (39), we set the two highest energy equal to each other (when the average occupation numbers are small) in order to account for the fall-off of the matrix element when one of the energies becomes very high. Using equations (13)-(14), this gives the following MBLDT criterion in the thermodynamic limit:

[TABLE]

We have taken the average occupation number to be large at energies smaller than the chemical potential, and checked that this is a good approximation.