P Lugan

L Sanchez-Palencia

Localization of Bogoliubov quasiparticles in interacting Bose gases with correlated disorder

Keywords: numbers: 05.30.Jp, 03.75.Hh, 64.60.Cn, 79.60.Ht

We study the Anderson localization of Bogoliubov quasiparticles (elementary many-body excitations) in a weakly interacting Bose gas of chemical potential µ subjected to a disordered potential V . We introduce a general mapping (valid for weak inhomogeneous potentials in any dimension) of the Bogoliubov-de Gennes equations onto a single-particle Schrödinger-like equation with an effective potential. For disordered potentials, the Schrödinger-like equation accounts for the scattering and localization properties of the Bogoliubov quasiparticles. We derive analytically the localization lengths for correlated disordered potentials in the one-dimensional geometry. Our approach relies on a perturbative expansion in V /µ, which we develop up to third order, and we discuss the impact of the various perturbation orders. Our predictions are shown to be in very good agreement with direct numerical calculations. We identify different localization regimes: For low energy, the effective disordered potential exhibits a strong screening by the quasicondensate density background, and localization is suppressed. For high-energy excitations, the effective disordered potential reduces to the bare disordered potential, and the localization properties of quasiparticles are the same as for free particles. The maximum of localization is found at intermediate energy when the quasicondensate healing length is of the order of the disorder correlation length. Possible extensions of our work to higher dimensions are also discussed.

I. INTRODUCTION

Disorder in ultracold quantum gases is attracting a growing interest due to unprecedented possibilities of controlling disorder and atom-atom interactions for bosons, fermions or mixtures of atomic species created in one-, two-or threedimensional traps [1][2][3][4]. So far much attention has been devoted to studies of the disorder-induced damping of motion in Bose [5][6][7][8][9][10] and Fermi [11] gases, classical localization [12][13][14][15][16][17][18][19][20], spatial diffusion [19][20][21][22][23][24][25] and Anderson localization in regimes where interactions can be neglected [18,21,24,[26][27][28][29][30][31][32][33][34][35][36][37][38][39]. The effects of disorder in interacting quantum systems have also been studied in a variety of contexts, such as transport in weakly interacting Bose-Einstein condensates [40][41][42][43][44][45][46][47][48], interacting Bose gases at equilibrium [26,, strongly interacting Fermi gases [73][74][75][76] and coupled two-component gases [77][78][79][80][81][82][83].

The interplay of disorder and interactions in quantum systems is an issue of fundamental importance to understand the behaviors of superfluid 4 He in porous media [84][85][86][87][88], dirty superconductors [89][90][91][92][93] and ultracold gases in optical disorder [55,[66][67][68]70]. Although a number of questions are open, in particular regarding the fate of Anderson localization, general behaviors can be found in various situations. For instance, weak repulsive interactions in a Bose gas at zero temperature in a disordered potential generally lead to delocalization [53,60,63,94]: In the absence of interactions, all bosons condense into the single-particle ground state, which is localized [95]. This kind of N -body Fock state is highly unstable in an infinite system where an infinity of spatially-separated single-particle states coexist at arbitrarily close energies [START_REF] Lifshits | Introduction to the Theory of Disordered Systems[END_REF]. Then, for weak repulsive interactions, the Bose gas fragments into a number of low-energy, localized single-particle states, so as to minimize the interaction en-ergy [53,60,63]. The Bose gas forms a gapless compressible insulating phase, known as the Bose glass [94,[START_REF] Ma | [END_REF][98][99][100][101][102][103]. For increasing mean-field interactions, the fragments merge and form a single extended condensate, which restores superfluidity [53,60,63,[104][105][106][107][108]. Finally, in the strongly interacting regime, repulsive interactions can finally destroy again superfluidity, forming Tonks-Girardeau gases in 1D [99,[109][110][111] or disordered Mott insulators in lattice gases [100][101][102].

The above results lead to the conclusion that, at zero temperature where only the ground state plays a role, repulsive interactions destroy Anderson localization in Bose gases for moderate interaction strengths that are compatible with the mean-field approach [53,94,104]. At non-zero temperature however, important properties, such as the correlation functions, phase coherence and long-range order [105,112], are determined by the excitations of the Bose gas, which are populated thermally. It is thus of prime importance to determine how disorder affects the behavior of the latter. These excitations of the many-body system are of collective nature. They can be viewed as quasiparticles, scattering on the disordered potential. In contrast to the mean-field background, which is extended, the Bogoliubov quasiparticles of a weakly interacting Bose gas can be localized in the presence of a disordered potential [51,54,58,71,72].

In this paper, following the approach of Ref. [54], we present a detailed theory of the Anderson localization of Bogoliubov quasiparticles in weakly interacting Bose gases subjected to correlated disordered potentials. On the one hand, this approach applies to any kind of weak disordered potentials with short or long-range correlations and offers a unified theory valid all along the crossover from the phonon regime to the free particle regime. On the other hand, it permits a straightforward interpretation of the effect of repulsive interactions in terms of a screening of the disorder by the density background. In Sec. II, the grand-canonical Hamiltonian of the interacting many-body system is expanded up to second order in phase and density fluctuations. The reduced Hamiltonian is quadratic and can be diagonalized by following the standard Bogoliubov-Popov approach. The excitations are the solutions of the Bogoliubov-de Gennes equations, which consist in a set of two coupled differential equations of order two. In Sec. III, we then introduce a general mapping, valid for weak inhomogeneous potentials in any dimension, of the Bogoliubov-de Gennes equations onto a single-particle Schrödinger-like equation with a screened potential. For disordered potentials, the effective Schrödinger-like equation describes the scattering and localization properties of the Bogoliubov quasiparticles. In Sec. IV, we apply this approach to study the Anderson localization of Bogoliubov quasiparticles in one-dimensional disorder. We derive analytical formulas for the localization length of the Bogoliubov quasiparticles up to third order in perturbation theory and compare our predictions to direct numerical calculations. Our results exhibit different localization regimes: For low energy, the effective disordered potential accounts for a strong screening by the quasicondensate density background and Anderson localization is suppressed. For high energy excitations, the screening is small; the effective disordered potential reduces to the bare disordered potential and the localization properties of quasiparticles are the same as for free particles. The maximum of localization is found at intermediate energy when the quasicondensate healing length is of the order of the disorder correlation length. Finally, in Sec. V, we summarize our results and discuss possible extensions of our work, in particular towards higher dimensions.

II. ELEMENTARY (BOGOLIUBOV) EXCITATIONS IN A BOSE GAS WITH WEAK DENSITY FLUCTUATIONS

We consider a d-dimensional, ultracold, dilute gas of bosons with weak repulsive interactions, in a potential V (r). The system is described by the grand-canonical Hamiltonian

K = Ĥ -µ N = dr 2 2m (∇ θ) 2 n + (∇ √ n) 2 +V (r)n + g 2 n2 -µn , (1)
where the short-range atom-atom interactions are modeled by a contact potential with coupling constant g > 0, m is the atomic mass, n and θ are the density and phase operators, which satisfy the commutation relation [n(r), θ(r ′ )] = iδ(rr ′ ) [START_REF] Barnett | Methods in Theoretical Quantum Optics[END_REF], and µ is the chemical potential. In the full form (1), solving K for eigenstates is difficult in general. Yet, for small density fluctuations around n c = n (i.e., for |δn| ≪ n c , where |δn| is the typical value of δn = n -n c in the state of the system) and for small phase gradient (2 |∇ θ| 2 /2m ≪ µ), the operator K can be expanded around the classical field n c , ∇ θ = 0, following the Bogoliubov-Popov approach [START_REF] Bogolyubov | [END_REF][115][116][START_REF] Popov | Functional Integrals in Quantum Field Theory and Statistical Physics[END_REF][START_REF] Mora | [END_REF][119]. 1 In the zeroth-order expansion the ground-state density profile is found by minimizing the grand-canonical energy functional associated with the Hamiltonian (1) with respect to the variation of n c (r). This yields the Gross-Pitaevskii equation (GPE)

- 2 2m ∇ 2 + V (r) + gn c (r) -µ n c (r) = 0. (2)
Then, retaining only the leading terms in the expansion of the density fluctuations δn and phase gradients ∇ θ, the Hamiltonian ( 1) is cast into the form

K = E 0 + ν ǫ ν b † ν bν
, where b † ν and bν are the bosonic creation and annihilation operators of an excitation (Bogoliubov quasiparticle, BQP) of energy ǫ ν . 2 The phase and density operators are expanded as

θ(r) = -i 2 n c (r) ν f + ν (r) bν -H.c. ( 3
)
δn(r) = n c (r) ν f - ν (r) bν + H.c. , (4)
where the Bogoliubov wavefunctions f ± ν (r) are solutions of the Bogoliubov-de Gennes equations (BdGEs) [START_REF] De Gennes | Superconductivity of Metals and Alloys[END_REF] -

2 2m ∇ 2 + V (r) + gn c (r) -µ f + ν (r) = ǫ ν f - ν (r) (5)
- 2 2m ∇ 2 + V (r) + 3gn c (r) -µ f - ν (r) = ǫ ν f + ν (r) (6)
with the normalization

dr f + ν (r)f - ν ′ * (r) + f - ν (r)f + ν ′ * (r) = 2δ ν,ν ′ . ( 7
)
Equations ( 3) and (4) reveal the simple physical meaning of the functions f + ν (r) and f - ν (r). Up to the factor n c (r), they describe the spatial dependence of the phase and density fluctuations associated with the BQPs, respectively. Notice that, as first pointed out by Popov [116,[START_REF] Popov | Functional Integrals in Quantum Field Theory and Statistical Physics[END_REF], the above derivation of the BdGEs in the phase-density representation provides an extension of the usual Bogoliubov-de Gennes theory. In the latter, the starting point consists in applying the usual Bogoliubov shift to the field operator, Ψ = √ n c +δ Ψ, and expanding K up to quadratic terms in the fluctuation δ Ψ [START_REF] Bogolyubov | [END_REF]115,[START_REF] De Gennes | Superconductivity of Metals and Alloys[END_REF]. This approach assumes weak phase and density fluctuations around a unique classical field √ n c which breaks the U (1) phase symmetry of the Hamiltonian. On the contrary, the phase-density picture used in this work does not rely on this assumption. In particular, it provides a satisfactory description of the mean-field ground state n c (r) and the excitations of Bose gases in the quasicondensate regime [START_REF] Petrov | [END_REF]124].

Within the above formalism, Eqs. ( 2), ( 5) and (6) form a closed set which describes non-interacting, bosonic quasiparticles. Interactions between these quasiparticles only arise with higher-order terms in the expansion of K, which we neglect here. Therefore, to study the low-temperature properties of the Bose gas in the external potential V (r), we are left with the sole modes defined by the GPE (2) and BdGEs ( 5) and ( 6). Still, the set of equations ( 2), ( 5) and ( 6) remains difficult to solve in general, as the GPE (2) is nonlinear, and the two coupled second-order BdGEs ( 5) and ( 6) themselves amount to a differential problem of order four. In the following, we develop a perturbative approach, valid in the limit of a weak potential V (r), which enables us to solve Eqs. ( 2), ( 5) and ( 6) rigorously, and to interpret the underlying physics in simple terms.

III. A SCHR ÖDINGER-LIKE EQUATION FOR BOGOLIUBOV EXCITATIONS IN WEAK POTENTIALS

From now on, we assume that V (r) is a weak external potential with a vanishing average (V = 0), 3 and a typical amplitude V R , 4 such that V 2 R = V 2 . While a less stringent weakness criterion can be derived (see below and Ref. [104]), |V R | ≪ µ is a sufficient assumption to start with. Note that, although we will focus on the case of a disordered potential in the following, the perturbative approach introduced here is general, and V need not be a disordered potential. In any case, we write the autocorrelation function of V , C 2 (r ′r) = V (r)V (r ′ ), in the dimensionless form

C 2 (r) = V 2 R c 2 (r/σ R ), (8)
where σ R is a characteristic length scale of V , which will be precisely defined below when needed. In the following paragraphs, we solve the GPE (2) for the ground-state density n c (Sec. III A), and we use the result to reduce the BdGEs ( 5) and ( 6) to a single Schrödinger-like equation (Secs. III B and III C).

A. The (quasi-)BEC density background

In the regime where the repulsive interactions are strong enough (while remaining compatible with the mean-field 3 When V is a disordered potential, we assume spatial homogeneity, so that the spatial average of V coincides with its statistical average [START_REF] Lifshits | Introduction to the Theory of Disordered Systems[END_REF]. For non-disordered potentials, V denotes the spatial average of V . 4 The sign of VR becomes relevant in the description of disordered potentials with asymmetric single-point probability distribution (see e.g. Ref [1]).

regime), i.e., when the healing length

ξ = √ 4mµ (9)
is much smaller than the size of the system L, the density profile n c is homogeneous in the absence of an external potential, and remains extended (delocalized) for a weak potential V , owing to the nonlinear term in the GPE (2) [53,94,104]. It then proves useful to write the mean-field density term in the form

n c (r) = µ + ∆ -V (r) g , (10)
where

V (r) = gn c -gn c (r) (11)
contains the inhomogeneous part such that V = 0, and

∆ = gn c -µ (12)
represents the mean deviation from the mean-field equation of state µ = gn c that holds in the homogeneous case (V = 0). The quantities V and ∆ both vanish for V = 0, and are expected to remain small for a weak external potential V and repulsive interactions that are strong enough. We then write the perturbation expansions of these quantities in increasing powers of V R /µ:

V (r) = V (1) (r) + V (2) (r) + . . . (13)
∆ = ∆ (1) + ∆ (2) + . . . .

The various terms can be calculated by generalizing the approach of Ref. [104] beyond the first order. Details of these calculations are presented in Appendix A (see also Ref. [72]). Below, we only discuss the main results.

First correction to the mean-field equation of state

The first-order term ∆ (1) in the deviation ∆ vanishes [104]. The leading term is thus provided by ∆ (2) , which depends explicitly on the potential V and on the healing length ξ through [see Eq. (A16) in Appendix A]:

∆ (2) = V 2 R σ d R 2(2π) d/2 µ dq (|q|ξ) 2 [1 + (|q|ξ) 2 ] 2 ĉ2 (qσ R ), ( 15
)
where ĉ2 is the Fourier transform5 of the reduced autocorrelation function c 2 defined in Eq. (8). As ĉ2 is a positive function by virtue of the Wiener-Khinchin theorem, we always have

0.3 0.2 0.1 0 ( 2µ / V R 2 ) ∆ 4 3 2 1 0 σ R / ξ speckle ∆ speckle ∆ (2) monochromatic ∆ monochromatic ∆ (2)
Figure 1. (color online) Comparison of the leading-order correction ∆ (2) to the mean-field equation of state with numerical computations of ∆ for a 1D speckle potential with a reduced autocorrelation function c2(u) = sin(u) 2 /u 2 as used in experiments [12,15,29] and for a 1D periodic lattice V (z) = VR √ 2 cos(z/σR). Here VR = 0.1µ.

∆ (2) > 0, i.e., µ < gn c in the presence of an external potential (see also Ref. [94]). In Fig. 1, Eq. ( 15) is compared to exact numerical calculations of ∆ for a disordered potential and a monochromatic lattice potential, with various values of the ratio σ R /ξ, in a 1D geometry. As expected, the agreement is good for values of V R /µ as low as in Fig. 1. We checked that the small discrepancy between ∆ (2) and ∆ in the disordered case is due to contributions of the order of V 3 R , which are absent in a monochromatic lattice. This validates the perturbative approach.

Let us examine some limiting cases. In the Thomas-Fermi limit, ξ ≪ σ R , we find

∆ (2) ≃ V 2 R 2(2π) d/2 µ ξ σ R 2 dκ |κ| 2 ĉ2 (κ), (16)
so that ∆ (2) /µ ∝ (V R /µ) 2 (ξ/σ R ) 2 ≪ (V R /µ) 2 .
The opposite limit, ξ ≫ σ R , corresponds in principle to the white-noise limit which is obtained by letting the ratio σ R /ξ vanish while keeping the product

V 2 R σ d R constant.
Then, the Fourier transform of the reduced autocorrelation function may be approximated by a constant, ĉ2 (κ) ≃ ĉ2 (0). In fact, we find that the white-noise limit of expression ( 15) is correctly defined only in 1D, for which we obtain

∆ (2) ≃ √ πV 2 R σ R 4 √ 2µξ ĉ2 (0), (17)
so that ∆ (2) /µ ∝ (V R /µ) 2 (σ R /ξ) ≪ (V R /µ) 2 .
In dimension higher than one, this limit cannot be defined because the integrand dq (|q|ξ) 2 / 1 + (|q|ξ) 2 2 scales as q d-3 dq for high momenta (q ≫ ξ -1 ). The integral in Eq. ( 15) would thus be plagued by an ultraviolet divergence for d ≥ 2 and a constant ĉ2 (κ) = ĉ2 (0). In other words, for d ≥ 2, the quantity ∆ (2) depends crucially on the precise form of the reduced autocorrelation function c 2 .

Inhomogeneous part of the (quasi-)BEC density

In contrast to the mean deviation ∆, the leading contribution to the inhomogeneous part of the density profile is provided by the first-order term, which reads

V (1) (r) = dr ′ G ξ (r -r ′ )V (r ′ ), (18)
where G ξ is the Green function associated with the differential operator -ξ 2 ∇ 2 + 1 (see Ref. [104] and Eq. (A15) in Appendix A). In Fourier space we have

G ξ (q) = (2π) -d/2 1 + (|q|ξ) 2 , (19)
so that

V (1) (q) = V (q) 1 + (|q|ξ) 2 . ( 20
)
The healing length ξ clearly appears as a threshold length scale in the response of the density n c (r) to the external potential V (r). Indeed, we have

V (1) (q) ≃ V (q) for |q| ≪ ξ -1 , whereas V (1) (q) ≪ V (q) for |q| ≫ ξ -1 .
In other words, the potential V (1) (r) follows the spatial modulations of V (r) while evening out the high-frequency components. It also follows from the Parseval-Plancherel theorem that | V

R | ≤ |V R |, where | V (1) (1)

R

| is the standard deviation of V (1) and the sign of V

R is chosen to be the same as that of V R . The potential V (1) (r) is thus termed a smoothed potential [104].

If V (r) is a homogeneous disordered potential, that is, a disordered potential whose statistical properties do not depend on the position r [START_REF] Lifshits | Introduction to the Theory of Disordered Systems[END_REF], then so is V (1) (r). If V (r) is a periodic potential, V (1) (r) is also periodic with the same period, but a smoothed Bloch amplitude in each periodic cell, and simply rescaled Fourier components, as shown in Fig. 2. In either case, Eq. ( 10) implies that, if V (r) is a homogeneous potential, the density profile n c is extended [53,104]. Note also that the first-order term V (1) is indeed a small perturbation of the homogeneous density profile whenever | V

(1) R | ≪ µ, which loosens the initial weakness criterion |V R | ≪ µ.
The smoothing solution (10) at first order (∆ = 0 and V = V (1) ) is compared to exact numerical computations of the density profile in Fig. 2, both for a disordered and a bichromatic periodic potential, in coordinate and Fourier spaces. For the sake of clarity, the power spectra for the disordered potential [Fig. 2(c)] have been smoothed by a running average of width 0.1 in qξ units. The agreement between the first-order smoothing solution and the numerical results is good, especially when compared to the Thomas-Fermi limit [for which V (1) (z) = V (z)]. For the bichromatic periodic potential, the appearance of red dots around qξ = ±1 in Fig. 2(d) results from the admixture of the two components of the bichromatic lattice potential due to the nonlinearity of the GPE, as described by the first nonlinear correction V (2) (see appendix A). The other higher-order components, which are several orders of magnitude smaller, are not shown.

1.0 0.5 0.0 |(gn c -gn c )(q)| 2 [arb.units] -2 -1 0 1 2 q ξ 1.0 0.5 0.0 |(gn c -gn c )(q)| 2 [arb.units] -2 -1 0 1 2 q ξ 1.

B. Bogoliubov-de Gennes equations in the decoupling basis

With the BEC background solution n c (r) at hand, we can now solve the BdGEs ( 5) and (6). Using Eq. ( 10), these equations can be cast into

- 2 2m ∇ 2 + V + ∆ -V f + ν = ǫ ν f - ν (21)
- 2 2m ∇ 2 + 2µ + V + 3∆ -3 V f - ν = ǫ ν f + ν , ( 22
)
where ∆, V (r) and V (r) are small compared to µ. As n c (r) is the ground-state solution of the GPE (2), the (quasi-) condensate is dynamically stable (see e.g. Refs. [START_REF] Castin | Proceedings of the Les Houches Summer School of Theoretical Physics[END_REF]126]), and we need only consider real-valued, strictly positive eigenvalues of the BdGEs. Now, given such an eigenvalue ǫ ν , we are interested in the properties of the corresponding mode {f + ν (r), f - ν (r)}. Following the approach developed in Ref. [54], we take advantage of the structure of the eigenmodes before perturbation by a weak potential and introduce an alternative representation of the BQPs in terms of the components {g + ν (r), g - ν (r)}, where the functions g ± and f ± are related by the linear transformation

g ± ν (r) = ±ρ ±1/2 ν f + ν (r) + ρ ∓1/2 ν f - ν (r), (23)
with

ρ ν = µ ǫ ν + 1 + µ ǫ ν 2 . (24)
Expression (24) defines ρ ν as a function of the eigenvalue ǫ ν which does not depend on the details of the specific mode under consideration. The linear transformation (23), derived in Appendix B, is chosen is such a way that the coupling between g + ν (r) and g - ν (r) vanishes in the homogeneous case V = 0. As such, it offers a convenient starting point for a perturbation expansion in the case of weak potentials, as shown below.

In the basis of the g ± ν (r) functions, the BdGEs take the exact form [see Eqs. (B13)-(B14) in Appendix B]:

2 k 2 ν 2m g + ν = - 2 2m ∇ 2 g + ν - 2ρ ν 1 + ρ 2 ν ( V -∆)g - ν + V - 3 + ρ 2 ν 1 + ρ 2 ν ( V -∆) g + ν ( 25
)
- 2 β 2 ν 2m g - ν = - 2 2m ∇ 2 g - ν - 2ρ ν 1 + ρ 2 ν ( V -∆)g + ν + V - 1 + 3ρ 2 ν 1 + ρ 2 ν ( V -∆) g - ν , (26)
where

2 k 2 ν 2m = µ 2 + ǫ 2 ν -µ (27)
2 β 2 ν 2m = µ 2 + ǫ 2 ν + µ. ( 28
)
Both k and β are real-valued functions of the energy ǫ. As a consequence, the associated g + ν and g - ν functions are essentially of the oscillating and the evanescent type, respectively, owing to the signs of the l.h.s. terms in Eqs. (25) and (26). This is consistent with the limit of a vanishing external potential (V = 0, and thus V = 0, ∆ = 0), where the equations for g + ν and g - ν are decoupled. It the latter case, the quantity k can be identified with the wave number |k| of an oscillating, plane-wave BQP mode of energy ǫ k , and Eq. ( 27) is equivalent to the usual Bogoliubov dispersion relation

ǫ k = 2 k 2 2m 2 k 2 2m + 2µ . (29)
The β coefficients, on the other hand, characterize a subset of solutions to the BdGEs that are forbidden when V = 0, as the boundary conditions imposed on the system (e.g. periodic or homogeneous Dirichlet boundary conditions) preclude monotonously growing or decreasing BQP components.

C. Effective Schrödinger equation

While the g - ν function vanishes identically in the absence of an external potential, this is no longer true when V couples g - ν to g + ν via Eqs. ( 25) and (26). For a weak external potential, however, all the terms introduced by g - ν in Eq. ( 25) are at least of second order in V R [see Eq. (B16) and the discussion below], so that we can neglect the second term on the r.h.s. of Eq. (25). Besides, the terms proportional to ∆ in Eqs. ( 25) and ( 26) are at least of second order in the disorder amplitude V R (see Sec. III A 1), and can also be neglected in a first-order approach. We are thus left with the following closed equation for g + ν , which is valid to first order in V R :

- 2 2m ∇ 2 g + ν (r) + V kν (r)g + (r) ≃ 2 k 2 ν 2m g + ν (r), (30)
where

V kν (r) = V (r) - 3 + ρ 2 ν 1 + ρ 2 ν V (1) (r). (31)
Equation ( 30) is formally equivalent to a Schrödinger equation for a bare particle of energy 2 k 2 ν /2m in an effective potential V kν (r). 6 The potential V kν (r) differs from both the bare potential V (r) and the smoothed potential V (r) and explicitly depends on the BQP energy ǫ ν via the parameter ρ ν . These features are illustrated in Fig. 3 where we plot a given realization of a 1D speckle potential V (z), together with the corresponding smoothed potential V (z) and effective potential V kν (hereafter named screened potential on grounds explained below), for two values of the BQP energy.

For further convenience, the dependence of the effective potential on the BQP energy ǫ ν is expressed by the subscript k ν , with the understanding that k ν is defined by Eq. ( 27). Combining Eqs. ( 24) and ( 27), we find

ρ ν = 1 + 1 (k ν ξ) 2 , ( 32
)
and hence

V kν (r) = V (r) - 1 + 4(k ν ξ) 2 1 + 2(k ν ξ) 2 V (1) (r). (33)
To gain more insight into the properties of V kν (r), let us turn to Fourier space where, by virtue of Eq. ( 20), the effective potential reads

V kν (q) = V (q) 1 - 1 + 4(k ν ξ) 2 1 + 2(k ν ξ) 2 1 1 + (|q|ξ) 2 . ( 34
)
Upon inspection, this expression shows that we have |V kν (q)| ≤ |V (q)| for any Fourier component q and any BQP energy ǫ ν . Note also that, by construction of V (1) , the potential V kν has a vanishing average. Hence, keeping in mind that V kν results from the competition of the bare potential V and the BEC background n c in the BdGEs, we term V kν a screened potential. The screening thus affects all Fourier components of the external potential in any dimension. Once Eq. ( 30) has been solved (possibly self-consistently) for g + ν , the function g - ν can be computed from g + ν . Indeed, since 2 β 2 ν /2m > µ [see Eq. ( 28)] and for |V R | ≪ µ, the last term in the r.h.s. of Eq. ( 26) can be neglected, and we find

g - ν (r) ≃ 2m 2 β 2 ν 2ρ ν 1 + ρ 2 ν dr ′ G 1/β (r -r ′ ) V (1) (r ′ )g + ν (r ′ ), (35) where
G 1/β (q) = (2π) -d/2
1+(|q/β|) 2 is the Green function associated with the differential operator -(1/β) 2 ∇ 2 + 1, written in Fourier space.

For analytical purposes, the g ± ν functions hence usefully replace the physically meaningful quantities f ± ν , which can readily be recovered by inverting transformation (23). In particular, as far as asymptotic localization properties in disordered potentials are concerned, Eq. ( 35) tells us that the typical amplitude of g - ν evolves parallel to the amplitude of g + ν on intermediate to long length scales, if V is homogeneous. In this respect, the benefit of the mapping of the exact BdGEs onto Eqs. ( 30) and ( 31) is that we can apply standard techniques for bare Schrödinger particles in weak disordered potentials, in any dimension, as long as these are consistent with the lowest-order approximation used to derive the effective equation ( 30). Yet, BQPs differ substantially from usual bare particles in their scattering and localization properties, because of the peculiar features of the screened potential V kν (r).

IV. LOCALIZATION OF BOGOLIUBOV QUASIPARTICLES IN ONE DIMENSION

The formalism developed in Sec. III is valid for any weak potential. It is particularly fruitful when applied to disordered potentials as known theories developed for single (non-interacting) particles can be directly applied to the effective Schrödinger-like equation (30). In this section, we focus on the one-dimensional geometry using the so-called phase formalism [START_REF] Lifshits | Introduction to the Theory of Disordered Systems[END_REF], which allows for an exact perturbative calculation of the Lyapunov exponent (inverse localization length). This procedure can be straightforwardly extended to higher dimensions, applying appropriate single-particle theories to Eq. ( 30), for instance the self-consistent theory of localization [127,128].

A. Phase formalism in the Born approximation

In 1D, the Lyapunov exponent γ k of a bare particle of energy E k in a disordered potential V (z) is simply related to the backscattering amplitude of the particle from the inhomogeneities of V . For a weak disorder, the Lyapunov exponent can be extracted from a perturbation expansion γ k = γ

(2) k + γ (3) k + • • • in powers of the disorder amplitude V R .
For γ k ≪ k, this approach yields the following result in the lowest-order (Born) approximation [START_REF] Lifshits | Introduction to the Theory of Disordered Systems[END_REF]:

γ (2) k = √ 2π 8k 2 2m 2 2 C 2 (2k), (36)
where C 2 (q) is the Fourier transform of C 2 (z), evaluated at q = 2k. In this formulation, V (z) has a vanishing average, and k, explicitly defined as k = √ 2mE k / , stands for the typical wave vector of the particle under consideration. As such, the parameter k, rather than the related energy E k , is the meaningful quantity in the interference effect that causes Anderson localization. Applying this result to the Schrödingerlike equation (30), we derive the Lyapunov exponent Γ k of a BQP in the Born approximation

Γ k ≃ Γ (2) k : Γ (2) k = √ 2π 8k 2 2m 2 2 C 2,k (2k), (37)
where C 2,k (q) is the Fourier transform of the two-point correlator C 2,k (z) = V k (z ′ )V k (z ′ + z), and k depends on the energy ǫ through Eq. ( 27) as we are now dealing with BQPs.

From the Wiener-Khinchin theorem, we have C 2,k (q) ∝ |V k (q)| 2 , so that, according to Eq. ( 31), the relevant Fourier component of V k for the calculation of Γ

(2)

k is V k (2k) = V (2k) - 1 + 4(kξ) 2 1 + 2(kξ) 2 V (1) (2k). (38)
Then, inserting Eq. ( 20) into Eq. ( 38), we obtain

V k (2k) = S(kξ)V (2k), (39)
where

S(kξ) = 2(kξ) 2 1 + 2(kξ) 2 . (40)
Finally, Eq. ( 37) can be rewritten as where

Γ (2) k = [S(kξ)] 2 γ (2) k , (41)
γ (2) k = √ 2π 32 V R µ 2 σ R k 2 ξ 4 ĉ2 (2kσ R ). (42)
Equation ( 41), together with Eqs. ( 27), (40), and ( 42), completely determines the Lyapunov exponent of a BQP of energy ǫ in a weak, correlated, 1D disordered potential. Remarkably, Eq. ( 41) shows that the Lyapunov exponent of a BQP can be simply related to that of a bare Schrödinger particle with the same average wave vector k [54]. The effects of interactions, in particular, are entirely absorbed in the term S(kξ) defined in Eq. ( 40), which we call screening function on the basis of the analysis presented in section IV B below. This function is shown in Fig. 4 and displays two regimes which can be traced back to the nature of the elementary excitations of the interacting Bose gas in the absence of disorder. In the homogeneous case (V = 0), the elementary (BQP) excitations of the Bose gas undergo a crossover from a regime of pair excitations with linear dispersion relation for ǫ ≪ µ, ǫ k ≃ ck with c = µ/m the speed of sound (phonon regime; PH), to a regime of single-particle excitations with a quadratic dispersion relation ǫ k ≃ E k = ( 2 /2m)k 2 for ǫ ≫ µ (free-particle regime; FP). Hence, while the localization of single particles merely results from the competition of their kinetic energy with the disorder amplitude and correlation [33], the localization of BQPs also crucially depends on interaction-induced particle correlations.

B. Localization regimes

Summarizing the results of the previous section, we find that the Lyapunov exponent of BQPs in weak 1D disordered potentials is given by the expression

Γ (2) k = √ 2π 8 V R µ 2 k 2 σ R [1 + 2(kξ) 2 ] 2 ĉ2 (2kσ R ), (43)
obtained by combining Eqs. ( 40), ( 41) and (42). In this expression, the quadratic dependence on the potential amplitude V R is characteristic of the Born approximation. The scaled Lyapunov exponent (µ/V R ) 2 Γ

(2)

k ξ depends only on the two parameters kξ and σ R /ξ. Expression (43) nevertheless contains contributions of three distinct physical origins, which appear more clearly in Eqs. ( 41) and ( 42): i) a 1/k 2 term which is representative of the kinetic energy of a bare particle, ii) the squared screening function S(kξ) 2 which accounts for particle interactions, and iii) the spectral density of disorder V 2 R σ R ĉ2 (2kσ R ) at the wave vector 2k. The role of these various contributions is discussed below.

Screening in the phonon regime

The interplay of the first two contributions is best understood by studying the case of a white-noise potential, which is obtained in the limit σ R → 0, |V R | → ∞, with V 2 R σ R = const. In this limit, the spectral density ĉ2 uniformly approaches a flat distribution with an amplitude of the order of one. Then, in the free-particle regime ǫ ≫ µ (i.e., kξ ≫ 1), we have S(kξ) ≃ 1 and Γ

(2) k ∼ 1/k 2 ∼ 1/ǫ. In other words, BQPs localize exactly like bare Schrödinger particles in this regime, as expected. In the phonon regime, on the contrary, the kinetic term is dominated by the S(kξ) 2 factor, which is approximately quartic in k (or ǫ). We then get the scaling 54], which is consistent with known results on the localization of acoustic phonons in 1D [129][130][131][132].

Γ (2) k ∼ k 2 ∼ ǫ 2 [51,
Interestingly, Eq. ( 41) combines the two limiting models in a unified picture, and provides a physical interpretation for the decreasing localization of phonon modes with decreasing energy. The function S(kξ) reflects the competition of the bare external potential V and the interaction of the BQPs with the quasi-BEC density background gn c , which appears here as V . In particular, the strong decay of S(kξ) in the phonon regime can be interpreted as an increasing screening of the external potential by the static quasi-BEC background, which adapts to the long-wavelength modulations of the disordered potential (see Sec. III A 2).

Correlated disordered potentials

To analyze the role of the correlation length σ R in Eq. ( 43), we consider optical speckle potentials, which are now widely used with ultracold atoms for their tunability and truly random properties, as a model of correlated disorder [5,6,9,12,13,15,19,29,66,67,133]. In the simplest case where the speckle pattern is obtained at the back focal plane of a lens with rectangular aperture and uniform illumination (see e.g. Ref. [15]), the reduced autocorrelation function reads

c 2 (u) = sin(u) 2 /u 2 , ( 44
)
where u = z/σ R . The corresponding Fourier spectrum is where Θ is the Heaviside step function. Then, Eq. ( 43) reads

ĉ2 (kσ R ) = π 2 1 - kσ R 2 Θ 1 - kσ R 2 , (45)
Γ (2) k = π 8 V R µ 2 k 2 σ R (1 -kσ R ) [1 + 2(kξ) 2 ] 2 Θ(1 -kσ R ), (46)
which is shown in Fig. 5. Such disorder correlations introduce several features, which we discuss below.

Effective mobility edge -Equation (46) shows that the Lyapunov exponent in the Born approximation, Γ

k , vanishes identically for kσ R > 1 (see also Figs. 5 and6). This feature originates from the special correlation properties of speckle potentials, the power spectrum of which has a highmomentum cutoff [see Eq. ( 45)], i.e., contains no 2k component able to backscatter a wave travelling with wavevector k > 1/σ R according to the elastic process +k → -k at the level of the Born approximation [18,27]. The Born approximation consists in truncating the perturbation series in powers of V R used to derive Lyapunov exponents after the leading order, proportional to V 2 R . In agreement with the understanding of localization in 1D, whereby single particles are localized (i.e., γ k > 0) at all energies under fairly general assumptions [95], higher-order terms in the perturbation expansions of γ k and Γ k are expected to carry contributions which do not all vanish identically for kσ R > 1. However, the larger power dependence on the small parameter V R /µ makes these terms negligible in the limit of weak disorder. In this limit, the Lyapunov exponents thus experience a sharp crossover and drop by orders of magnitude when k is varied through the value 1/σ R . Such a crossover characterizes an effective mobility edge which strongly affects localization properties in finite-size systems. As a matter of fact, the thirdorder contribution to Γ k , proportional to V 3 R , can be shown to vanish abruptly for momenta above the same cutoff at 1/σ R (see Sec. IV E), so that corrections to Eq. ( 46) beyond that cutoff scale as V 4 R at least 7 . This behavior is specific to potentials with cutoffs in their Fourier-transformed correlation functions [27,32,33,134,135].

Localization maxima -In the white-noise limit, the BQPs localize best for kξ = 1/ √ 2 (i.e. ǫ = √ 3µ), that is, in the cross-over region between the phonon and the free-particle regime [51]. This behavior results from the competition of bare kinetic energy and mean-field interactions via the screening effect, as discussed above. With correlations, however, the detailed statistical properties of the disorder play a role as well, and the wave vector of maximum localization k max decreases with increasing correlation length σ R . More generally, it can be checked from Eq. ( 43) that for any correlated disorder with monotonously decreasing power spectrum akin to model (45), the wave vector of maximum localization is shifted to lower values than the corresponding white-noise value. The locus of k max as a function of the correlation and healing lengths is plotted in green solid line in Fig. 5. For each σ R /ξ ratio, we indeed find a unique maximum k max with

k max ≃ 1 √ 2 ξ 1 - σ R /ξ 2 √ 2 , σ R ≪ ξ, ( 47
)
k max ≃ 2 3 σ R , σ R ≫ ξ. ( 48
)
These asymptotic expressions show that k max is controlled by the longest length scale in the problem. Finally, we find an absolute maximum at fixed ξ for σ R = 3/2 ξ and kξ = 1/ √ 6, which yields a localization length

L max (ξ) = Γ -1 max (ξ) = 512 √ 6 9π µ V R 2 ξ. (49)
Current experiments with ultracold atoms implement disordered potentials with correlation lengths of the order of σ R ≃ 0.25µm [19,29,66,67], which yields L max ≃ 230µm for V R = 0.2µ. Since this value can be of the order of or even smaller than the system size, we conclude that localization of BQPs in ultracold Bose gases is relevant for present-day experiments.

C. Validity of the leading-order result

Before turning to some numerical tests, let us review the conditions of validity of the results discussed above. The Born approximation for BQPs, which yields

Γ k ≃ Γ (2) k = [S(kξ)] 2 γ (2)
k , requires (i) the first-order smoothing solution [Eq. (10) with V replaced by V (1) as given by Eq. ( 18)], (ii) the first-order decoupling of the g + and g -modes that leads to Eq. (30), and (iii) the Born approximation γ k ≃ γ

(2) k to 7 For a complete discussion in the framework of the pure Schrödinger particles in speckle potentials, see for instance Refs. [32,33].

be valid. The weak disorder condition | V R | ≪ µ alone ensures (i) and (ii). Note that this criterion of weak disorder appears less stringent on the amplitude V R of the bare potential, since smoothing reduces the amplitude of V with respect to V . The regime of validity of the Lyapunov exponent derived for Schrödinger particles in a weak-disorder expansion is in itself a subtle issue, as the successive terms in the perturbation series all depend on the disorder amplitude and the energy of the particle. The resulting asymptotic series is wellbehaved in the high-energy limit. A precise inspection of the low-energy limit, where the terms of the series blow up, is necessary to exhibit a rigorous criterion for the validity of a truncated perturbation expansion (see e.g. Ref. [136]). For single particles, γ k ≪ k is usually retained [START_REF] Lifshits | Introduction to the Theory of Disordered Systems[END_REF]. In physical terms, the localization length should exceed the typical wavelength of the particle. This sets a V R -dependent lower bound on the single-particle energies for which the perturbative result is meaningful. Translating the above criterion to BQPs (Γ k ≪ k), we obtain

|V R | µ σ R ξ ĉ2 (2kσ R ) ≪ (kξ) 3/2 + 1 2(kξ) 1/2 . ( 50
)
This condition of validity resembles the corresponding one for Schrödinger particles:

(|V R |/µ) (σ R /ξ) ĉ2 (2kσ R ) ≪ (kξ) 3/2 .
As expected, the two coincide in the FP regime (k ≫ 1/ξ). However, they differ significantly in the PH regime (k ≪ 1/ξ). Indeed, for free particles, perturbation theory always breaks down at low energy (i.e., for k → 0). Conversely, for BQPs in the PH regime, the strong screening of the disordered potential leads to a completely different condition:

(|V R |/µ) (σ R /ξ) ĉ2 (2kσ R ) ≪ 1/(kξ) 1/2
. The latter is always valid at low energy, with the assumption that ĉ2 (2kσ R ) is of the order of unity at most. We thus find that the validity condition ( 50) is easily satisfied on the whole spectrum by a potential that is weak enough, i.e., for

(|V R |/µ) (σ R /ξ) ≪ 1.

D. Numerical calculations

In order to test the accuracy of our perturbative approach, we performed numerical calculations of the Lyapunov exponent of BQPs in a 1D speckle potential, for various ratios σ R /ξ. The first step consisted in determining the ground-state solution n c of the GPE (2), using a propagation scheme in imaginary time. As a precise determination of n c (z) is required for a correct calculation of the low-energy eigenmodes of the BdGEs ( 5) and ( 6), we compared the result of this procedure with the smoothing expansion including up to ten perturbation orders (see appendix A). The values of ∆ computed with the two methods for σ R = 3/2 ξ agreed within a relative difference δ∆/∆ of 0.3% for V R = 0.05µ, while the r.m.s. difference of the computed density profiles typically amounted to a few 10 -4 V R /g. Homogeneous Dirichlet boundary conditions were imposed in the calculations, i.e., both the density profile n c and the BQP components f + and f -were constrained to vanish at the system boundaries. 8 We used system sizes of the order of a few 10 5 healing lengths, so that the corrections to the equation of state µ = gn c solely due to kinetic terms at the boundaries of the system were negligible, 9 in particular as compared to the corrections introduced by the disordered potential, described in section III A 1.

In a second step, the profile n c (z) obtained as described above was included into the exact BdGEs ( 5) and ( 6), and the Bogoliubov modes were computed by solving the associated discretized eigenvalue problem. Eigenvalues and eigenvectors of the resulting large, non-Hermitian band matrices were obtained for a limited set of target energies using standard ARPACK routines. 10 The Lyapunov exponent Γ k was then obtained by computing estimators of

-lim z→∞ ln[r(z)/r(z 0 )] |z -z 0 | , (51)
where r = (g + ) 2 + (∂ z g + /k) 2 defines the envelope of g + (see e.g. Refs. [33,[START_REF] Lifshits | Introduction to the Theory of Disordered Systems[END_REF]), and the abscissa z 0 refers to the localization center of each eigenstate. In order to obtain accurate estimates of the average value of the logarithm of r(z) at infinity for a wide range of σ R /ξ and kξ parameters, the numerical calculations were carried out in a large box of size L = 3.2 × 10 5 ξ (i.e., 10 5 -10 6 σ R ), and disorder averaging was performed over 200 randomly generated 1D speckle patterns. We checked that, as expected, using r(z) |g + (z)| in expression (51) produces the results but requires more extensive disorder averaging, owing to the divergence of ln |g + (z)/g + (z 0 )| at the nodes of g + . We also found that replacing g + by one of the functions f + , f -or g - in the expression of r leaves our numerical estimates of the Lyapunov exponent (51) unchanged, within a relative difference of less than 5% for all the parameters in this study.

On the whole spectral range which spans the phonon and free-particle regimes, the numerical data (filled points) shown in Fig. 6 are in excellent qualitative and fair quantitative agreement with the analytic prediction of Eq. ( 46) (thick solid lines), which corresponds to second-order perturbation theory. The choice of the value of V R /µ was motived by experimental relevance and, as aforesaid, by numerical tractability for an entire set of σ R /ξ and kξ parameters.

Still, for this intermediate value of V R /µ, we find slight deviations between the numerical data and Eq. ( 46). As discussed in Sec. IV E, this small difference can be attributed 8 Periodic boundary conditions were found to complicate the analysis of the asymptotic localization properties of BQP modes, due to the periodic regrowths of the localized states in space. 9 In a 1D box of length L, the relative correction to the equation of state associated with the Dirichlet boundary conditions scales as 2E L /gnc ≃ 2ξ/L for gnc/E L ≫ 1, where E L = 2 /2mL 2 . 10 This approach was preferred over other traditional methods for boundaryvalue problems, such as shooting algorithms [START_REF] Press | Numerical Recipes[END_REF], since the propagation of initial data by a differential operator ξ∂z diverges exponentially on the length scale of a few healing lengths due to the coupling into evanescent modes in the case of Bogoliubov excitations (see appendix B).

k given by Eq. ( 46). The thick purple dashed lines correspond to the next-order expansion

Γ k ≃ Γ (2) k +Γ (3) k (see section IV E). The term Γ (3) k scales as V 3
R , and results from the addition of three contributions:

Γ (3) V k ,Vn (green dash- dotted), Γ (3)
V k ,V -(red dashed), and Γ

V k ,V k ,V k (blue dotted). The dots are the numerical data obtained with the procedure described in section IV D. mostly to third-order terms which contribute to the exact Lyapunov exponent Γ k . Hence, the observed deviations are expected to be negligible for lower values of V R /µ. We also find that, even if these deviations cannot be completely neglected for V R ≃ 0.05µ, they do not change the qualitative behavior of Γ k . Finally, note that for comparable parameters the deviations appear smaller in the data of Ref. [54] compared to those of the present work. The present results are actually more accurate as only a lowest-order smoothing expansion was used to compute the density n c in Ref. [54].

the Bogoliubov-de Gennes equations (BdGE), which govern the (quasi-)condensate background and elementary excitations (Bogoliubov quasiparticles) of the Bose gas, respectively. In the first step, we calculate the mean-field density profile using a perturbative expansion of the GPE in V /µ. In the second step, the result is incorporated into the BdGEs. Turning to an adapted basis for the Bogoliubov wavefunctions, we then show that the BdGEs can be approximately mapped onto a Schrödinger-like equation, with an effective potential V k which depends on the bare potential V , the condensate density background, and the quasiparticle wave vector k.

Our approach is well suited to study the effects of disorder in interacting Bose gases and, in particular, to examine the Anderson localization of Bogoliubov quasiparticles. On the one hand, it applies to any kind of weak, correlated disordered potentials. We stress that it is not limited to i) Gaussian disorder, ii) (uncorrelated) white-noise potentials or iii) models of non-overlapping impurities. In particular, it applies to speckle potentials as used in many experiments with ultracold atoms. On the other hand, the only small parameter of the perturbative expansion is the ratio of the disorder amplitude V R over the chemical potential µ. Our approach differs in this respect from standard approximations used in various other works: (i) the approximation of the BdGEs by hydrodynamical equations, which confine the theory to excitations of typical wavelength λ ≫ ξ, where ξ is the healing length of the condensate [51,[START_REF] Castin | Proceedings of the Les Houches Summer School of Theoretical Physics[END_REF], (ii) the Thomas-Fermi approximation of the background density profile, which is valid only under the assumption ξ ≪ σ R , where σ R is the typical (minimal) length scale on which the external potential varies [51], (iii) the white-noise approximation which requires at least σ R ≪ λ [51,[START_REF] Giorgini | [END_REF]. Conversely, the approach developed here holds for any ordering of the length scales ξ, σ R and λ.

Although our approach can be used to describe the scattering and localization properties of Bogoliubov quasiparticles in any dimension, we focused in this work on the onedimensional case, which leads to the strongest localization effects [54]. In the low-energy limit and at the leading and next-leading orders in the disorder amplitude, we found a quadratic scaling of the Lyapunov exponent with the quasiparticle energy, Γ k ∼ ǫ 2 ∝ k 2 . This finding is consistent with known results on the localization of phonons [129][130][131][132] and studies in the white-noise limit [51]. The effective potential V k derived in our approach provides a physical interpretation of this suppression of localization in the low-energy phonon regime (kξ ≪ 1) in terms of an efficient screening of the longwavelength modulations of the external potential by the background density of the Bose gas. In the free-particle regime (kξ ≫ 1), the Lyapunov exponent asymptotically approaches the exponent of a bare Schrödinger particle, Γ k ≃ γ k . For uncorrelated potentials, the Lyapunov exponent of Bogoliubov quasiparticles thus falls off as Γ k ∼ 1/ǫ ∝ 1/k 2 in the highenergy limit. For correlated disorder, the high-energy decay of the Lyapunov exponent strongly depends on the large k behavior of the disorder power spectrum, Γ k ∼ C 2 (2k)/k 2 . If C 2 (2k) has a finite support as for speckle potentials for instance [18,27,33], effective mobility edges arise as for bare-particles [32,33,134,135].

Most importantly, our approach covers the crossover between the phonon and free particle regimes. We find that localization (Γ k ) is maximum at a given energy ǫ. For uncorrelated potentials, this maximum lies around E ≃ µ (i.e., kξ ≃ 1), i.e., at the crossover between the phonon and free-particle regimes [51,54]. For correlated potentials, the strength of localization is also determined by the detailed power spectrum of the potential and the energy of strongest localization depends on both the quasicondensate healing length ξ and the disorder correlation length σ R .

Finally, let us discuss some possible extensions of our work. On the experimental side, the observation of localized Bogoliubov excitations appears as a challenge. It would be interesting to search for evidence of localization in the broadening of the dynamic structure factor, as measured in Bragg spectroscopy experiments [139][140][141]. Such broadenings have been measured to characterize coherence lengths of a couple hundred micrometers in elongated quasicondensates [141]. We infer therefrom that Bragg spectroscopy should allow the measurement of localization lengths of the same order of magnitude. On the theoretical side, the localization properties of Bogoliubov quasiparticles in two and three dimensions are expected to exhibit an even richer phenomenology. In particular, as pointed out previously [72,142], the screening of disorder by interactions is expected to lead to the possible occurrence of two mobility edges in three dimensions. In such a scenario, a first delocalization transition would occur at high energy, as for bare particles, and a second one would occur at low energies, as the effect of disorder is suppressed in the limit of vanishing quasiparticle energy [131]. The localized states would then reside around kξ ≃ 1. Below a critical amount of disorder however, no quasiparticle states should be localized at all. As the screened potential derived here accurately describes scattering in higher dimensions as well, it may offer a simple avenue for the description of such a phenomenology. Second, Eq. (A8) can be used to show by induction that, whenever the potential V is extended, then so are all the perturbation orders φ (n) and V (n) . In the disordered case, these simple considerations show how localization can be destroyed in a regime of weak interactions compatible with the meanfield approach.

Appendix B: Decoupling basis for Bogoliubov-de Gennes equations in weak potentials

In this appendix, we motivate the introduction of the functions g + and g -(see Sec. III B) to solve the BdGEs ( 5) and ( 6), and justify the use of the Schrödinger-like equation (30) for weak potentials. Throughout the paper, we assume that the system lies in a box whose dimensions eventually tend to infinity to emulate the continuum limit, and we impose periodic or homogeneous Dirichlet boundary conditions on the functions f + and f -. In the latter case, the density n c and the BQP components f + and f -vanish at the system boundaries. However, in the limit ξ/L → 0, where L is the system size, and the absence of an external potential, the system can be regarded as homogeneous. Together with these boundary conditions, the BdGEs form the differential problem to be solved. The system of coupled equations ( 21) and ( 22) can be rewritten as a differential problem in matrix form:

ξ 2 ∇ 2 F (r) = H ǫ (r)F (r), (B1)
where

F (r) = f + (r) f -(r) , (B2)
and H ǫ (r) = H (0)

ǫ +W (r) is a real-valued, symmetric matrix, which depends on the energy ǫ and the position r, with11

H (0) ǫ = 0 -ǫ/2µ -ǫ/2µ 1 (B3) W (r) = 1 2µ V + ∆ -V 0 0 V + 3∆ -3 V , (B4)
where the position dependence of V and V was dropped for conciseness. The two differential equations on f + and f - associated with Eq. (B1) are strongly coupled via the offdiagonal terms in H

ǫ . Since W is small (at most of first order in V R ), it is worth working in the basis that diagonalizes H (0) ǫ . Indeed, although the change of basis may introduce off-diagonal (coupling) terms in W , these terms will remain small. We will then show that this approach is suitable for the set-up of a perturbation expansion.

Bogoliubov-de Gennes equations in the decoupling basis

In the absence of an external potential (V = 0), the matrix W vanishes identically. Then, the matrix H ǫ = H

F k ∝ √ ρ +1/ √ ρ and F β ∝ -1/ √ ρ √ ρ , (B7)
respectively, where

ρ = µ ǫ + 1 + µ ǫ 2 . ( B8
)
For simplicity, let us restrict our discussion to the 1D case.12 Since Eq. (B1) is of second order, each eigen-subspace corresponds to two possible solutions of the BdGEs. First, the solutions corresponding to +β 2 ξ 2 are e ±βz F β , which either grow or decrease exponentially. These modes are the evanescent modes discussed below Eq. ( 28), which thus appear naturally in this formulation. Since no solution of Eq. (B1) in the subspace spanned by these eigenvectors can satisfy the boundary conditions on both boundaries, these modes are forbidden in the case of vanishing V . Second, the solutions corresponding to -k 2 ξ 2 are e ±ikz F k , which are oscillating (plane wave) modes. Since solutions of Eq. (B1) that satisfy the boundary conditions can be built with linear combinations of these eigenvectors, these modes are allowed for V = 0. They correspond to the well-known physical solutions of the BdGEs (5) and (6).

The procedure above allowed us to decouple the BdGEs in homogeneous space. Let us introduce now the external potential V (r). Rewriting Eq. (B1) in the eigenbasis of H (0) ǫ , we find

ξ 2 ∇ 2 G(r) = H ′ ǫ (r)G(r), (B9)
where G(r) = P -1 F (r), and

P -1 = + √ ρ +1/ √ ρ -1/ √ ρ + √ ρ (B10)
is the inverse of the transformation matrix from the F -basis to the G-basis. The term H ′ ǫ = H ′ ǫ (0) + W ′ , which is the analog of H ǫ in the G-basis, contains a homogeneous part

H ′ ǫ (0) = -k 2 ξ 2 0 0 β 2 ξ 2 , (B11)
These are the expressions reproduced in Eqs. ( 30) and (31), which form the basis of our approach to calculate the BQP modes in leading-order perturbation theory. The advantage of the g ± representation is that, at this level of approximation, the coupled equations ( 21) and ( 22) reduce to a simple closed equation for g + , the solution of which also determines g -via Eq. (B16).

Figure 2 .

 2 Figure 2. (color online) Comparison of the first-order smoothing solution [Eq.(10) with ∆ = 0 and V = V(1) ] and exact numerical computations of the density profile in the cases of (a,c) a 1D speckle potential with reduced autocorrelation function c(u) = sin(u) 2 /u 2 and correlation length σR = ξ and of (b,d) a 1D bichromatic periodic potential. The periodic potential reads V (z) = VR[cos(k0z) + cos(k0z/7)] with k0ξ = 1. In both cases, VR = 0.1µ. (a,b) The analytical expression gnc -V(1) (z) (broad grey line) is hardly distinguishable from gnc(z) as obtained numerically (red solid line). The average density gnc (black dashed line) and the Thomas-Fermi limit gnc -V (z) (green dotted line) are shown as references. (c,d) Corresponding power spectra of the modulations of gnc(z) (red solid line / red dots), V(1) (z) (broad grey line / black crosses), and V (z) (green dotted line / green circles). The grey dashed line in (d) is a plot of 1/[1 + (qξ) 2 ] 2 .

ε / µ = 10 Figure 3 .

 103 Figure 3. (color online) Plot of the screened potential V k (z) for the same 1D speckle potential (σR = ξ, VR/µ = 0.1) as in Fig. 2(a), with ǫ/µ = 0.1 (i.e., kξ ≃ 0.05) and ǫ/µ = 10 (i.e., kξ ≃ 2). The bare potential V (z) and the smoothed potential V (z) are shown for comparison.

Figure 4 .

 4 Figure 4. Screening function S. The dashed lines show the asymptotic behaviors in the phonon regime [S(kξ) ≃ 2(kξ) 2 for kξ ≪ 1 or, equivalently, ǫ ≪ µ] and free-particle regime [S(kξ) ≃ 1 for kξ ≫ 1 or, equivalently, ǫ ≫ µ].

Figure 5 .

 5 Figure 5. Contour plot of the Lyapunov exponent of BQPs in a speckle potential, as given by Eq. (46). Beyond kσR = 1 (black dashed line), the Lyapunov exponent vanishes completely in the Born approximation, due to the finite support of the speckle power spectrum. The green solid line represents the wave vector of maximum localization for each ratio σR/ξ.

Figure 6 .

 6 Figure 6. Lyapunov exponent of BQPs in 1D speckle potentials with the statistical properties described in section IV B 2, with VR/µ = 0.05, and for a) σR = 3.7 ξ, b) σR = 3/2 ξ, and c) σR = 0.4 ξ. The thick solid lines correspond to the Born term Γ

The definition of a phase operator requires special care. A suitable definition can be found in Ref.[START_REF] Mora | [END_REF], where a lattice model is used for a rigorous formulation of the Bogoliubov-Popov theory for quasicondensates. The equations derived in this lattice model coincide with the continuous formulation of Eqs. (2), (5) and (6) upon replacement of the coupling constant g by an effective coupling constant that depends on the lattice spacing l, and converges to g in the limit l → 0 in 1D[START_REF] Mora | [END_REF]119

].2 We discard in the canonical form of K and in Eqs. (3) and (4) the contribution of the P and Q operators which arise in Bogoliubov approaches without particle number conservation[120,121], as these operators play no role in the elementary excitations of the Bose gas and vanish in numberconserving approaches. Note also that the orthogonal projection of the Bogoliubov modes, as used in conserving approaches[START_REF] Mora | [END_REF]121], does not alter the results presented here.

Throughout the paper, the Fourier transform is defined as f (q) = (2π) -d/2 drf (r)e -iq•r . The notation f is used for the Fourier transform of functions f with dimensionless arguments.

Note that the effective energy 2 k 2 ν /2m appearing in the Schrödinger-like equation(30) differs from the actual energy ǫν of the BQP under consideration, as shown by Eq. (27).

In an equivalent formulation, Hǫ can be regarded as an operator acting on two-vectors of functions f + and f -. The position representation of Eq. (B1) is adopted here for simplicity.

The conclusions are naturally extended in higher dimensions.

ACKNOWLEDGMENTS

We thank A. Aspect, P. Bouyer, D. Clément, T.-L. Dao, L. Fontanesi, C. Gaul, V. Gurarie, C. A. Müller, V. Savona, and M. Wouters for fruitful discussions. This research was supported by CNRS, the European Research Council (FP7/2007-2013 Grant Agreement No. 256294), Agence Nationale de la Recherche (Contract No. ANR-08-blan-0016-01), RTRA-Triangle de la Physique, and the Alexander von Humboldt Foundation (Fellowship No. 1139948). We acknowledge the use of the computing facility cluster GMPCS of the LUMAT federation (FR LUMAT 2764). Laboratoire Charles Fabry de l'Institut d'Optique (LCFIO) is a member of the Institut Francilien de Recherche sur les Atomes Froids (IFRAF).

E. Beyond the Born approximation

While the screened potential V k of Eq. ( 30) and the Lyapunov exponent Γ

(2) k of Eq. (41) accurately account for the scattering and localization properties of BQPs in the limit of weak potentials, going beyond the leading-order (Born) approximation used to derive them is interesting in several respects. First, it should allow us to address the question of localization beyond cutoffs in the Lyapunov exponent Γ (2) k such as the one arising in speckle potentials at wave number 1/σ R , described in section IV B 2. Second, studying the localization properties of BQPs for stronger disorder or weaker interactions is of particular importance, since with increasing V R /µ ratio the interacting Bose gas moves away from the deep superfluid (quasi-)BEC regime and into the weakly interacting fragmented Bose-glass phase [53,98,99,105,112]. While recent studies [58,105,112] have shown that the low-energy scaling of the inverse participation length of BQPs (which characterizes their short-range localization properties) is modified through the phase transition, the impact of strong disorder on the long-range localization properties of BQPs remains an open issue.

We now briefly address the localization properties of the Bogoliubov quasiparticles beyond the Born approximation by an inspection of the next-order terms of the weak-disorder expansion. While details of the derivation can be found in appendix C, here we just outline the approach. The starting point is again the set of BdGEs ( 25) and (26) in the decoupling basis g ± . Retaining terms up to second order in the potential amplitude V R , we find a new approximate Schrödinger-like equation for g + :

where V k is the screened potential of Eq. ( 30), proportional to V R . The two additional terms V n and V -scale as V 2 R , and account for the second-order modulations of the density profile n c and for the coupling of g + and g -, respectively. Their expressions are given in Eqs. (C3) and (C4). Although Eq. ( 52) is valid a priori only in the regime k ≪ min(1/σ R , 1/ξ) (see Appendix C), we found it to provide a rather good approximation over a wider range of parameters (see Fig. 6 andbelow). On the basis of the Schrödinger-like equation ( 52), we calculate the third-order contributions to the Lyapunov exponent following the approach of Ref. [33]. We find

where Γ

V k ,Vn (resp. Γ

V k ,V-) stems from the cross-correlator of V k with V n (resp. V -), and the last term involves the three-point autocorrelation function of V k [see Eqs. (C9) to (C11)]. These various contributions are plotted in Fig. 6 for a speckle potential and various σ R /ξ ratios. We find an excellent agreement between the numerical data and the analytical result

k . For the smallest value of σ R /ξ however, the third order term Γ k and the numerical data [see Fig. 6(c)]. This may be due to the fact that the criterion k ≪ 1/ξ is not met or to contributions of higher order in Γ k . Let us now discuss the properties of the terms appearing in Eq. (53).

Note first that the third-order contributions scale as V 3 R . They are thus non-zero only for potentials which do not possess symmetric statistics under the transformation V → -V . This holds for speckle potentials as used in experiments with ultracold atoms [33].

A generic feature of the contributions to Γ

k is that they are all of the form

where i is an index labeling any of the terms in Eq. ( 53), F i is some function, and ĉ3 (q, q ′ ) is the Fourier transform of the reduced three-point autocorrelation function of the bare potential V . Thus, we find that if ĉ3 has a compact support, Γ

k vanishes over an extended part of the spectrum, as in the single-particle case [32,33]. In particular, in the case of speckle potentials, one finds a high-momentum cutoff that is identical to the cutoff in ĉ2 , whereby

k does. This result is consistent with the argument that no odd power of V R can provide the leadingorder term in a given part of the spectrum, since the Lyapunov exponent is a non-negative quantity [33]. Note that this feature emerges from the analysis of the terms in Eq. ( 53) although the region k ≈ 1/σ R a priori lies outside the regime of validity of Eq. (52).

Figure 6 shows that the magnitude of Γ

V k ,V-, which originates from the coupling of g + and g -, is small compared to the other terms in Γ

(3) k . This further legitimizes the use of the g ± basis and suggests that the difference between the analytical and the numerical results in Fig. 6(c) are likely to be due to higher order terms rather than the violation of criterion k ≪ ξ. We also note that, remarkably, the contributions

to constant values at zero energy which turn out to be opposite and thus cancel out in the calculation of Γ (3) k . This cancellation between two terms that seem to have different origins in the perturbation expansion is certainly not accidental, and must be due to the g ± representation chosen to set it up. Finally, an inspection of the low-energy limit of expression (53) in the case of the above speckle potential yields the scaling

a feature which is likely to be generic for disordered potentials with ĉ3 (0, 0) = 0.

V. CONCLUSION

In this work, we presented a general approach, valid in any dimension, to describe a weakly interacting Bose gas of chemical potential µ subjected to a weak inhomogeneous potential V (r). This approach relies on a two-step perturbative expansion of the Gross-Pitaevskii equation (GPE) and

APPENDICES Appendix A: Perturbation series of the smoothing solution

In this appendix, we calculate explicitly the leading terms of the modulations V (r) = gn c -gn c (r) of the mean-field interaction term gn c (r) and the deviation ∆ = gn c -µ from the homogeneous equation of state. We consider the limit of a weak external potential V and strong interactions which is relevant to our study (see Sec. III A). Working along the lines of Ref. [104], the weakness of V suggests an expansion of the square root of the density in powers of V R /µ, which we write as

A1) where φ (0) (r) = 1 is the solution in the absence of disorder [START_REF] Pitaevskii | Bose-Einstein Condensation[END_REF] and the functions φ (n) (r) are real-valued. We thus have

and the quantities of interest at any order in the expansion series are readily obtained by collecting the terms at the corresponding order:

and

for n ≥ 1.

The functions φ (n) (r) are determined by inserting the perturbation series (A1) into the GPE (2), which is equivalently written as

and by collecting the terms of order n. The explicit calculation of φ (n) becomes increasingly tedious as the number of terms involved grows like n . Yet, the above perturbation hierarchy produces a simple recursion formula which can be used in analytical or numerical calculations. Following the procedure outlined above, for all n ≥ 1, we obtain

where G ξ (r) is the Green function associated with the operator -ξ 2 ∇ 2 + 1, which is best written in Fourier space as

and the convolution product is defined as

Applying the recursive procedure up to second order, we find

Then, using Eqs. (A4) and (A6), we find

and

Finally, let us make two comments on the above perturbative solution of the GPE. First, for the perturbation expansion be valid, the mean-field density profile n c must be weakly perturbed around the homogeneous value µ/g. While the original small parameter is V R /µ, Eq. (A8) shows that

is a somewhat looser criterion for the successive terms of the expansion to be small. The r.m.s. amplitude | V

R | [see Eq. (A19)] can be calculated explicitly from Eq. (A15):

(A20) and a potential-dependent part

Note that the BQP energy ǫ is here embedded in the dependence of k, β and ρ on ǫ [see Eqs. (B5), (B6) and (B8)].

For the sake of clarity, let us write explicitly the two coupled equations associated to Eq. (B9) in terms of the components g ± of G(r) = (g + (r), g -(r))

T :

where

These equations are equivalent to the BdGEs, without any approximation. The benefit of the transformation we have used is that, for a weak external potential, the terms appearing in W ′ are all small. In particular, the terms coupling g + and g -in Eqs. (B13) and (B14) are at most of the order of V (as

. The G-basis thus offers a suitable starting point, which takes into account the full structure of the BdGEs, and which allows for a perturbative approach in the regime of weak disorder.

Leading-order terms: mapping the Bogoliubov-de Gennes equations onto a Schrödinger-like equation

Let us now develop the perturbation expansion of the BdGEs in the G-basis. Since Eqs. (B13) and (B14) are weakly coupled, we can resort to the following self-consistent approach. Assuming that g -is vanishingly small compared to g + for small V R , we neglect the third term on the right-hand side of Eq. (B14). Then, solving for g -and retaining only the leading-order term in V R , we obtain 13

1+(|q|/β) 2 is the Green function associated with the differential operator -(1/β) 2 ∇ 2 + 1, written in 13 Note that ∆ is of second order in VR (see Sec. III A or Appendix A). Figure 7. (color online) BQP mode in the g ± basis. This eigenmode was computed at energy ǫ ≃ 1.1µ for a speckle potential with VR = 0.05µ and σR/ξ = 3/2. The mode is displayed here over a few healing lengths, while the size of the box used for the calculation is L = 3.2×10 5 ξ. The g + and g -components obtained numerically are given by the black dashed line and the blue solid line, respectively. The thick red line represents the convolution formula (B16), and the green dotted line the somewhat cruder approximation (C1) to g -. Note that the g + component has been rescaled by VR/µ in this figure for better comparison with the various representations of g -. With the above parameters, V

Fourier space. The positive smoothing function G 1/β satisfies dr ′ G 1/β (r ′ ) = 1, and decays on the length scale 1/β, which is smaller than the healing length ξ and than 1/k, i.e., the typical length scale over which g + varies (see Fig. 7). Thus, owing to the fact that 2m/( 2 β 2 ) < 1/µ and 2ρ/(1 + ρ 2 ) < 1, we can safely write

Then, in terms of orders of magnitude,

which is consistent with our initial assumption, i.e., g -is small compared to g + . In Fig. 7 we show numerical results which corroborate expression (B16) and the fact that g -is a term of order V

R /µ at most to g + . From the upper bound (B18), we infer that the second term on the right hand side of Eq. (B13) is of the order of ( V 2 R /µ)|g + | at most, while the third term contains terms scaling as |V R g + |. Hence, we neglect the former contribution, and obtain a closed equation for g + which is valid up to first order in V R /µ:

where

Appendix C: Lyapunov exponent of Bogoliubov quasiparticles beyond the Born approximation

In this appendix, the perturbation expansion of the Lyapunov exponent of BQPs is extended one order beyond the leading-order (Born) approximation, so as to include terms scaling as V 3 R . To this aim, we consider explicitly the terms in Eq. (B13) that are of second order in V R .

The third term on the right hand side of Eq. (B13) contains both an inhomogeneous term V (2) and an offset ∆ (2) that are proportional to V 2 R , and for which explicit expressions are given in appendix A. Elements of order V 2 R are also introduced into Eq. (B13) by the cross-term V g -, where g -may be replaced by expression (B16). If g + varies on a length scale larger than that of the other quantities in the integrand of Eq. (B16), we can use the approximation

While Eq. (C1) is justified for k ≪ min(1/σ R , 1/ξ), we found that it is actually a good approximation on a broader range of parameters (see Sec. IV E). For instance, this approximation shows good agreement with direct numerical results for g -for the parameters of Fig. 7. Hence, we get a new closed equation for g + , which now comprises all the terms up to order V 2 R and is legitimate in the low-energy limit:

) where V k is the screened potential (B20), and the terms V n and V -are potentials proportional to V 2 R :

The potential term V n follows from a second-order expansion of the ground-state density profile, and V -originates from the coupling between g + and g -. Both V n and V -have a nonvanishing average. These non-vanishing averages suggest an evaluation of the correlation functions at a wave vector off the energy shell (29) in the fourth-order Lyapunov exponent Γ (4) . However, these averages play no role in the correlation functions contributing to Γ

(3) k (see below), and can thus be disregarded at this level of approximation.

For Schrödinger particles of energy E = 2 k 2 /2m in a 1D disordered potential V , the leading orders of the weakdisorder expansion of the Lyapunov exponent read [33]:

where P denotes a Cauchy principal value, C 3 (z, z ′ ) = V (0)V (z)V (z ′ ) is the three-point correlation function and C 3 (q, q ′ ) is its Fourier transform. Replacing V in these formulas by the sum of the potential terms appearing in Eq. (C2) and collecting the different terms according to their order in V R , up to V 3 R , we find

V k ,Vn + Γ

V k ,V-+ Γ

(3)

where

corresponds to result (37), obtained in the Born approximation, and the remaining terms are the three contributions to the third-order Lyapunov exponent of BQPs:

Note that, since V k = 0, the non-vanishing mean values of V n and V -play no role in the various contributions to Γ (3) , and can substracted from V n and V -in Eqs. (C9) and (C11).

As the expressions of contributions (C9) to (C11) in Fourier space are quite involved, we do not reproduce them here. We refer to Fig. 6 and to section IV E, which provide a discussion of the behavior of these terms.