Quantum transport of atomic matterwaves in anisotropic 2D and 3D disorder

The macroscopic transport properties in a disordered potential, namely diffusion and weak/strong localization, closely depend on the microscopic and statistical properties of the disorder itself. This dependence is rich of counter-intuitive consequences. It can be particularly exploited in matter wave experiments, where the disordered potential can be tailored and controlled, and anisotropies are naturally present. In this work, we apply a perturbative microscopic transport theory and the self-consistent theory of Anderson localization to study the transport properties of ultracold atoms in anisotropic 2D and 3D speckle potentials. In particular, we discuss the anisotropy of single-scattering, diffusion and localization. We also calculate a disorder-induced shift of the energy states and propose a method to include it, which amounts to renormalize energies in the standard on-shell approximation. We show that the renormalization of energies strongly affects the prediction for the 3D localization threshold (mobility edge). We illustrate the theoretical findings with examples which are revelant for current matter wave experiments, where the disorder is created with a laser speckle. This paper provides a guideline for future experiments aiming at the precise location of the 3D mobility edge and study of anisotropic diffusion and localization effects in 2D and 3D.


Introduction
Transport in disordered media is a fascinatingly rich field, which sparks a broad range of phenomena such as Brownian motion [1], electronic conductivity [2,3], superconductivity [4], superfluid flows of 4 He on Vycor substrates [5], as well as localization of classical (electromagnetic or sound) waves in dense media [6,7] and of ultracold atoms in controlled disorder [8,9,10,11,12]. In the case of a matter particle for instance, two regimes should be distinguished. In the classical regime, where the de Broglie wavelength is vanishingly small, transport leads to normal or anomalous diffusion [13,14]. The dynamics is characterized by the appearance of a percolation transition, which separates a trapping regime -where the particle is bound in deep potential wells -from a diffusion regime -where the particle trajectory is spatially unbounded [15,16]. In the quantum regime, the wave nature of the particle determines its transport properties, in close analogy with those of a classical wave [17,18]. In this case, interference effects can survive disorder averaging, leading to striking effects such as weak localization [6], the related coherent back-scattering effect [19], and strong (Anderson) localization [20,21,22].
Localization shows a widely universal behaviour [23], but observable features significantly depend on the details of the system. It shows a renewed interest in the context of ultracold matter waves [8,9,10,11,12]. On the one hand, the microscopic parameters in these systems are precisely known and, in many cases, tunable, which paves the way to unprecedented direct comparison between experiments and theory [24,25]. This is a great advantage of ultracold atoms, compared to traditional condensed-matter systems. On the other hand, these systems offer new situations, which can induce original effects [26] and provide new test-grounds in non-standard disorder [27,28,29,30,31]. Major advances in this field were the observation of one-dimensional (1D) Anderson localization of matterwaves [32,33] and studies of the effects of weak [34,35,36,37,38,39,40,41,42,43] and strong [44,45,46,47] interactions in disordered gases. Presently, a major challenge is the study of quantum transport in dimensions higher than one. While localization is the dominant effect in one dimension [48,49], higher dimensions show a richer phenomenology where regimes of diffusion, weak localization and Anderson localization can appear [23]. Recent experiments reported the observation of an Anderson transition in momentum space using cold-atom kick-rotor setups [50,51,52], study of classical diffusion in twodimensional (2D) speckle potentials [53,54], coherent back-scattering [55,56], and evidence of Anderson localization in noninteracting Fermi [57] and Bose [58] gases in three-dimensional (3D) speckle potentials.
From a theoretical viewpoint, diffusion and localization of noninteracting matter waves have been thoroughly studied for disordered potentials with zero-range correlations [59,60] and isotropic correlation functions [61,62,63,64,65,66]. However, transport experiments in dimensions higher than one are most often performed with speckle potentials which are anisotropic, either effectively in 2D setups [53,54], or for fundamental optical constraints in 3D [57,58]. Moreover, correlations in speckle potentials can be tailored in a broad range of configurations [67], which offers scope for investigation of localization in nonstandard models of disorder [29,30]. Taking into account anisotropic effects is of fundamental importance because they can strongly affect coherent transport and localization properties. This was demonstrated in various stretched media [68,69,70,71,72,73,74,75,76,77,78]. Optical disorder, relevant to ultracold-atom experiments [57,58], can show significantly more complex anisotropic correlation functions, the effect of which has been addressed only recently [79].
In this paper, we study quantum transport and Anderson localization of matter waves in 2D and 3D anisotropic speckle potentials. We first introduce the basics of quantum transport of matter waves in disordered media (Sec. 2) and the models of disorder we focus on in 2D and 3D (Sec. 3). We then present a detailed description of the theoretical framework pioneered in Refs. [75,80], which intends to be pedagogical. We study single-scattering (Sec. 4), Boltzmann diffusion (Sec. 5), and localization (Sec. 6), as a function of the particle energy, and discuss in particular the different anisotropies of these quantities. From a technical viewpoint, while the scattering allows for analytic expressions as for isotropic models of disorder [62], diffusion and localization are more involved and require in general numerical diagonalization of a certain operator. Some analytic expressions are however found in some limits for anisotropic disorder. In Secs. 4, 5 and 6, we focus on the 2D case, which contains most of the anisotropy effects discussed in the paper. The 3D cases are discussed in the next sections, where we study the same quantities as above (Sec. 7). We also show that energy-dependent quantities calculated in the usual on-shell approximation should be renormalized in strong disorder, and propose a method to do it. It does not strongly alter the overall energy-dependence of the quantities calculated in the previous sections, and in particular their anisotropies. However, it may be important when comparing to energy-resolved experimental measurements. Most importantly, we show that it strongly affects the calculation of the 3D mobility edge. Finally, we summarize our results and discuss their impact on recent and future experiments on ultra-cold atoms in speckle potentials in the conclusion (Sec. 8).

Basics of quantum transport
Before turning to a more formal description, it is worth recalling the basic picture of coherent transport in a disordered medium, which is genuinely understood in a microscopic approach [81,80]. Consider a wave of momentum k and velocity υ = ∇ k ǫ/ [ǫ(k) is the dispersion relation] propagating in a disordered medium. We assume for the moment that the medium is isotropic and will drop this assumption in the following sections. The wave propagation is governed by scattering from the random impurities. Three typical energy-dependent length scales can be identified, which characterize three basic effects induced by the disorder (see Fig. 1). First, single scattering from impurities depletes the k-wave states, which can be seen as quasiparticles in the disordered medium, with a finite life-time τ s (k). Single scattering hence defines the first length scale, namely the scattering mean-free path, l s = υτ s , which characterizes the typical length travelled by the wave before it loses the memory Figure 1. (Color online) Schematic view of the coherent transport of a matter wave in a disordered medium, with special emphasis on the characteristic length scales. The figure shows a trajectory of a particle (solid multicolor line) in a twodimensional disordered landscape (blue surface). Along its trajectory, the wave loses the memory of its phase (encoded in the various colors along the trajectory) on the characteristic length ls (scattering mean-free path). Multiple scattering then deflects the trajectory and the wave loses the memory of its direction on the characteristic length lB (transport mean-free path). Interference between the multiple-scattering paths can finally cancel diffusion (strong or Anderson localization). The wave then acquires an exponentially decaying probability profile (orange-green surface) of characteristic length Lloc (localization length). of its initial state, and primarily the memory of its initial phase. Then, multiple scattering defines the second length scale, namely the transport (Boltzmann) meanfree path, l B , which characterizes the typical length travelled by the wave before it loses the memory of its initial direction. In general, several scattering events are necessary to significantly deflect the trajectories so that l B ≥ l s . The two length scales are found to be equal only in the white-noise limit (if it exists), where the wavelength is smaller than the typical size of the impurities. In this case the scattering is isotropic and the wave loses the memory of its phase and initial propagation direction at the same time. Within the distance l B , the transport crosses over from ballistic to diffusive. The average squared size of the wavepacket increases linearly in time, [2,3]. Finally, diffusive transport allows the wave to return to its initial position via loop paths, and interference effects enter the game. Each loop can be traveled in one way or the other, which generates two multiple-scattering paths along which exactly the same phase is accumulated during the successive scattering events. This coherent effect holds for any specific realization of the disordered potential and thus survives disorder averaging. Moreover, since these two paths are in phase, it gives rise to a constructive interference of the matter wave, which significantly enhances its return probability. This effect induces coherent back-scattering and weak localization, which leads to diffusive transport with a reduced diffusion coefficient, D * < D B [6]. For strong enough disorder, the diffusion can completely cancel, an effect known as strong, or Anderson, localization [22]. Then, the probability distribution of the wave decays exponentially in space, hence defining the third characteristic length, L loc , the so-called localization length.
The picture above shows that localization relies on two characteristics of the medium: coherence along the multiple-scattering paths and return probablity to the origin. One then understands that the strength of localization should be governed by the interference parameter kl B [82] (since the more the coherence length exceeds the typical length of a loop path, the more significant interference terms are) and by the dimension of space d (since the return probability decreases when d increases). As a matter of fact, in 1D and 2D, any state is localized, although disorder correlations may lead to strong energy-dependence of the localization length [83,84,28,27]. In 1D, one finds that L loc ∼ l B so that diffusion is strictly absent. In 2D, one finds l B < L loc , and diffusion shows up at intermediate distances and times. In 3D, the return probability is finite and localization appears only for sufficiently low values of kl B . A mobility edge shows up for kl B ∼ 1, which separates localized states (for kl B 1) from diffusive states (for kl B 1) [23,85].
The microscopic description outlined above offers a comprehensive picture of transport and localization effects for coherent waves in disordered media. The next subsections give mathematical support to this picture within a formalism adapted to anisotropic disorder.

Green functions
Consider a quantum particle in a given homogeneous underlying medium and subjected to some static randomness. Its dynamics is governed by the singleparticle Hamiltonian H = H 0 + V (r), where H 0 is the disorder-free, translationinvariant, Hamiltonian of the underlying medium, and V (r) is the time-independent (conservative) disordered potential. Without loss of generality, the disordered potential can be assumed to be of zero statistical average b , V = 0. The evolution of the wave function between t 0 and t > t 0 is determined by the retarded singleparticle propagator G(t, t 0 ) ≡ exp[−iH(t − t 0 )/ ] Θ(t − t 0 ), where the Heaviside step function Θ(t − t 0 ) accounts for temporal ordering. In the energy domain c , G is the retarded Green operator where E is the particle energy. It is the solution of the equation since the first order term, G 0 V G 0 , vanishes It is convenient to represent this equation diagrammatically: where a plain line is a Green function (grey for G 0 and black for G), the vertices (black dots) are scattering events and the dashed lines recall that they are correlated. The Dyson equation [86] with Σ(E) the self energy, can be developped in powers of V thanks to Eq. (3) so as to determine Σ. The average Green function then reads If the disorder is homogeneous, i.e. if its statistical properties are translationinvariant [87], then the disorder-averaged Green function is diagonal in k-space d : where ǫ(k) is the dispersion relation associated to H 0 and d the space dimension. In addition, if the statistical properties of the disorder are isotropic, then G(E, k) ≡ G(E, |k|). This features an effective homogeneous (i.e. translation-invariant) medium, which contains all necessary information to determine the disorder average of any quantity linear in G. It is the case of the spectral function A(E, k) defined by [81]: It contains all the information about the spectrum of the disordered medium. Using The spectral function can be interpreted (up to a numerical factor) as the (normalized) probability density for an excitation of momentum k to have energy E and dE 2π A(E, k) = 1. It is also the unnormalized probability, per unit energy, to find a particle of energy E with momentum k and is the density of states per unit volume. For a particle in disorder-free space, it is given by In the presence of disorder, Eqs. (7) and (9) yield (Color online) Schematic representation of the spectral function A(E, k) of a particle of energy E = 2 k 2 E /2m, as a function of the particle momentum k. The vertical red line is the spectral function for the disorder free particle A0(E, k) = 2πδ [E − ǫ(k)] with ǫ(k) = 2 k 2 /2m. In the presence of disorder the spectral function is shifted and broadened (black line).
with Σ ′ and Σ ′′ the real and imaginary parts of Σ, respectively. As represented schematically in Fig. 2, for a particle in free space [ǫ(k) = 2 k 2 /2m, where m is the mass of the particle] with a weak disordered potential, the spectral function has a Lorentzian-like shape as a function of k. It is centered in k 0 , solution of E − ǫ(k 0 ) − Σ ′ (E, k 0 ) = 0. The quantity Σ ′ (E, k 0 ) thus describes the shift in energy of the free-particle modes when they are dressed by the disorder. The quantity Σ ′′ (E, k) is the energy width of the spectral function, which defines the scattering mean free time or equivalently the scattering mean free path l s (E, k) = |υ|τ s (E, k). It accounts for the depletion of the free particle mode at E = ǫ(k) due to scattering from the disordered medium. The spectral function, which contains all the information about the relative weight, the energy, and the life time of the quasi-particles, will be the key ingredient in the following calculations. In addition, in ultracold atomic systems, a broad range of energies are involved, but only the momentum distribution is usually measured by time-of-flight techniques. The spectral function relates the energy distribution (D E ) and the momentum distribution (D k ) of the stationary particles in the disorder via which is normalized by dE 2π D E (E) = 1. The exact calculation of the spectral function requires the knowledge of the real and imaginary parts of the self energy Σ [see Eq. (10)], or, according to Eq. (8), the direct diagonalization of the disordered Hamiltonian and an average over disorder realizations. This is, in general, a complicated task, especially in dimensions larger than one and for anisotropic disorder. In Secs. 5 to 7, we work within the usual on-shell approximation [80], in which one neglects the real-part of the self energy Σ ′ (E, k) and the structure of the spectral function (see schematic dashed blue line in Fig. 2). In Sec. 7.4, we describe a method to go beyond the on-shell approximation, which amounts to renormalizing the energies in a self-consistent way [79].

Propagation of the Wigner function
Some quantities are not simply related to the averaged Green function G and require a more elaborate treatment. It is for instance the case of the spatial density and the momentum distribution. More generally, consider the time evolution of the one-body density matrix ρ(t) [81] or equivalently of the Wigner function [88] The spatial density probability is given by n(r, t) = dk (2π) d W (r, k, t) and the momentum distribution by D k (k, t) = dr W (r, k, t). It is fruitful to rewrite Eq. (13) in a form indicating explicitly the initial conditions, using the relation When averaging over the disorder, if there is no correlations between the initial state and the disorder, one finds [63] where W 0 (r, k) ≡ W (r, k, t 0 ) is the initial Wigner function and F k,k ′ (R; t) is the phase-space propagation kernel, defined by (if t > 0) and with k ± ≡ k ± q/2, k ′ ± ≡ k ′ ± q ′ /2, E ± ≡ E ± ω/2, and (q, ω) the Fourier conjugates of the space and time variables e . As discussed above, disorder averaging features a translational invariance in space so that Eq. (15) depends only on the difference R = r − r ′ . For the same reason, translational invariance, or equivalently momentum conservation, imposes that the sum of the in-going wavevectors (k + and k ′ − ) on one hand, and out-going wavevectors (k ′ + and k − ) on the other hand, are equal. It leads to the condition on momentum transfer, q = q ′ , in Eq. (16).
As can be seen in Eqs. (14) and (15), the building block to describe wave propagation in random media is the density propagator Φ, which can be represented as a four-point vertex with k ± and k ′ ± the left and right entries [see left-hand side e We use the Fourier transformf (q, ω) ≡ drdt f (r, t) exp[−i(q · r − ωt)].
of Eq. (18)]. The skeleton of this vertex is made by retarded and advanced Green functions (respectively G, represented by the top line, and G † , represented by the bottom line). It contains all possible correlations between the scattering events of these Green functions. Following the same approach as used for the average field propagator G [leading to the Dyson equation (5)], the vertex Φ = G ⊗ G † is formally constructed from the uncorrelated-average vertex G⊗G † . Without any approximation, Φ is then governed by the so-called Bethe-Salpeter equation (BSE) [81] represented diagrammatically as where U is the vertex function including all irreducible four-point scattering diagrams: The first term in the BSE (17)-(18) describes uncorrelated propagation of the field and its conjugate in the effective medium. The second term accounts for all correlations in the density propagation. Analogously to Eq. (6), the solution of the BSE (17)- (18) can be formally obtained from the inverse, if it exists, of the four-point operator Λ ≡ 1 − G ⊗ G † U f [89]: More explicitly, the (k, k ′ ) component of a four-point vertex Λ which fulfills momentum conservation is Λ k,k ′ (q, ω, E), such that k + , and Therefore Eq. (20) reads f In this context, the inverse of an operator Λ is defined by and can be expressed as a geometric series The operator Λ −1 (ω, E) can be expressed in terms of the eigenvectors and associated eigenvalues of the operator Λ(ω, E) which was used in Refs. [75,90] to solve the BSE. It then gives access, via Eq. (23) to Φ, which is the quantity of interest [see Eqs. (13) to (15)].
In the following we will see that the intensity kernel Φ has a diffusion pole, which takes the form where D is the so-called dynamic diffusion tensor. The average spatial density distribution is then given by where represents the initial joint positionenergy density and P (r − r ′ , t − t 0 |E) is the probability of quantum transport, i.e. the probability for a particle of energy E originating from point r ′ at time t 0 to be in r at t. It can be expressed thanks to Eqs. (14), (15) and (25) as the space-time Fourier Transform of the diffusion pole 1/[i ω − q·D(ω, E)·q]. In the long-time limit, we will encounter two different situations. First, if lim ω→0 D(ω, E) = D(E) is a real definite positive tensor, the diffusion pole of the intensity kernel (25) describes normal diffusion with the anisotropic diffusion tensor D(E), and the probability of quantum transport reads Second, if D(ω, E) ∼ 0 + − iωΛ(E) in the limit ω → 0 + with Λ(E) a real positive definite tensor, the pole describes localization. It leads to exponentially localized phase-space propagation kernel and probability of quantum transport at long distance. In 2D, where K 0 is the modified Bessel function, and in 3D, In both 2D and 3D, the fonction P (R) decays exponentially g over the characteristic length L u loc (E) along the eigenaxis u of the localization tensor L loc (E) ≡ Λ(E).

Conductivity and Einstein's relation
Finally, another quantity of interest for our problem -in parallel of those studied in sections 2.3 and 2.4 -is the conductivity. In complete analogy to the usual conductivity of charge in condensed-matter systems [86], we here define the conductivity tensor σ in our system as proportional to the current-current correlation function, via the Kubo formula h [6]: where Thanks to Einstein's classical argument, it was realized that, at thermal equilibrium, in a gas submitted to a force, the diffusion and drift currents have to be equal. This relation holds in general for quantum systems in the linear response regime (see e.g. Ref. [81]). In particular, here we expect the DC conductivity and diffusion tensors to be proportional : σ(ω = 0) ∝ D. Calculating σ B (ω = 0) in the Boltzmann and Born approximations for anisotropic disorder permits us to find the proportionality factor (see details in appendix Appendix C.1), which in our system yields

Disorder correlation function
Having recalled the general theory of quantum transport in disordered media, we now specify the framework of our study. We will consider ultracold matter waves in speckle potentials as realized in several experiments [91,92,93,94,95,96,32,46,47,53,57,58]. In brief, a speckle pattern is created when a coherent light beam is shone through a diffusive plate and focused by an optical lens of focal distance f (see Fig. 3 and Ref. [97]). At each point of its surface, the diffusive plate imprints a random phase on the electric field. The resulting electric field in the right-hand side of the lens is then the summation of many complex independent random components, and is (Color online) Schematic of the apparatus used to create an optical speckle pattern. A laser beam is diffracted by a ground-glass plate diffuser of pupil function ID(ρ), where ρ ≡ (ρx, ρy) spans the diffuser, which imprints a random phase on the various light paths. The intensity field, I(r), observed in the focal plane of a converging lens, is a speckle pattern, which creates a disordered potential V (r) for the atoms. therefore a Gaussian random variable according to the central limit theorem. The potential acting on the atoms is proportional to the intensity pattern (i.e. the square modulus of the electric field). It is thus a spatially (non Gaussian) random variable. It is mainly characterized by the two-point correlation function C(r) = V (r)V (0). For a fine-grain diffuser, the two-point correlation function C(r) is determined by the pupil function I D (ρ) (i.e. the intensity pattern just after the diffusive plate) [97]. For Gaussian laser beams of waists w x,y and plates with homogeneous transmission i , we have I D (ρ x , ρ y ) = I 0 e −2(ρ 2 . For the configuration of Fig. 3, in the paraxial approximation, we find x,y and σ ⊥x,y = λ L f /πw x,y where λ L is the laser wavelength. Here x and y are the coordinates orthogonal to the propagation axis z, and z = 0 corresponds to the focal plane. We chose V R ≡ C(r = 0) as definition of the amplitude of the disorder.

Anisotropic Gaussian speckle (2D)
If the atoms are confined in a 2D geometry by a strong trapping potential along z centered on z = 0, they experience a disordered potential with correlation function i Here we assume that the diffuser covers the full area lit by the Gaussian beam. If it is not the case, a cut-off has to be introduced in the pupil function, which results into some oscillations in the wings of the correlation function. In experiments, if the diffusive plate is sufficiently large, this effect is small, and we disregard it in the following.
, with σ ⊥ = σ ⊥x and ξ = σ ⊥x /σ ⊥y the configuration anisotropy factor. The Fourier transform gives the power spectrumC Without loss of generality, we assume that ξ ≥ 1. When |k| ≪ σ ⊥ −1 ⊥ ξ and we recover the power spectrum of white noise disorder, the only relevant parameter being V 2 R σ ⊥x σ ⊥y . The power spectrum (35) is obtained by shining an anisotropic Gaussian beam on the diffusive plate. It also approximately holds in the case of Ref. [53] where a quasi-2D Bose gas of width l z is subjected to a speckle created by an isotropic Gaussian laser beam shone with an angle θ with respect to the plane of atoms, if l z ≪ σ ⊥ ≪ σ . In this case ξ ≃ 1/ sin θ (θ ≃ π/6 for the experiment of Ref. [53]).

Single speckle (3D)
In the 3D case, the disorder correlation function C(r) is given by Eq. (33) with w x = w y = w. The resulting speckle pattern has correlation lengths σ in the propagation axis (z) and σ ⊥ in the orthogonal plane (x, y). In general 4f > w, and C(r) is elongated along z . The corresponding disorder power spectrum reads where k ⊥ is the projection of k in the (k x , k y ) plane. It is isotropic in the (k x , k y ) plane but has a significantly different shape along the k z axis. In particular, it shows a strong algebraic divergence when k z = 0 and k 2 x + k 2 y → 0. It features absence of white-noise limit, which reflects the long-range correlations of the potential j . The consequences of this property, obtained in the paraxial approximation, will be further discussed in the following.

Single-scattering
We now focus on the first time scale introduced in Sec. 2.1: The scattering mean free time.  whereC(k) is the disorder power spectrum. Using Eq. (11) and the disorder-free Green function, we thus have

Scattering mean-free time
where represents the integration over the k-space shell defined by ǫ(k) = E. In the following we discuss anisotropic properties of the scattering time for the 2D case (the 3D cases are presented in Sec. 7.1).
the scattering time does not depend on the direction of the incoming wave vector k. In general, the scattering is however anisotropic, i.e. the probability that the particle acquires a direction k ′ depends on the direction of k ′ relative to k. Isotropic scattering is found only for δ-correlated disorder In the case of anisotropic disorder we are interested in, not only the scattering depends on the relative direction of k ′ and k, but it also depends on the direction of the incoming wave k.

Anisotropic Gaussian speckle (2D)
Let us consider the 2D anisotropic speckle potential of geometrical anisotropy factor ξ introduced in Sec. 3.1. ReplacingC(k) by Eq. (35) in Eq. (40) and using the disorder-free dispersion relation of the vacuum in Eq. (41), we obtain the scattering mean free time wherek ≡ k/|k| is the unit vector pointing in the direction of k, Ωk is the k-space solid angle, k E ≡ √ 2mE/ is the momentum associated to energy E in free space and E σ ⊥ ≡ 2 /mσ 2 ⊥ is the correlation energy of the disorder. The scattering time (42) is plotted in Fig. 4 as a function of energy along the two main axes, for |k| = k E and for a fixed geometrical anisotropy ξ = 4. Let us discuss some limiting cases and use the notation In the low-energy limit, k E σ ⊥ ≪ 1, we have which is displayed in Fig. 4 (left-hand-side black lines). In this limit the de Broglie wavelength of the particle (2π/k E ) exceeds the correlation lengths of the disorder (σ ⊥x and σ ⊥y ) and the speckle can be approximated by a white-noise (uncorrelated) disordered potential. Equation (35) and τ E,k is isotropic, constant, and it only depends on the product V 2 R σ ⊥x σ ⊥y (up to corrections of relative order E/E σ ⊥ ).
In the opposite, high-energy limit, k E σ ⊥ ≫ ξ, the de Broglie wavelength of the particle is much smaller than the smallest correlation length of the disorder and the particle behaves 'classically'. SinceC(k) has a wider extension in thek y direction than in thek x direction (for ξ > 1), there are more scattering channels for particles travelling along x so that τ E,kx < τ E,ky . More precisely, we find which is shown in Fig. 4 (right-hand-side black lines). In particular, we find that in the high-energy limit τ E,k ∝ √ E. It is also interesting to study the anisotropy factor of the scattering time which is shown in Fig. 5 as a function of E/E σ ⊥ and ξ. As already mentioned τ E,k is isotropic in the white-noise limit, so that ξ s ≃ 1 for k E σ ⊥ ≪ 1 (left-hand-side red line in Fig. 5). When increasing the energy, the scattering time first increases along the direction with the largest correlation length, i.e. the direction in whichC(k) is narrower (x for ξ > 1). Therefore, ξ s increases with E, for sufficiently small values of E/E σ ⊥ , and we have ξ s > 1. Using Eq. (43), an explicit calculation yields For k E σ ⊥ ≫ ξ, using Eq. (44), we obtain  43) and (47)].
which shows that the anisotropy factor of scattering is proportional to the inverse of the geometrical anisotropy (right-hand-side red line in Fig. 5). Note that the classical limit relation (47) is universal provided that the configuration anisotropy factor is well defined, i.e. that the disorder correlation function can be obtained by the anisotropic homothety of an isotropic one, C(x, y) = C iso (x, ξy). In this high-energy limit, ξ s < 1 (contrary to the low-energy limit case). Therefore, for any value of ξ, τ E,k exibits an inversion of anisotropy when the energy increases, typically at E ∼ E σ ⊥ .
As described in section 2.3 the scattering time is the width of the spectral function. It can be measured in a 2D experiment such as that of Ref. [53] by monitoring the momentum distribution of an almost energy-resolved wavepacket [66]. To illustrate this, a plot of the spectral function as a function of momentum and at fixed energy is shown in Fig. 6. In each directionk the spectral function peaks at 4τ E,k / and has a width proportional to 1/τ E,k . The anisotropy of the scattering time is revealed in the angle-dependence of both these quantities. It is more apparent in the angular dependence of the amplitude, which shows marked peaks. At low energy, the maxima are located on the k x axis, while at high energy, they are located on the k y axis, which signals inversion of the scattering anisotropy.

Boltzmann diffusion
We now turn to the behaviour of the spatial density in the incoherent diffusive regime, which is characterized by the Boltzmann diffusion tensor D B (E). We first give an explicit formula for the diffusion tensor, in the framework of the usual on-shell approximation, and then apply it to 2D disorder (3D cases are discussed in Sec. 7.2).

Solution of the Bethe-Salpeter equation
In the independent scattering (Boltzmann) and weak disorder (Born) approximation, only the first term in Eq. (19) is retained and the irreducible vertex function U equals the disorder structure factor [81]: or equivalently Then, incorporating Eq. (48)-(49) into the BSE (17)- (18) and expanding it in series of U, one finds where the diffuson Γ reduces to ladder diagrams: It describes an infinite series of independent scattering events, which leads to Drudelike diffusion.
In appendix Appendix A, explicit calculations are detailed. In brief, in the longtime (ω → 0) and large-distance (|q| → 0) limit the vertex Φ is the sum of a regular term and a singular term [75,90]: The regular part is given by where (22)] and φ n E,k (λ n E ) are the eigenvectors (eigenvalues) of an integral operator involving the disorder correlation function and k f E,k : The regular part contributes to the finite time and finite distance propagation of the density, which we disregard here. The existence of the singular part is a direct consequence of the Ward identity [98] which expresses the conservation of particle number, and which guarantees that one of the eigenvalues of Eq. (54) is equal to one, λ n=1 E = 1 (see appendix Appendix A). In the framework of the on-shell approximation, such that ǫ(k) = ǫ(k ′ ) = E, in the long time and large distance limit (|q|, ω) → 0, the vertex Φ is given by with N 0 (E) the disorder-free density of states, and is the disorder-free spectral function. Equation (55) shows that the vertex Φ is dominated by the diffusion pole (i ω − q·D B (E)·q) −1 .
In the anisotropic case, harmonics couple, and the φ n E,k are no longer cylindrical/spherical harmonics.

Anisotropic Gaussian speckle (2D)
Consider again the 2D anisotropic speckle potential of Sec. 3.1. The first step in the calculation of D B is to determine the eigenfunctions φ n E,k and the associated eigenvalues λ n E of Eq. (58). We solve Eq. (58) numerically, by a standard algorithm of diagonalization, with 2 9 = 512 points, regularly spaced on the k-space shell |k| = k E . The diffusion tensor is diagonal in the basis made by the symmetry axes of the correlation function (35): {û x ,û y }.
The eigenvalues and some eigenfunctions obtained numerically are shown in Fig. 7 for various values of E/E σ ⊥ . As discussed above, we find λ n=1 E = 1. For E ≪ E σ ⊥ , only the first term in the right-hand side of Eq. (57) contributes to the diffusion tensor since all λ n>1 E are vanishingly small. When the energy increases, the values of the coefficients λ n>1 E increase. It corresponds to an increase of the weight of the terms associated to the orbitals with n > 1 in Eq. (57), and a priori all the orbitals with n > 1 might have an increasing contribution. However, the symmetry properties of the functions φ n E,k cancel the contributions of most of them, and only the orbitals with n = 2 and 3 do contribute (see below).  In the low energy limit, one can develop Eq. (35) in powers of |k|. Up to order O(E 2 /ξ 4 E 2 σ ⊥ ), the first three eigenfunctions are given by: with eigenvalue λ 1 E = 1; where B 2 and B 3 are constant values that do not intervene in the following. In this limit the numerical results agree very well with the analytical findings (which for clarity are not shown on Fig. 7). In the very low energy limit, the disorder power spectrum becomes isotropic and constant, 3.1 and Eq. (35)]. The orbitals φ n E,k are thus proportional to the cylindrical harmonics, which are exact solutions of Eq. (58) in the isotropic case (see appendix Appendix B, and use the parametrizationk x = cos θ andk y = sin θ). In contrast to the isotropic case where the values of λ n E are degenerated in a given l-level, here we find that the degeneracy inside a l level is lifted for any anisotropy ξ = 1 [see the values of λ 2,3 E below Eqs. (60) and (61)]. When the energy further increases, the anisotropy plays a more important role and the harmonics are more and more distorted (see Fig. 7). However their topology remains the same, and in particular the number of nodal points and their positions are unchanged. In the following, we thus refer to Z ±1 l -like orbitals.
; a fit of the numerical data gives the prefactors D The inset shows the transport anisotropy factor Incorporating the values of λ n E , φ n E,k and τ E,k in Eq. (57), we can determine the Boltzmann diffusion tensor. Figure 8 shows the resulting eigencomponents of the diffusion tensor. In the low energy limit (E ≪ E σ ⊥ ), using Eqs. (43), (60) and (61), we find that the first term in the right-hand side of Eq. (57) gives the leading contribution to D B (E) (of order E/E σ ⊥ ). This contribution is isotropic owing to the isotropy of τ E,k at low energy and of the underlying medium. At very low energy, in the white-noise limit, we recover an isotropic diffusion tensor D is universal for 2D disorder in the white-noise limit (when it exists). The Z +1 1 -like orbital φ 2 E,k contributes to the next order of D x B and the and which are displayed on Fig. 8 (left-hand-side solid lines). When the energy increases, the anisotropy first comes from the anisotropic contribution of the scattering time τ E,k , and from the lift of the degeneracy between λ 2 E and λ 3 E . When the energy further increases, the harmonics are distorted, -but their symmetries (i.e. periodicity and parity) are preserved (see Fig. 7). Hence, for the same reasons as in the isotropic case (see appendix Appendix B) only the Z ±1 1 -like orbitals couple to υ in Eq. (57) and contribute to D B while the others don't. The associated λ n E increase (see Fig. 7), the weight of the second term in Eq. (57) increases, and the components of the diffusion tensor show a very different behavior in the large-E limit. For k E σ ⊥ ≫ ξ, we found τ E,k ∝ k E (see Sec. 4.2). In addition, we find numerically a weak topological change of the orbitals with energy for E/E σ ⊥ 10 2 . Therefore the evaluation of D B with E is mainly determined by the normalization condition [see formula below Eq. (57)], which yields φ n E,k ∝ 1/ √ k E . Then, assuming the scaling 1 − λ n E ∝ 1/E, also verified numerically, we obtain D u B (E) ∝ E 5/2 , which matches the numerical results (see dotted black lines in Fig. 8). This scaling is similar to that found for isotropic disorder [62]. As shown in Fig. 8, the change of slope between the low-and high-energy regimes is different in the two directions. For this reason, the anisotropy factor of the diffusion tensor, ξ B = D x B /D y B shows a nonmonotonous behaviour versus E, with a marked peak (see inset of Fig. 8).
The Boltzmann transport anisotropy factor ξ B is shown in Fig. 9 for various configuration anisotropies ξ. As it is well-known, the scattering and transport mean free times are different quantities in correlated disorder, due to angle-dependent scattering [99,100,81]. In particular, in the 2D speckle we consider, we do not find any inversion of the anisotropy of the diffusion, contrary to the scattering time, i.e. the component D x B (E) of the diffusion tensor is always larger than the component D y B (E). For large values of E/E σ ⊥ , the Boltzmann transport anisotropy ξ B reaches a constant value (see the inset of Fig. 8 for a cut at ξ = 4), which increases with the geometrical anisotropy ξ (see inset of Fig. 9). This asymptotic value is larger than ξ for small ξ and smaller for larger ξ. Therefore the anisotropy of the diffusion in the classical regime is not simply related to the configuration anisotropy.
The two distinct regimes found in the behaviour of D B and the non-trivial anisotropy effects make the Boltzmann diffusion regime in anisotropic 2D potentials very interesting for future experiments. Those properties could be probed by imaging directly the atoms in the 2D speckle (as in Ref. [53]) and controlling the width of the atomic energy distribution.

Weak and strong localization
We now consider interference effects, which lead to weak and strong localization. We first describe the quantum corrections (Sec. 6), then the self-consistent theory (Sec. 6.1), and apply it to the 2D speckle potential (Sec. 6.2). The 3D case, which follows the same route, is discussed in Sec. 7.3.
subsectionWeak localization correction We calculate corrections to Boltzmann diffusion by taking into account quantum interference terms between the multiple-scattering paths. Those interferences appear when the correlated scattering events do not occur in the same order in the propagation of the field and its conjuguate. This is diagrammatically translated into crossing correlation lines as in the second term of Eq. (19) for example. In the weak scattering regime only the two-point correlations are retained in the scattering diagrams and the leading scale-dependent corrections to the classical conductivity are given by the maximally crossed diagrams [101,75,90,6]: the cooperon [Eq. (64)] and the first two Hikami boxes [Eqs. (65) and (66)].
where the cooperon X is the sum of maximally crossed diagrams and is the renormalized vertex function (see appendix Appendix C.2). Using time-reversal invariance [23,102,98,81], the cooperon X can be expressed in terms of the diffuson Γ [defined in Eq. (51)] The diffusion pole carried by Γ in the limit (ω, q) → 0 leads to a divergence of X when ω, k + k ′ → 0. In appendix Appendix C.3 we translate those diagrams into equations, and show that Using Einstein's relation (32) we then obtain the dynamic diffusion tensor D * (ω, Note that the quantum corrections ∆D(ω, E) do not explicitly depend on the disorder [i.e. onC(k)], but only on the Boltzmann diffusion tensor D B (E) [75]. In other words, in this approach, Boltzmann incoherent diffusion sets a diffusing medium, which contains all necessary information to compute coherent terms l . In particular, it follows from Eq. (71) that the weak localization quantum correction tensor ∆D(ω, E) has the same eigenaxes and anisotropies as the Boltzmann diffusion tensor D B (E). Thus the anisotropy can be removed by rescaling distances along the transport eigenaxes u by D u B /D av B (i.e. momenta are rescaled by D av B /D u B ) with D av B ≡ det{D B } 1/d the geometric average of the Boltzmann diffusion constants. Since ∆D is always negative in the limit ω → 0 + , the weak localization correction features slower diffusion than the one obtained from incoherent diffusion. Equivalently, as long as the correction (71) is small, one can write which is the lowest-order term of a perturbative expansion of 1/D * (ω, E).

Strong localization
The quantum interference correction (71) has been derived perturbatively and is therefore valid as long as the correction itself is small. In order to extend this approach and eventually describe the localization regime where D * vanishes, Vollhardt and Wölfle [102,98] proposed to self-consistently replace D B (E) by the dynamic diffusion tensor D * (ω, E) in the right-hand side of Eq. (72). For isotropic scattering this procedure amounts to resumming more divergent diagrams than the cooperon (which contain a square of a diffusion pole), thus contributing to localization [98,80]. Generalizing this standard approach to anisotropic disorder yields l This property is a consequence of the on-shell approximation.
In dimension d ≥ 2 the integral in the right-hand side of Eq. (73) features ultraviolet divergence. Since the diffusive dynamics is relevant only on length scales larger than the Boltzmann mean free path l u B (E) ≡ d m/2E D u B (E) along each transport eigenaxis, we regularize this divergence by setting an upper ellipsoidal cut-off of radii 1/l u B in the integral domain m . It corresponds to an isotropic cut-off in the space rescaled according to the anisotropy factors of D B as described above.

Anisotropic Gaussian speckle (2D)
We now solve the self-consistent equation (73) for the inverse dynamic diffusion tensor in the 2D case. In the long time limit ω → 0 + , the unique solution of Eq. (73) is of the form D * (ω, E) ∼ 0 + − iωL 2 loc (E), where L loc (E) is a real positive definite tensor. As described in Sec. 2.4, it leads to the exponentially decreasing propagation kernel (28). Solving Eq. (73) then yields the anisotropic localization tensor, where l av B (E) ≡ d m/2E D av B (E). The eigenaxes of the localization tensor are thus the same as those of the Boltzmann diffusion tensor and its anisotropy factor is the square root of that of D B (E), i.e. ξ loc ≡ L x loc /L y loc = √ ξ B . We now apply the self-consistent theory to our running example: the 2D anisotropic speckle potential with correlation function (35). Including the results for the Boltzmann diffusion tensor D B (E) obtained in Sec. 5.2 into Eq. (74), we find the localization tensor L loc (E). Figure 10 presents the eigencomponents of L loc in its eigenbasis {û x ,û y } as a function of energy, for a configuration anisotropy of ξ = 4 and two different amplitudes of the disorder, V R /E σ ⊥ = 0.2 and 2. At low energy (E ≪ E σ ⊥ , V R , V 2 R /E σ ⊥ ), using Eqs. (62) and (63), we find where the upper sign holds for direction x, and the lower sign for direction y. Equation (75) corresponds to the solid black lines in Fig. 10. As D B is almost isotropic for E/E σ ⊥ 1 (see Fig. 8), L loc is also almost isotropic in the whole range presented in Fig. 10. Equation (75) describes an isotropic localization tensor with an anisotropic correction which is significant only if V R /E σ ⊥ ξ 3/2 / ξ 2 − 1 (≃ 2 for ξ = 4). At higher energy, when k E l av which is plotted as dotted black lines in Fig. 8. According to Eqs. (62) and (63) (retaining only the lowest-energy term), this regime appears for E/E σ ⊥ m Although somewhat arbitrary the factor unity between the cut-off radius and 1/l u B (E) is justified by the agreement we find with another approach in the isotropic case, provided that the real part of the self energy is included, see Sec. 7.4.2.  (π/2ξ)(V R /E σ ⊥ ) 2 . When ξ = 4 (as in Fig. 10), it gives E/E σ ⊥ 0.015 for V R /E σ ⊥ = 0.2 and E/E σ ⊥ 1.5 for V R /E σ ⊥ = 2. As predicted by the scaling theory of Anderson Localization [23] and explicitely seen in Eq. (76), the 2D localization length increases exponentially at large energy (hence the limited energy range in Fig. 10). Therefore measuring it experimentally with ultracold atoms [54,103,104] is very challenging and can be done in a rather narrow energy window, in which L loc is larger than the resolution of the imaging system (L res ) but smaller than the size of the sample (L sys ). This is illustrated for σ ⊥ = 0.25 µm on Fig. 10 by the grey dashed lines corresponding to L res ≃ 15 µm and L sys ≃ 2 mm, which are typical values extracted from Refs. [53,58].
One can finally note that 2D speckle potentials bear a classical percolation threshold at energy E p ≃ −V R /2 [54]. In the classical regime (1/k < σ ⊥x , σ ⊥y ), genuine Anderson localization has to be distinguished from classical trapping, which happens for E < E p . However, classical percolation is not relevant for the parameters used in Fig. 10. Indeed, for |V R | ≤ 2E σ ⊥ (as in the figure) and for E < E p , we have E |V R |/2 ≤ E σ ⊥ , so that kσ ⊥y ≤ kσ ⊥x = kσ ⊥ 1, which is not in the classical regime.

Three-dimensional anisotropic disorder
In this section we apply the formalism introduced in Secs. 4 to 6 to the 3D speckle potential of Sec. 3.2. We discuss single-scattering (Sec. 7.1), Boltzmann diffusion (Sec. 7.2) and localization (Sec. 7.3) properties and the position of the mobility edge (Sec. 7.4).

Single-scattering
Inserting Eqs. (36) and (37) into Eq. (40), we find the scattering mean free time which is shown in Fig. 11 for |k| = k E [we use the definition τ E,k ≡ τ s (E, k Ek )]. SinceC(k) is isotropic in the (k x , k y ) plane, τ E,k only depends on the polar angle θ between k andk z and not on the azimutal angle φ. We find that the scattering time is an increasing function of energy. It is also shorter for particles travelling along the z direction (τ E,kz < τ E,k ⊥ ) for all values of E. As for the 2D case analyzed in Sec. 4.2, this is due to the wider extension ofC(k) in the plane (k x , k y ), which offers more scattering channels to particles travelling along z. In contrast to the 2D speckle case however, τ E,k shows no inversion of anisotropy.
In the low energy limit (k E σ ⊥ ≪ 1), τ E,k converges to a constant value. In contrast to the 2D case, it signals the absence of a 3D white-noise limit n . This can be attributed to the strong anisotropic divergence ofC(k) when |k| → 0, which reflects the long-range correlations of the disorder (see Sec. 3.2). More precisely, for |k|σ ⊥ ≪ 1, n In the case of a white-noise limit in 3D, the scattering time is isotropic with the scaling τ E,k ∝ 1/ √ E (i.e. ls E,k is constant). This can be found by inserting a constantC(k) in Eq. (40).
we havec Replacing this expression into Eq. (77) we then find which is independent of E. Equation (79) is plotted as solid black lines on the lefthand side of Fig. 11. Note that τ E,k does not become strictly isotropic in this limit. However, the residual anisotropy of the scattering time, found from Eq. (79) and from the anisotropy ofc(k) in Eq. (78), is very small, and practically unobservable (τ E,k ⊥ /τ E,kz ≃ 1.002). When the energy increases, the scattering time in the (x, y) plane is the first to deviate significantly from the low-energy behaviour at E ∼ E σ (= 3 × 10 −2 E σ ⊥ for the parameters of Fig. 11), while the scattering time in the z direction increases only at E ∼ E σ ⊥ . This can be understood again by the narrower width of the power spectrumC(k) in the k z direction. In the high-energy limit (k E σ ⊥ ≫ 1) the k-space shell integral of Eq. (77), which is done on a sphere of radius k E containing the origin, can be reduced to integrating c 1sp on the plane which is tangent to the sphere at the origin. We then find In particular, we find τ E, (both shown as the right-hand-side solid black lines in Fig. 11). The anisotropy of the scattering then becomes significant for the parameters of Fig. 11, τ E,k ⊥ /τ E,kz = √ πσ /2σ ⊥ in this limit. The high-energy scaling τ E,k ∝ k E , which was also found in our 2D speckle, is quite universal: as long as the power spectrum is of finite integral in all the planes (lines in 2D) crossing the origin, the procedure described above can be applied to Eq. (40). Then τ E,k only depends on the dispersion relation ǫ(k) and, in particular, it is independent of the space dimension.

Boltzmann diffusion
The Boltzmann diffusion is obtained, as in the 2D case analyzed previously, by solving Eq. (58) numerically and incorporating the results into Eq. (57). For the diagonalization of the integral operator (58) we use 2 7 × 2 7 = 128 × 128 points regularly spaced on the k-space shell |k| = k E . We have studied several values of the configuration anisotropy σ /σ ⊥ , which all show the same behaviour discussed below.
The eigenvalues λ n E of Eq. (58) for different energies, as well as the topography of the eigenvectors of Eq. (58) that dominate D x B (bottom row), D y B (2 nd row), and D z B (3 rd row) are shown in Fig. 12 for σ /σ ⊥ = 5.8. Similarly as for the 2D case, we fond that λ n E decays from 1 to 0 when n increases, more sharply for low energy. The φ n E,k are topologically similar to the spherical harmonics at all energies, i.e. they show similar nodal surfaces, but the associated λ n E are not degenerated in a given l-like level. More precisely, due to the cylindrical symmetry of the power spectrum, the value of λ n E associated to the Y +m l -like and Y −m l -like orbitals are the same for a given m, but the degeneracy between the different values of |m| is lifted. Figure 13(a) shows the resulting eigencomponents of the diffusion tensor in the 3D case for σ /σ ⊥ = 5.8. It is isotropic in the (x, y) plane, because of the cylindricalinvariance of the correlation functionC(k) around the axisk z . For the same symmetry reasons as in the isotropic case (see appendix Appendix B) and as in the 2D case, only the p-level-like orbitals couple to υ. For k E σ ⊥ ≪ 1, we find that D x,y B is dominated by the first term in Eq. (57) and D z B by the Y 0 1 -like orbital (n = 2 at all energies). For k E σ ⊥ ≫ 1, the situation changes: while D z B is still dominated by the Y 0 1 -like orbital, D x B is now dominated by the Y +1 1 -like orbitals and D y B by the Y −1 1 -like orbitals (respectively n = 6 and 5 at E = 50E σ ⊥ in Fig. 12) with a contribution of the Y ±1 3 -like orbitals increasing with E. At high energy, we find that the nodal lines of the Y ±1 3like orbitals calculated numerically are displaced compared to the associated spherical harmonics. Therefore their contribution in Eq. (57) does not cancel out for symmetry reasons. Those properties explain the main features of D B .
Firstly, we find that the diffusion tensor is larger along axis z (D z B > D x,y B ) for all values of E [see Fig. 13(a)], and the anisotropy of D B is thus reversed with respect to that of τ E,k (we recall that we found τ E,kz < τ E,k ⊥ for any E, see Sec. 7.1). This is due to the fact that the (Y 0 1 -like) orbitals contributing to D z B are associated to values of λ n E larger than those contributing to D x,y B (in Fig. 12, the φ n E,k are numbered by decreasing eigenvalues).   Secondly,C(k) shows a strong anisotropic, infrared divergence in the paraxial approximation (see Secs. 3.2 and 7.1). Following-up with the scaling ofc 1sp (k), Eq. (78), used to show that τk ,E is independent of energy for k E σ ⊥ ≪ 1, and inserting it into Eq. (58) and the associated normalization, we find that λ n E does not depend on E, and φ n E,k is of the form ϕ n (k)/ √ k E . Then, all terms in Eq. (57) are topologically unchanged and scale as E at low energy. The anisotropy of D B thus persists down to arbitrary low values of E and D u B ∝ E, as observed in the left-hand side of Fig. 13(a) for k E σ ⊥ ≪ 1 (i.e. E ≪ E σ ⊥ ). This is another manifestation of the absence of whitenoise limit o .
Thirdly, for k E σ ⊥ ≫ 1, we found τ E,k ∝ √ E. Then, assuming weak topological change of the orbitals and the scaling 1 − λ n E ∝ 1/E (confirmed numerically), we get φ n E,k ∝ 1/k E and D u B (E) ∝ E 5/2 . This scaling is confirmed in Fig. 13(a) by fits to the data for E ≫ E σ ⊥ (right-hand-side dotted lines). Remarkably, in spite of the different contributing terms in Eq. (57) at low and high values of E, the transport anisotropy is nearly independent of E with D z B /D x,y B ≃ 10 [see inset of Fig. 13(a)]. We have repeated the same study for various values of the configuration anisotropy, ξ = σ /σ ⊥ . They all show a similar behaviour as a function of energy as reported in Fig. 13(a) for ξ = 5.8. In particular, we found the same scalings with energy and a diffusion anisotropy that is nearly independent of energy. In Fig. 13(b), we plot ξ B versus ξ for three values of the energy. We find that the diffusion anisotropy monotonously increases with the configuration anisotropy, as could be intuited. In order to guess a fitting function for ξ B , one may rely on a simplified model of random walk in an anisotropic lattice of anisotropy factor ξ. If the transition time is governed by the travel duration between two wells, one expect o A 3D white-noise limit would lead to the scaling D u B (E) ∝ √ E and an isotropic limit at low energy.
ξ B ∝ ξ. If it is governed by the trapping time, one expects ξ B ∝ ξ 2 . In our continuous model of disorder, the situation may be expected to be somehow intermediate. For our considered range of ξ, we find that the fit ξ B = 0.59ξ + 0.21ξ 2 reproduces well our results as shown in Fig. 13(b).

Localization
In order to analyze strong localization effects, we now solve the self-consistent equation (73) for the 3D case in the long time limit (ω → 0). A threshold energy E c (mobility edge) appears, solution of D av for ω → 0, where L loc (E) is a real positive definite tensor. It characterizes exponential localization within the propagation kernel (29) with the anisotropic localization tensor L loc (E). The localization tensor is diagonal in the same basis as the Boltzmann diffusion tensor D B . Explicitely, we have where L av loc = det{L loc (E)} 1/3 is the unique solution of L av For E > E c , D * (ω, E) converges to a real definite positive tensor when ω → 0. It describes anisotropic normal diffusive dynamics, characterized by the propagation kernel (27) where D(E) is replaced by the quantum-corrected diffusion tensor  83)]. Therefore, as for D B , for the 3D configuration, the anisotropy factors of L loc and D * are nearly independent of E. The localization anisotropy ξ loc = L z loc /L x,y loc is plotted versus the configuration anisotropy on Fig. 13(b). At low energy, using the scaling of D u B (E) obtained previously we . When E increases, L u loc (E) grows and finally diverges at E c . In the diffusive regime the quantum corrections are significant only close to E c , while for higher values of E, D * (E) ≃ D B (E). Therefore, in the high E limit we have D u * (E) ∝ (D u B /D av B )E 5/2 as found previously (see Sec. 7.2).

About the 3D mobility edge
The self-consistent approach used above is expected to fairly describe the quantum transport properties [75,80,62]. It gives some quantitative estimates consistent with numerical calculations [105] and experimental data [106,56]. However it has two main flaws.
On the one hand, it predicts that, just below the mobility edge, the localization length diverges as L u loc (E) ∝ (E c − E) −ν with ν = 1 and, just above the mobility edge E c , the corrected diffusion tensor increases as D u * (E) ∝ (E − E c ) s with s = 1. Those values of the critical exponents ν and s are consistent with the prediction s = ν(d − 2) of the scaling theory [23,107] and they are independent of the choice of cut-off that we made. However, it is known, from advanced numerical calculations on the Anderson model [108,109] and from experiments [52], that they are not correct. The correct value of the critical exponents in 3D is ν = s = 1.58 ± 0.01 [108,109]. In order to reproduce this value, it seems necessary to take into account the fractal nature of the wave functions at the critical point [110], which is beyond the self-consistent theory of AL.
On the other hand, in contrast to critical exponents, the mobility edge, E c is a non-universal quantity and should be determined from microscopic theory. In this respect, the on-shell approximation is questionnable because it neglects the strong modification of the spectral function induced by the disorder. This renormalizes energies and may thus strongly affect the value of E c . 7.4.1. Energy renormalization In order to improve the self-consistent method, one could in principle use the more sophisticated approach of Ref. [105], which does incorporate the spectral function, and provides values of E c in agreement with numerical calculations in the Anderson model. For continuous disorder, one may rely on the approach of Refs. [60,64], which has been applied to several standard models of disorder. However, since we are interested in continuous disordered potentials with fine and anisotropic structures, these methods are hardly practicable. From a numerical point of view, estimates of necessary ressources seem out of present-day possibilities. In order to overcome this issue, we have proposed in Ref. [79] an alternative method based on the assumption that the leading term missing in the on-shell approximation is the real part of the self energy, where P is the Cauchy principal value, see Eq. (39). A quasi-particle of momentum k has an energy E, solution of E − ǫ(k) − Σ ′ (E, k) = 0. Here, we incorporate Σ ′ (E, k) into the theory self-consistently and by averaging, in first approximation, its k-angle dependence. It amounts to replace the on-shell prescription by Within this approach, all previous quantities [τ s (k), D B , L loc , D * ] are now regarded as functions of E ′ instead of E. It does not change the overall energy dependence of the quantities discussed above, but may be important for direct comparison to energyresolved experimental measurements. In the following we concentrate on the 3D mobility edge E c . It is the solution of E c − ∆(E c ) = E ′ c , where E ′ c is determined using the on-shell approach and ∆ can be regarded as an energy shift, which renormalizes the energies.

Isotropic disorder
Here, we validate the above approach by a direct comparison to an alternative method applicable to isotropic disorder. Consider a   [64], that we reproduce here), with the on-shell method (E ′ c , red crosses) and with the renormalized selfconsistent approach (corrected Ec, thick blue circles), for an isotropic 3D speckle potential. When comparing to Fig. 8 of Ref. [64], note that in Ref. [64] the reference of energy is the minimum value of the disorder and that we have the correspondences E ξ = Eσ/2 and U = V 2 R .
speckle disorder obtained inside an integrating sphere lit with a laser beam, the realspace correlation function of which reads [62,64] with σ the correlation length. The associated power spectrum (see appendix Appendix B) is isotropic and bears the same infrared divergence as the anisotropic model of 3D disorder considered in this work as well as other configurations [79]:C(k) ∝ 1/|k| when |k| → 0. Figure 14 shows the on-shell mobility edge E ′ c calculated as in Sec. 7.3 (see also Ref. [62]), the renormalized mobility edge E c calculated by our method (see Sec. 7.4.1), and the mobility edge found using the selfconsistent Born approximation (SCBA) in Ref. [64]. As it is clearly seen in Fig. 14, the disorder-induced modification of the spectral function plays a major role for the prediction of the mobility edge. While the on-shell mobility edge, E ′ c , is above the statistical average of the potential (V = 0 for our choice of energy reference), the corrected mobility edge, E c , as calculated either by the method of Ref. [64] or by our self-consistent renormalized approach, is below the statistical average of the potential. In addition, we find that the renormalized self-consistent approach predicts values of E c in very good agreement (within 5 − 7%) with those of Ref. [64]. These results support our method to estimate E c .

Anisotropic disorder
We now apply our method to anisotropic disorder in the 3D single-speckle configuration. The mobility edge is found by searching the root of  Figure 16. (color online) On-shell (E ′ c ) and renormalized (Ec) mobility edges versus the disorder amplitude VR for the 3D (single-speckle) case with σ /σ ⊥ = 5.8.
the self-consistent equation (85). Note that the averaging of the angular dependence of Σ ′ in Eq. (85) is justified a posteriori by the weakk-angle variations of Σ ′ found around its mean value at E c (with standard deviations less than 10 − 15%). This is illustrated in Fig. 15, which presents the angular variations obtained numerically in the calculation of ∆(E c ), for typical values of V R and for an anisotropy of σ /σ ⊥ = 5.8.
The on-shell (E ′ c ) and renormalized (E c ) mobility edges are shown in Fig. 16. As for isotropic disorder, it is eye-catching that the shift of the energy states completely changes the behavior of the mobility edge. While the on-shell mobility edge, E ′ c , is above the statistical average of the potential, the renormalized mobility edge, E c , is below. This behaviour seems very robust for 3D speckle disorder. It was found for isotropic 3D speckles (see Ref. [64] and Sec. 7.4.2), as well as other models of speckle potentials with structured correlations [79].

Conclusions
Disordered potentials with finite-range correlations are often characterized by a counter-intuitive and interesting behaviour [32,28,29,25,30,79]. These are directly related to the microscopic statistical properties of the potential, hallmarked by the disorder correlation function. In this paper we have focused on anisotropy effects in 2D and 3D correlated disorder. We have quantitatively studied the transport and localization of matter waves by using mesoscopic transport theory [81] and a standard on-shell self-consistent perturbative approach [75]. The latter, first pioneered by Vollhardt and Wölfe [102,98], remains the most powerful, quantitative, microscopic approach to Anderson localization in dimension higher than one (d 2), in spite of the unavoidable problem of describing the physics inside the critical region in d > 2. Within this approach, we have characterized incoherent diffusion, quantum corrected diffusion and localization tensors versus the particle energy and found rich diffusion and localization properties. We have supported the general theory with application to speckle potentials in 2D and 3D.
In the 2D case, we have considered an anisotropic Gaussian correlation function as used in Refs. [53,54]. The energy-dependences of relevant quantities are studied: For E ≪ E σ ⊥ , in the white-noise limit, we find τ E,k ∝ 1 for the scattering time and D B ∝ E for the Boltzmann diffusion tensor, which are both isotropic. For E ≫ E σ ⊥ , we find τ E,k ∝ √ E and D B ∝ E 5/2 . As a general rule, the anisotropy of the disorder (ξ), of the scattering time (ξ s ) and of Boltzmann diffusion (ξ B ) are all different. The scattering time shows an inversion of anisotropy from ξ s > 1 (for ξ > 1) at low energy to ξ s = 1/ξ (< 1) at high energy. In contrast, the transport anisotropy is always ξ B > 1 (for ξ > 1) but shows a strongly nonmonotonic behaviour as a function of energy with a marked maximum at E ∼ E σ ⊥ . The anisotropy of localization is simply the square root of that of transport. For typical experimental parameters, we found that it is very small in observable regimes, except for very strongly anisotropic disorder. So far, experiments have only studied the classical regime [53,54] and our study offers scope for future studies of quantum transport and localization in 2D speckle potentials.
In the 3D case we have considered the strongly anisotropic correlation function of speckle potentials obtained with a single laser. Here, the energy dependence of relevant quantities are the following: For E ≪ E σ ⊥ , we find τ E,k ∝ 1 and is slightly anisotropic, while D B ∝ E and is significantly anisotropic, which is due to anisotropic suppression of the white-noise limit in the model we used. For E ≫ E σ ⊥ , we find τ E,k ∝ √ E and D B ∝ E 5/2 , both being anisotropic. We have also analyzed the anisotropy of transport as a function of the configuration anisotropy. We found that it is almost independent of the energy, and has a the behaviour ξ B = 0.59ξ + 0.21ξ 2 . In our approach, the anisotropy of the localization tensor is the square root of that of the Boltzmann diffusion tensor. We have also studied the behaviour of the 3D mobility edge. To do so, we have extended the on-shell approach and proposed a way to renormalize energies. We have found a striking agreement of our method with the more involved method based on SCBA developed in Ref [64] for isotropic disorder. The effect of renormalizing energies does not alter the overall energy dependence of the quantities discussed above, but may be important for direct comparison to energyresolved experimental measurements. As regards the mobility edge, we have found that the renormalization of energies has both a quantitative and qualitative impact. In particular, we find that, as for isotropic disorder, the renormalized mobility edge is below the average value of the disorder.
Finally, our results and method may provide a guide line to future experiments investigating the so-far unexplored effect of anisotropy in quantum transport of matter waves. In the case of ultracold atoms, to which our study directly applies, the transport properties can be probed by direct imaging of the atoms and control of the energy. First experimental studies of Anderson localization of 3D matter waves in anisotropic speckle potentials have been reported [57,58]. Our study is directly relevant to these experiments. For a detailed comparison of theoretical predictions and experimental observations, see Ref [79]. In addition, the effects discussed in this manuscript can be expected for other kinds of waves and/or other models of disorder, and are particularly relevant to new systems where the disorder correlations can be controlled [10,57,58,67,112,113,114].
(i) The eigenvalues λ n E and the eigenvectors φ n E,k of Eq. (A.1) are real.
Proof. By multiplying Eq. (A.1) by G † (E, k), we obtain . The latter is Hermitian since G † (E, k) * = G(E, k) and U E k,k ′ is real and symmetric. Therefore all the eigenvalues λ n E are real. By taking the complex conjugate of Eq. (A.3), dividing by G(E, k) and comparing it to Eq. (A.1), we obtain that the functions φ n E,k are real. If U E k,k ′ is positive-definite, the eigenvalues λ n E are positive. In particular, this is always true in the Born approximation p . When U E k,k ′ is symmetric and positive-definite, we can write it as where d k ′′ > 0 and Q is an orthogonal operator. For any vector of components It shows that M E k,k ′ is positive definite. Its eigenvalues λ n E are therefore positive. (ii) The eigenvectors φ n E,k can be chosen to satisfy the orthonormalization condition Proof. This is an immediate consequence of the fact that, according to Eq. (A.3), the functions G † (E, k) φ n E,k are eigenfunctions of the Hermitian operator M E k,k ′ . (iii) The eigenvectors φ n E,k satisfy the completeness relation Proof. This follows from the fact that the eigenfuntions G † (E, k) φ n E,k of the matrix M E k,k ′ , Eq. (A.3), form a complete basis.
(iv) The irreducible vertex function U E k,k ′ can be expressed as Proof. We multiply both terms of Eq. (A.1) by φ n E,k ′ and sum over n. Equation (A.6) is recovered by using the completeness relation Eq. (A.5).
p In this case, U E k,k ′ =C(k − k ′ ) is symmetric and positive-definite. This latter property is assured for any disordered potential by the fact that the power spectrumC(k), being the Fourier Transform of the autoconvolution product of the potential, is positive for any k. and the corresponding eigenvector is proportional to the inverse scattering mean free time: Proof. This is a direct consequence of the Ward identity [98]: Proof. This is a consequence of the parity of the vertex U E k,k ′ , in particular, In this case no diffusion would be observed. As noticed above, however, the conservation of particle number, through the Ward identity, imposes that there is one eigenvalue equal to one. As there is no other conserved quantity in the system we are considering, we can assume that the eigenvalue λ = 1 is not degenerated and that there is a finite gap between this eigenvalue and the rest of the spectrum when (q, ω) → 0 [115,116]. This suggests the following ansatz for the solution of the BSE (17)- (18) [see Eq. (20)], in the small (but non-zero) q and ω limit: where φ 1 k (q, ω, E) and 1 + λ(q, ω, E) are solutions of the eigenequation 14) The latter are the first eigenvalue and eigenvector at small (q, ω), and reduce to Eqs. (A.7) and (A.8) when (q, ω) = (0, 0), respectively. We then write f k (q, ω, E) = f E,k + F k (q, ω, E) the expansion of f k (q, ω, E). Making the ansatz φ 1 k (q, ω, E) = n a n (q, ω, E)φ n E,k , we find λ(q, ω, E) = n a n (q, ω, E) Finally, the coefficients a n (q, ω, E) are found by imposing that Eq. (A.13) solves the BSE. After some algebra one finds a 1 (q, ω, E) = 1 and a n (q, ω, E) =

Appendix A.4. On-shell approximation
We now proceed to the on-shell (weak disorder) approximation, and we neglect the effect of disorder on the spectral function. Equation (A.2) becomes f E,k ≈ τ E,k A 0 (k, E), (A. 16) where τ E,k is the on-shell scattering mean free time [τ E,k ≡ τ s (E, k Ek )], A 0 (k, E) = 2π δ[E−ǫ(k)] and ǫ(k) are, respectively, the disorder-free particle spectral function and dispersion relation. An explicit calculation of the small (q, ω) expansion of f k (q, ω, E), Then, making use of the parity properties of the functions φ n E,k [Eqs. (A.11) and (A.12)], τ E,k (even function ofk) and ∇ k ǫ(k) (odd function of k), we finally obtain φ 1 k (q, ω, E)f E,k = 2πγ k (q, E)/ τ −1 E,k where γ k is given by Eq. (56) and λ(q, ω, E) = 2N 0 (E) [i ω − q·D(E)·q] / τ −1 E,k with the diffusion tensor of Eq. (57). The solution of the BSE is thus given by Eq. (52) with Eqs. (53) and (55). Note that this expression for the diffusion constant is quite general (only the on-shell approximation has been made), provided that the full irreducible vertex function U is considered in the eigenequation (A.1). In Sec. 5.1 the Born and Boltzmann approximations are made U = U B [see Eq. (58)].
Let us consider for instance the first order term in ω. We find F k (q, ω, E) ≈ ω 2 [f E,k G † (E, k) − f E,k G(E, k)]. In the on-shell approximation this equation appears to go as the square of a δ-function, and one has to handle this divergence correctly [86]: we assume that f E,k G(E, k) ∼ 2π c δ(E − ǫ(k)), where the factor c is calculated by imposing that the integral over energy of f E,k G(E, k) remains invariant, i.e. c = dE 2π f E,k G(E, k). With this method, we find f E,k G(E, k) = i(τ 2 E,k / 2 )A 0 (k, E) and therefore F k (q, ω, E) ≈ ω i(τ 2 E,k / 2 )A 0 (k, E), as in Eq. (A.17). Following the same method, we can calculate the other terms in Eq. (A.17). Finally note that Eq. (A.17) also assumes that τs(E, k) is a smooth function of k, such that ∇ k τs(E, k) ≈ 0.
to the orthonormal cylindrical harmonics, with the prefactor determined by the normalization condition (A.4): In the calculation of the diffusion constant, it is actually possible to see that only the first term plus the l = 1 terms (with m = −1, +1) in the summation of the right-hand side of Eq. (57), contribute to the diffusion coefficient. More precisely the on-shell scattering mean free time τ E,k does not depend onk, υ x (respectively υ y ) is a 2πperiodic and even (resp. odd) function of θ, and Z +1 l (resp. Z +1 l ) is 2π/l-periodic and even (resp. odd). Therefore, when performing the angular averaging of the product τ E,k υ i φ n E,k in Eq. (57), one finds that only the term with l = 1 and m = +1 (resp. m = −1) couples to υ x (resp. υ y ) and contribute to D In the calculation of the diffusion constant, using the same type of symmetry arguments as in the 2D case, we find that only the l = 1 (with m = −1, 0, 1) terms couple to υ and contribute in the summation of Eq. (57). We thus find which agrees with the expression found in Ref. [62].

Appendix B.3. Three-dimensional isotropic speckle
A simple model of 3D speckle with isotropic correlation properties, is found when considering the light pattern obtained inside an integrating sphere lit by a laser beam of wavevector k L . The real-space correlation function is given in Eq. (86) and the associated power spectrumC is isotropic. Although this isotropic model is unrealistic from an experimental point of view, it is useful here in two respects. First, it bears the same divergence as the anisotropic 3D models of disorder considered in Sec. 3:C(k) ∝ 1/|k| when |k| → 0. Second, several properties of this model are analytical and known [61,62], and therefore provides a test for our numerical methods. As done previously, for the diagonalization of the integral operator (58) we use 2 7 × 2 7 points regularly spaced on the k-space shell |k| = k E . Some eigenfunctions and eigenvalues of Eq. (58) are presented in Fig. B1. We indeed find spherical harmonics [see Eq. (B.6)], and the eigenvalues λ n E agree well with theory [Eq. (B.5) withC given by Eq. (B.8), not shown on the figure]. We further incoporate these results in Eq. (57). Figure B2 presents the numerical results for the Boltzmann diffusion constant (red dots) which agree very well with the analytic formula (solid black line) found when incorporating Eq. (B.8) into Eq. (B.7). Note that we recover the same asymptotic behaviours as for our anisotropic cases: D B (E) ∝ E for E/E σ < 1/2 and D B (E) ∝ E 5/2 for E/E σ ≥ 1/2. In particular, those tests show that the discretization used here correctly treats the |k| → 0 divergence.

Appendix C. Conductivity
Appendix C.1. Einstein relation As presented in Sec. 2.5, we expect σ(ω = 0) ∝ D in the linear response regime.
Here we calculate σ B (ω = 0) in the Boltzmann approximation and verify this relation explicitly, which enables us to find the proportionality factor in Eq. (32). Let us first rewrite the Boltzmann diffusion tensor as where J k is the renormalized current vertex : We want to calculate the conductivity σ B in the ladder approximation. We have to evaluate