Quantum transport of atomic matter waves in anisotropic two-dimensional and three-dimensional 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 in 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 two-dimensional (2D) and three-dimensional (3D) speckle potentials. In particular, we discuss the anisotropy of single-scattering, diffusion and localization. We also calculate disorder-induced shift of the energy states and propose a method to include it, which amounts to renormalizing 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 relevant for current matter wave experiments, where the disorder is created with 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 widely universal behaviour [23], but observable features significantly depend on the details of the system. It has renewed interest in the context of ultracold matter waves [8][9][10][11][12]. On 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 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 matter waves [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 coldatom kick-rotor setups [50][51][52], study of classical diffusion in two-dimensional (2D) speckle potentials [53,54], coherent back-scattering [55,56] and evidence of Anderson localization in non-interacting Fermi [57] and Bose [58] gases in three-dimensional (3D) speckle potentials.
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 (section 2) and the models of disorder we focus on in 2D and 3D (section 3). We then present a detailed description of the theoretical framework pioneered in [75,80], which is intended to be pedagogical. We study single-scattering (section 4), Boltzmann diffusion (section 5) and localization (section 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 sections 4-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 (section 7). We also show that energydependent 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 compared 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 ultracold atoms in speckle potentials in the conclusion (section 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 [80,81]. Consider a wave of momentum k and velocity υ = ∇ k /h ( (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 random impurities. Three typical energy-dependent length scales can be identified, which characterize three basic effects induced by the disorder (see figure 1).
First, single scattering from impurities depletes the k-wave states, which can be seen as quasiparticles in the disordered medium, with a finite lifetime τ 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 of its initial state, and primarily the memory of its initial phase. Then, multiple scattering defines the second length scale, namely the transport (Boltzmann) mean-free 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, r 2 ∼ 2d D B t with D B = υl B /d the Boltzmann diffusion constant (d is the space dimension) [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 followed in either direction, 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]. 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 multicolour line) in a 2D disordered landscape (blue surface). Along its trajectory, the wave loses the memory of its phase (encoded in the various colours along the trajectory) on the characteristic length l s (scattering mean-free path). Multiple scattering then deflects the trajectory and the wave loses the memory of its direction on the characteristic length l B (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 L loc (localization length).
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 the localization relies on two characteristics of the medium: coherence along the multiple-scattering paths and return probability 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 [27,28,83,84]. 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, subjected to some static randomness. Its dynamics is governed by the single-particle Hamiltonian H = H 0 + V (r), where H 0 is the disorder-free, translation-invariant, 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 4 , V = 0. The evolution of the wave function between t 0 and t > t 0 is determined by the retarded single-particle propagator , where the Heaviside step function (t − t 0 ) accounts for temporal ordering. In the energy domain 5 , G is the retarded Green operator where E is the particle energy. It is the solution of the equation where G 0 = (E − H 0 + i0 + ) −1 is the disorder-free retarded Green function associated with the unperturbed Hamiltonian H 0 .

Properties of the disordered medium
In a disordered medium, meaningful observable quantities correspond to statistical averages over realizations of the disorder. When averaging over disorder realizations, some quantities can be written in terms of the average Green function G(E), for instance the spectral function (see below). The Born series of equation (2), averaged over the disorder, reads 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 developed in powers of V thanks to equation (3) so as to determine . The average Green function then reads If the disorder is homogeneous, i.e. if its statistical properties are translation invariant [87], then the disorder-averaged Green function is diagonal in k-space 6 : A(E,k) k k E k 0 Schematic representation of the spectral function A(E, k) of a particle of energy E =h 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 In the presence of disorder the spectral function is shifted and broadened (black line).
where (k) is the dispersion relation associated with 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] 2π It contains all the information about the spectrum of the disordered medium. Using equation (1), it yields 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, equations (7) and (9) yield with and the real and imaginary parts of , respectively. As represented schematically in figure 2, for a particle in free space ( (k) =h 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 centred at 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 lifetime 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 knowledge of the real and imaginary parts of the self-energy (see equation (10)), or, according to equation (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 sections 5-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 figure 2). In section 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. That 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 equation (13) in a form indicating explicitly the initial conditions, using the relation ρ(t) = (t − t 0 ) e −iH (t−t 0 )/h × ρ(t 0 ) e +iH (t−t 0 )/h . When averaging over the disorder, if there is no correlation 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) e iq·R e −iωt k,k (q, ω, E) and with k ± ≡ k ± q/2, k ± ≡ k ± q /2, E ± ≡ E ±hω/2 and (q, ω) the Fourier conjugates of the space and time variables 7 . As discussed above, disorder averaging features a translational invariance in space so that equation (15) depends only on the difference R = r − r . For the same reason, translational invariance, or equivalently momentum conservation, imposes that the sums of the in-going wavevectors (k + and k − ) on the 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 equation (16).
As can be seen in equations (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 of equation (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.
Analogous to equation (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 [89] 8 : 7 We use the Fourier transformf (q, ω) ≡ dr dt f (r, t) exp[−i(q · r − ωt)]. 8 In this context, the inverse of an operator is defined by More explicitly, the (k, k ) component of a four-point vertex which fulfils momentum conservation is k,k (q, ω, E), such that k + , k − | |k + , k − ≡ (2π) d δ(q − q ) k,k (q, ω, E), and and Therefore, equation (20) reads 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 [75,90] to solve the BSE. It then gives access, via equation (23) to , which is the quantity of interest (see equations (13)- (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 position-energy 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 equations (14), (15) and (25) as the space-time Fourier transform of the diffusion pole 1/[ihω −hq · D(ω, E) · q]. In the long-time limit, we will encounter two different situations. Firstly, 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 Secondly, 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 function P(R) decays exponentially 9 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 with 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 10 [6] where υ i =hk i /m is the velocity along axis i. As the structure of equation (30) is reminiscent of the definition of the four-point vertex (see equation (16)), calculations of the conductivity tensor can also be represented diagrammatically. The skeleton diagram, shown in equation (31), consists of the in-and out-going velocities υ and υ and of a bubble made of a retarded (top line) and an advanced (bottom line) Green function. As for , the scattering events of the top and bottom lines can be correlated (see e.g. equations (18) and (19)).
Thanks to Einstein's classical argument, it was realized that, at thermal equilibrium, in a gas subjected 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. [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 C.1), which in our system yields 9 Note that K 0 (x) ∼ e −x √ π/2x for x 1. 10 This corresponds to the more general definition σ i, j (ω, E) = ∞ 0 dt e iωt Tr{δ(E − H ) j i (x, t) j j (x)} ( j is the current operator) where the correlations between G and G † have been dropped (see e.g. [6]). Schematic view of the apparatus used to create an optical speckle pattern. A laser beam is diffracted by a ground-glass plate diffuser of pupil function I D (ρ), 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.
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 figure 3 and [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 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 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 11 . For the configuration of figure 3, in the paraxial approximation, we find with 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 (two dimensions)
If the atoms are confined in a 2D geometry by a strong trapping potential along z centred on z = 0, they experience a disordered potential with correlation function , with σ ⊥ = σ ⊥x and ξ = σ ⊥x /σ ⊥y the configuration anisotropy factor. The Fourier transform gives the power spectrum Without loss of generality, we assume that ξ 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 [53] where a quasi-2D Bose gas of width l z is subjected to a speckle created by an isotropic Gaussian laser beam shone at an angle θ with respect to the plane of atoms, if l z σ ⊥ σ . In this case ξ 1/sin θ (θ π/6 for the experiment of [53]).

Single speckle (three dimensions)
In the 3D case, the disorder correlation function C(r) is given by equation (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, 4 f > w and C(r) is elongated along z . The corresponding disorder power spectrum reads where k ⊥ is the projection of k on 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 the absence of white-noise limit, which reflects the long-range correlations of the potential 12 . 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 section 2.1: the scattering mean-free time.

Scattering mean-free time
In order to calculate the scattering mean-free time (11), we retain only the lowest order contribution to the self-energy (Born approximation). Within this approximation, the Born series of equations (3)-(4) is truncated after the first two terms, which, according to equation (6), yields For homogeneous disorder, k| whereC(k) is the disorder power spectrum. Using equation (11) and the disorder-free Green function, we thus have 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 section 7.1).
In the case of isotropic disorder (i.e.C(k − k ) =C(|k − k |)), 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 that the scattering depends on the relative direction of k and k, but that it also depends on the direction of the incoming wave k.

Anisotropic Gaussian speckle (two dimensions)
Let us consider the 2D anisotropic speckle potential of geometrical anisotropy factor ξ introduced in section 3.1. ReplacingC(k) by equations (35) in equation (40) and using the disorder-free dispersion relation of the vacuum in equation (41), we obtain the scattering meanfree time wherek ≡ k/|k| is the unit vector pointing in the direction of k, k is the k-space solid angle, k E ≡ √ 2m E/h is the momentum associated with energy E in free space and E σ ⊥ ≡h 2 /mσ 2 ⊥ is the correlation energy of the disorder. The scattering time (42) is plotted in figure 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 τ E,k ≡ τ s (E, k Ek ).  (42) for |k| = k E ) along thek x (solid red line) andk y directions (dotted blue line) for the 2D speckle potential defined in section 3.1 with ξ = 4. The solid black lines are the isotropic low-energy limits obtained for k E σ ⊥ 1 (equation (43)) and the high-energy limit obtained for k E σ ⊥ ξ (equation (44)). The insets show the angular dependence of τ E,k at two different energies (with the parametrization k ≡ (cos θ, sin θ )). The points on the lines are colour-and shape-coded to match those in the insets.
In the low-energy limit, k E σ ⊥ 1, we have which is displayed in figure 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) ⊥ ξ (see section 3.1) and τ E,k is isotropic, constant, and it only depends on the product V 2 R σ ⊥x σ ⊥y (up to corrections of relative order 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,k x < τ E,k y . More precisely, we find which is shown in figure 4 (right-hand side black lines). In particular, we find that in the high- Anisotropy factor of the scattering time, ξ s = τ E,k x /τ E,k y , as a function of E/E σ ⊥ and ξ , for the 2D speckle potential of section 3.1. The red lines are the low-(ξ s → 1) and high-energy limits (ξ s → 1 ξ ) (see equations (43) and (47)).
It is also interesting to study the anisotropy factor of the scattering time which is shown in figure 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 figure 5). With increasing 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 equation (43), an explicit calculation yields For k E σ ⊥ ξ , using equation (44), we obtain which shows that the anisotropy factor of scattering is proportional to the inverse of the geometrical anisotropy (right-hand side red line in figure 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 (in contrast to the lowenergy limit case). Therefore, for any value of ξ , τ E,k exhibits 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 [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 figure 6. In each directionk the spectral function peaks at 4τ E,k /h 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 section 7.2).

Solution of the Bethe-Salpeter equation
In the independent scattering (Boltzmann) and weak disorder (Born) approximations, only the first term in equation (19) is retained and the irreducible vertex function U equals the disorder structure factor [81]: or equivalently Then, incorporating equations (48)- (49) into the BSE (17)- (18) and expanding it in series of U , one finds where the diffusion reduces to ladder diagrams: It describes an infinite series of independent scattering events, which leads to Drude-like diffusion.
In appendix A, explicit calculations are detailed. In brief, in the long-time (ω → 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 f E,k ≡ f k (q = 0, ω = 0, E) (see equation (22)) and φ n E,k (λ n E ) are the eigenvectors (eigenvalues) of an integral operator involving the disorder correlation function and 13 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 equation (54) is equal to one, λ n=1 E = 1 (see 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 (55) 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 is the on-shell scattering meanfree time (see equation (40)), and . . . k|E represents integration over the k-space shell defined by (k) = E (see equation (41)). The functions φ n E,k and the real-valued positive numbers λ n E are the solutions of the integral eigenproblem (54), which becomes, in the on-shell approximation [75]. It follows from equation (57) that the incoherent (Boltzmann) diffusion tensor D B (E) is obtained from the two-point disorder correlation function C(r), which determines τ E,k (see equation (40)) as well as φ n E,k and λ n E (see equation (58)). In the isotropic case (for details see appendix B), equation (58) is solved by the cylindrical, Z ±1 l (2D; see appendix B) or spherical, Y m l (3D) harmonics, the same level harmonics (i.e. with the same l) being degenerate in λ n E . Then, it follows from the symmetries of the cylindrical/spherical harmonics that only the first term in equation (57) plus the p-level harmonics (Z ±1 1 in 2D and Y m 1 with m = −1, 0, 1 in 3D) couple to υ and contribute to D B (E). Incorporating the explicit formulae for φ n E,k and λ n E (see equations (B.1)-(B.7)), we then recover well-known expressions for isotropic disorder [61][62][63]65].
In the anisotropic case, harmonics couple, and the φ n E,k are no longer cylindrical/spherical harmonics.

Anisotropic Gaussian speckle (two dimensions)
Consider again the 2D anisotropic speckle potential of section 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 equation (58). We solve equation (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): The eigenvalues and some eigenfunctions obtained numerically are shown in figure 7 for various values of E/E σ ⊥ . As discussed above, we find λ n=1 only the first term in the right-hand side of equation (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 with the orbitals with n > 1 in equation (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 contribute (see below). In the low-energy limit, one can develop equation (35) in powers of |k|. Up to order O(E 2 /ξ 4 E 2 σ ⊥ ), the first three eigenfunctions are given by 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 in figure 7). In the very low-energy limit, the disorder power spectrum becomes isotropic and constant,C(k) V 2 R πσ 2 ⊥ /ξ (see section 3.1 and equation (35)). The orbitals φ n E,k are thus proportional to the cylindrical harmonics, which are exact solutions of equation (58) in the isotropic case (see appendix B, and use the parametrizationk x = cos θ andk y = sin θ ). In contrast to the isotropic case where the values of λ n E are degenerate 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 equations (60) and (61)). When the energy further increases, the anisotropy plays a more important role and the harmonics are more and more distorted (see figure 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.  (62) and (63)), with the isotropic white-noise limit D x Incorporating the values of λ n E , φ n E,k and τ E,k in equation (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 equations (43), (60) and (61), we find that the first term in the right-hand side of equation (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 x The scaling D u B (E) ∝ E 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 Z −1 and which are displayed in figure 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 figure 7). Hence, for the same reasons as in the isotropic case (see appendix B), only the Z ±1 1 -like orbitals couple to υ in equation (57) and contribute to D B , while the others do not. The associated λ n E increase (see figure 7), the weight of the second term in equation (57) increases, and the components of the diffusion tensor show very different behaviour in the large-E limit. For k E σ ⊥ ξ , we found τ E,k ∝ k E (see section 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 equation (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 figure 8). This scaling is similar to that found for isotropic disorder [62]. As shown in figure 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 non-monotonic behaviour versus E, with a marked peak (see inset of figure 8).
The Boltzmann transport anisotropy factor ξ B is shown in figure 9 for various configuration anisotropies ξ . As well known, the scattering and transport mean-free times are different quantities in correlated disorder, due to angle-dependent scattering [81,99,100]. In particular, in the 2D speckle we consider, we do not find any inversion of the anisotropy of the diffusion, in contrast 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 figure 8 for a cut at ξ = 4), which increases with the geometrical anisotropy ξ (see inset of figure 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 [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 (section 6.1), then the self-consistent theory (section 6.2), and apply it to the 2D speckle potential (section 6.3). The 3D case, which follows the same route, is discussed in section 7.3.

Weak 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 conjugate. This is diagrammatically translated into crossing correlation lines as in the second term of equation (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 [6,75,90,101]: the cooperon (equation (64)) and the first two Hikami boxes (equations (65) and (66)).
where the cooperon X is the sum of maximally crossed diagrams and is the renormalized vertex function (see appendix C.2).
Using time-reversal invariance [23,81,98,102], the cooperon X can be expressed in terms of the diffusion (defined in equation (51)) The diffusion pole carried by in the limit (ω, q) → 0 leads to a divergence of X when ω, k + k → 0. In 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 14 . In particular, it follows from equation (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 [98,102] proposed to selfconsistently replace D B (E) by the dynamic diffusion tensor D * (ω, E) in the right-hand side of equation (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 [80,98]. Generalizing this standard approach to anisotropic disorder yields In dimension d 2, the integral in the right-hand side of equation ( (75)) and the dotted ones for high values of E/E σ ⊥ (equation (76)). The dashed grey lines indicate typical values of the imaging resolution (L res ) and the system size (L sys ) in ultracold-atom experiments, see text at the end of section 6.3.
domain 15 . 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 (two dimensions)
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 equation (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 section 2.4, it leads to the exponentially decreasing propagation kernel (28). Solving equation (73) then yields the anisotropic localization tensor where l av 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 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 section 5.2 into equation (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 15 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 section 7.4.2. 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 equations (62) and (63), we find L x,y where the upper sign holds for direction x, and the lower sign for direction y. Equation (75) corresponds to the solid black lines in figure 10. As D B is almost isotropic for E/E σ ⊥ 1 (see figure 8), L loc is also almost isotropic in the whole range presented in figure 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 figure 8. According to equations (62) and (63) (retaining only the lowest-energy term), this regime appears for When ξ = 4 (as in figure 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 explicitly seen in equation (76), the 2D localization length increases exponentially at large energy (hence the limited energy range in figure 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 in figure 10 by the grey dashed lines corresponding to L res 15 µm and L sys 2 mm, which are typical values extracted from [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 figure 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 sections 4-6 to the 3D speckle potential of section 3.2. We discuss single scattering (section 7.1), Boltzmann diffusion (section 7.2) and localization (section 7.3) properties and the position of the mobility edge (section 7.4).

Single scattering
Inserting equations (36) and (37) into equation (40), we find the scattering mean-free time  Figure 11. Scattering mean-free time τ E,k in the 3D case (equation (77)) with σ /σ ⊥ = 5.8 with |k| = k E , in the (k x ,k y ) plane (solid red line) and along thek z direction (dotted blue line). The black lines are the low-energy (k E σ ⊥ 1, see equation (79)) and the high-energy (k E σ ⊥ 1, see equation (80)) limits. Note that in both limits τ E,k is anisotropic, although for k E σ ⊥ 1, the anisotropy is very small, ξ s 1.002. The insets show the angular dependence of τ E,k at different energies (with θ = (k;k z )). The points on the lines are colour-and shape coded to match those in the insets.
which is shown in figure 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 azimuthal 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,k z < τ E,k ⊥ ) for all values of E. As for the 2D case analysed in section 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 16 . This can be attributed to the strong anisotropic divergence ofC(k) when |k| → 0, which reflects the long-range correlations of the disorder (see section 3.2). More precisely, for |k|σ ⊥ 1, we havẽ Replacing this expression into equation (77), we then find which is independent of E. Equation (79) is plotted as solid black lines on the left-hand side of figure 11. Note that τ E,k does not become strictly isotropic in this limit. However, the residual anisotropy of the scattering time, found from equation (79) and from the anisotropy ofc(k) in equation (78), is very small, and practically unobservable (τ E,k ⊥ /τ E,k z 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 figure 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 equation (77), which is done on a sphere of radius k E containing the origin, can be reduced to integratingc 1sp on the plane which is tangent to the sphere at the origin. We then find In particular, we find πσ (both shown as right-hand side solid black lines in figure 11). The anisotropy of the scattering then becomes significant for the parameters of figure 11, τ E,k ⊥ /τ E,k z = √ πσ /2σ ⊥ in this limit. The highenergy 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 equation (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 analysed previously, by solving equation (58) numerically and incorporating the results into equation (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 equation (58) for different energies, as well as the topography of the eigenvectors of equation (58) figure 12 for σ /σ ⊥ = 5.8. Similarly as for the 2D case, we found 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 degenerate in a given l-like level. More precisely, due to the cylindrical symmetry of the power spectrum, the value of λ n E associated with the Y +m l -like and Y −m l -like orbitals is 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 cylindrical invariance of the correlation functionC(k) around the axisk z . For the same symmetry reasons as in the isotropic case (see 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 equation (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 B as a function of the configuration anisotropy ξ = σ /σ ⊥ , at E/E σ ⊥ = 6 × 10 −3 , 6 × 10 −1 and 60. The dotted line is a fit of all the data, which gives ξ B = 0.59ξ + 0.21ξ 2 .
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 figure 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 3 -like orbitals calculated numerically are displaced compared to the associated spherical harmonics. Therefore, their contribution in equation (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 figure 13(a)), and the anisotropy of D B is thus reversed with respect to that of τ E,k (we recall that we found τ E,k z < τ E,k ⊥ for any E, see section 7.1). This is due to the fact that the (Y 0 1 -like) orbitals contributing to D z B are associated with values of λ n E larger than those contributing to D x,y B (in figure 12, the φ n E,k are numbered by decreasing eigenvalues). Secondly,C(k) shows a strong anisotropic, infrared divergence in the paraxial approximation (see sections 3.2 and 7.1). Following up with the scaling ofc 1sp (k), equation (78), used to show that τˆk ,E is independent of energy for k E σ ⊥ 1, and inserting it into equation (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 equation (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 on the left-hand side of figure 13(a) for k E σ ⊥ 1 (i.e. E E σ ⊥ ). This is another manifestation of the absence of a white-noise limit 17 .
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 obtain φ n E,k ∝ 1/k E and D u B (E) ∝ E 5/2 . This scaling is confirmed in figure 13(a) by fits to the data for E E σ ⊥ (right-hand side dotted lines). Remarkably, in spite of the different contributing terms in equation (57) at low and high values of E, the transport anisotropy is nearly independent of 10 (see inset of figure 13(a)). We have repeated the same study for various values of the configuration anisotropy, ξ = σ /σ ⊥ . They all show similar behaviour as a function of energy as reported in figure 13(a) for ξ = 5.8. In particular, we found the same scalings with energy and a diffusion anisotropy x,y B that is nearly independent of energy. In figure 13(b), we plot ξ B versus ξ for three values of the energy. We find that the diffusion anisotropy monotonically 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 expects ξ 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 figure 13(b).

Localization
In order to analyse 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, 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 . 17 A 3D white-noise limit would lead to the scaling D u B (E) ∝ √ E and an isotropic limit at low energy.
Explicitly, we have where L av loc = det{L loc (E)} 1/3 is the unique solution of L av loc 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 As already mentioned in section 6.1 the behaviour of L loc and D * is completely determined by that of D B in our approach. The anisotropies of L loc (E) are the square roots of those of D B (E) (see equation (81)) and the anisotropies of D * (E) are the same as those of D B (E) (see equation (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 in figure 13(b). At low energy, using the scaling of D u 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 section 7.2).

About the three-dimensional mobility edge
The self-consistent approach used above is expected to fairly describe the quantum transport properties [62,75,80]. It gives some quantitative estimates consistent with numerical calculations [105] and experimental data [56,106]. 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 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 nonuniversal quantity and should be determined from microscopic theory. In this respect, the onshell approximation is questionable 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 .

Energy renormalization.
In order to improve the self-consistent method, one could in principle use the more sophisticated approach of [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 [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 resources seem out of present-day possibilities. In order to overcome this issue, we have proposed in [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 equation (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 replacing 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 energy-resolved 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 speckle disorder obtained inside an integrating sphere lit with a laser beam, the real-space correlation function of which reads [62,64] with σ the correlation length. The associated power spectrum (see 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 section 7.3 (see also [62]), the renormalized mobility edge E c calculated by our method (see section 7.4.1), and the mobility edge found using the self-consistent Born approximation (SCBA) in [64]. As clearly seen in figure 14, the disorderinduced 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 [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 [64]. These results support our method to estimate E c . Comparison of the mobility edges as calculated with the SCBA method (the full black squares are the results obtained by Yedjour and van Tiggelen in [64], which we reproduce here), with the on-shell method (E c , red crosses) and with the renormalized self-consistent approach (corrected E c , thick blue circles), for an isotropic 3D speckle potential. When comparing to figure 8 of [64], note that in [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 .

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 the selfconsistent equation (85). Note that the averaging of the angular dependence of in equation (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 figure 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 figure 16. As for isotropic disorder, it is eye catching that the shift of the energy states completely changes the behaviour 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 [64] and section 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 counter-intuitive and interesting behaviour [25, 28-30, 32, 79]. This is 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 [98,102], 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 [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 strongly non-monotonic 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 is as follows: 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 analysed the anisotropy of transport as a function of the configuration anisotropy. We found that it is almost independent of the energy, and has 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 [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 energy-resolved experimental measurements. As regards the mobility edge, we have found that the renormalization of energies has both 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 guideline 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 [79]. In addition, the effects discussed in this paper 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,[111][112][113].

A.1. Preliminary remark
First, let us notice that we have where A(E, k) is the spectral function defined in equation (9) and τ s (E, k) is the scattering mean-free time defined in equation (11).

A.2. Properties of equation (A.1)
The main properties of equation (A.1) and of its eigenfunctions are listed below: (i) The eigenvalues λ n E and the eigenvectors φ n E,k of equation (A.1) are real.
Proof. By multiplying equation (A.1) by G † (E, k), we obtain dk where M E k,k ≡ G † (E, k) U E k,k G(E, k ). 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 equation (A.3), dividing by G(E, k) and comparing it to equation (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 18 . 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 x k , we have  18 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.
Proof. This is an immediate consequence of the fact that, according to equation (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 eigenfunctions G † (E, k) φ n E,k of the matrix M E k,k , equation (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 equation (A.1) by φ n E,k and sum over n. Equation (A.6) is recovered by using the completeness relation (A.5). and the corresponding eigenvector is proportional to the inverse scattering mean-free time: Proof. This is a direct consequence of the Ward identity [98]  (vi) The eigenfunctions φ n E,k have the parity properties: Proof. This is a consequence of the parity of the vertex U E k,k , in particular,

A.3. Solution of the Bethe-Salpeter equation
Note first that, if equation (A.1) could be diagonalized with all eigenvalues different from one (λ n E = 1 for all n), it is straightforward to show, using equation (A.5), that we would 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 degenerate and that there is a finite gap between this eigenvalue and the rest of the spectrum when (q, ω) → 0 [114,115]. This suggests the following ansatz for the solution of the BSE (17)-(18) (see equation (20)), in the small (but non-zero) q and ω limit: where φ 1 k (q, ω, E) and 1 + λ(q, ω, E) are solutions of the eigenequation dk 14) The latter are the first eigenvalue and eigenvector at small (q, ω), which reduce to equations (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) a 1 (q, ω, E) Finally, the coefficients a n (q, ω, E) are found by imposing that equation (A.13) solves the BSE. After some algebra one finds a 1 (q, ω, E) = 1 and a n (q, ω, E) = where τ E,k is the on-shell scattering mean-free time (τ E,k ≡ τ s (E, k Ek )), 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), gives 19 Then, making use of the parity properties of the functions φ n E,k (equations (A.11) and (A.12)), τ E,k (even function ofk) and ∇ k (k) (odd function of k), we finally obtain where γ k is given by equation (56) and λ(q, ω, E) = 2N 0 (E) ihω −hq·D(E)·q /h τ −1 E,k with the diffusion tensor of equation (57). The solution of the BSE is thus given by equation (52) with equations (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 section 5.1 the Born and Boltzmann approximations are made, U = U B (see equation (58)).

B.1. Two-dimensional case
In the 2D isotropic case, inserting the cylindrical harmonics Z 0 = 1, Z +1 l = cos(lθ) and Z −1 l = sin(lθ) into equation (58), we find λ l,m E = 2π 0 dθ p(k E , θ ) cos(lθ) 2π 0 dθ p(k E , θ) , (B.1) 19 The small (q, ω) expansion of f k (q, ω, E) requires special attention in the on-shell approximation. Let us consider for instance the first-order term in ω. We find F k (q, ω, E) ≈¯h ω 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 /h 2 )A 0 (k, E) and therefore F k (q, ω, E) ≈hω i(τ 2 E,k /h 2 )A 0 (k, E), as in equation (A.17). Following the same method, we can calculate the other terms in equation (A.17). Finally note that equation (A.17) also assumes that τ s (E, k) is a smooth function of k, such that ∇ k τ s (E, k) ≈ 0.
where l 0 and m ∈ {−1, +1} are integer. In particular, we find λ l=0 E = 1 in agreement with equation (A.7). They are doubly degenerate for l > 0 and the corresponding normalized eigenfunctions are proportional 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 equation (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 equation (57), one finds that only the term with l = 1 and m = +1 (resp. m = −1) couples to υ x (resp. υ y ) and contributes to D x B (resp. D y B ). Then, inserting equations (B.1)-(B.3) into equation (57), we find .

(B.4)
This formula agrees with the result of [62], obtained by a different approach.

B.2. Three-dimensional case
In the 3D isotropic case, proceeding in a similar way, we find that the eigenvalues of equation (58)  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 equation (57). We thus find which agrees with the expression found in [62]. (second row) and D z B (third row) (with the parametrizationk = (k x ,k y ,k z ) ≡ (sin θ cos φ, sin θ sin φ, cos θ)), the red lines locate the nodal lines. From left to right E = 6.3 × 10 −3 E σ , E = 6.3 × 10 −1 E σ and E = 63E σ .

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 equation (86) and the associated power spectrum is isotropic. Although this isotropic model is unrealistic from an experimental point of view, it is useful here in two respects. Firstly, it bears the same divergence as the anisotropic 3D models of disorder considered in section 3:C(k) ∝ 1/|k| when |k| → 0. Secondly, several properties of this model are analytical and known [61,62], and therefore provide 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 equation (58) are presented in figure B.1. We indeed find spherical harmonics (see equation (B.6)), and the eigenvalues λ n E agree well with theory (equation (B.5) withC given by equation (B.8), not shown in the figure). We further incorporate these results in equation (57). . Note that we recover the same asymptotic behaviour 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.

C.1. Einstein relation
As presented in section 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 equation (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 where is defined in equation (51). It reads In the same way as before, and using the on-shell approximation formulae G(E, k) f E,k ∼ −i(τ E,k /h) 2 A 0 (E, k) and G † (E, k) f E,k ∼ i(τ E,k /h) 2 A 0 (E, k), we obtain σ (H 1 ) σ (H 2 ) and (C.11) C.3.3. Corrected conductivity tensor. We now consider the quantity J k − dk (2π ) d U Bk,k f E,k Jˆk . Using the relation U Bk,k = λ n E =1 λ n E φ n E,k φ n E,k and the parities of the functions φ n E,k (see equations (A.11) and (A.12)), one can show that (C.12) Therefore the Hikami contributions renormalize one of the J k /h back to the bare vertex υ, and we have (C. 13) which gives the final expression (70).