Position-dependent radiative transfer as a tool for studying Anderson localization: Delay time, time-reversal and coherent backscattering

Previous work has established that the localized regime of wave transport in open media is characterized by a position-dependent diffusion coefficient. In this work we study how the concept of position-dependent diffusion affects the delay time, the transverse confinement, the coherent backscattering, and the time reversal of waves. Definitions of energy transport velocity of localized waves are proposed. We start with a phenomenological model of radiative transfer and then present a novel perturbational approach based on the self-consistent theory of localization. The latter allows us to obtain results relevant for realistic experiments in disordered quasi-1D wave guides and 3D slabs.


I. INTRODUCTION
Among many other features, Anderson localization of waves is characterized by a halt of diffuse transport [1]. Since the formulation of the scaling theory of localization we understand that diffusion cannot entirely vanish in open media due to leakage of waves across the sample boundaries [2]. The result is a suppressed but scale dependent conductance, depending in a universal way on the size and dimensionality of the random medium. Later work established that many features of scaling theory can be understood from the self-consistent theory of localization [3], which adopts the constructive interference of time-reversed waves as the sole mechanism of the suppression of diffusion. More recently, this theory was extended to predict a position-dependent diffusion constant [4][5][6]. Microscopic derivations were provided using diagrammatic [7] and super-symmetric [8] approaches. Diffusion is suppressed deep inside the sample, yet hardly near the boundaries. This result is physically plausible, and consistent with scaling theory for macroscopic transport quantities such as, e.g., the conductance. It was tested against numerical simulations [9] and observed in an experiment [10].
If stationary transport is described by a spatially varying local diffusion constant, what does this imply for the dynamics, and in particular for the energy transport velocity v E ? In a weakly disordered three-dimensional (3D) medium the wave transport is diffusive with the diffusion constant given by D = v E /3 [11]. It was shown that v E is intrinsically a dynamic property, whereas the transport mean free path emerges by itself in stationary (DC) measurements, like in diffuse transmission through a slab of thickness L, T ∼ /L. This is important because the transport velocity can be very small due, for example, to strongly resonant scattering, and can thus lead to "small" diffusion constants. It can be then difficult to distinguish between situations in which D is small due to Anderson localization effects leading to small or due to small v E . Up to now, transport theories of local-ization essentially concentrated on and not on v E . So far we know that v E depends on many sample properties such as, e.g., the scattering crosssection of scatterers and their number density, but not on sample size or boundary conditions. Is v E well defined in the localized regime? What does the position-dependent imply for the transverse spreading of a wave packet in experiments similar to those of Ref. [12]? Similar questions arise for the Wigner delay time in reflection and transmission, and for coherent backscattering (CBS). In this work, we first develop general arguments for the delay time in a medium with position-dependent , valid under very broad conditions, and then present a perturbational approach to Anderson localization in the framework of the self-consistent theory of localization. The latter is applied to study wave dynamics, the transverse spreading of a wave packet, CBS, and time-reversal of localized waves.

II. FRIEDEL IDENTITY IN RADIATIVE TRANSFER
Depth-dependent extinction is very common in radiative transfer. Here we conjecture that the phenomenological equation of radiative transfer (EQRT) applies in the localized regime though with a depth-dependent scattering mean free path (z). This is clearly an oversimplified picture. In particular it disregards off-shell scattering (ω = kc) that becomes significant when the scattering mean free path becomes small. Nevertheless, it is consistent with the macroscopic picture of depth-dependent diffusion. In this section we make the connection between phase delay time, a highly mesoscopic wave property, and the EQRT, apparently a classical equation where all wave properties seem to have disappeared. The essential elements in this approach are the transport velocity and the link between stored energy and delay time established in condensed matter physics. We will reproduce some exact results from standard radiative transfer theory, such as the relation between incoming flux, source and energy density away from the boundaries, that will be necessary arXiv:1609.01556v2 [cond-mat.dis-nn] 13 Jan 2017 to elucidate the exact role of the incident flux for the delay time, and its scaling with the sample size.
We consider a slab of disordered medium confined between the planes z = 0 and z = L and made of isotropic, conservative scatterers. The incoming specific intensity is I(0, 0 < µ < 1, s) on the left and I(L, −1 < µ < 0, s) on the right. The EQRT can be written as [13][14][15] s v E I(z, µ, s) + µ∂ z I(z, µ, s) + 1 (z) I(z, µ, s) where s is the Laplace conjugate of time and µ = cos θ. All observables in this paper are ensemble-averaged and if no confusion exists, no explicit reference to this will be given. We have assumed the existence of an energy velocity v E independent of the direction of scattering µ and position (or depth) z. We can introduce the optical depth dτ = dz/ (z) and write with the energy density w(τ, s) = v −1 E 1 −1 dµI(τ, µ, s) (and equal to 2J/v E in terms of the source function J featuring in radiative transfer).
Let us first obtain a useful result relevant to the study of the delay time. Upon integrating Eq. (2) over the depth z and over all angles we see that The term in square brackets on the left-hand side contains transmitted and reflected fluxes F + (L, s) and F − (0, s), respectively, the right-hand side is the incident flux on both sides of the slab. The total energy is W = S B 0 dτ (τ )w(τ ) = S L 0 dzw(z), with the total optical thickness B = τ (L) and slab surface S. To discuss the delay time, we assume the incident flux independent of s (i.e., perfect delta functions in time): F − (L, s) + F + (0, s) = F in . The average transmission and reflection coefficients are then T (s) = F + (L, s)/F in and R(s) = F − (0, s)/F in , respectively. For s = 0 (time-integrated signal), we infer flux conservation R(0) + T (0) = 1. Taking the derivative of Eq. (3) with respect to s we obtain where the angular brackets · · · indicate that the ensemble averaging is to be carried out for the products T dφ T /dω and Rdφ R /dω. This equation makes the desired connection between stored energy and total (channel-summed) phase delay time. The second equality follows from the notion that the complex transmission coefficient is t = √ T exp(iφ T ), with φ T (ω) being the phase shift, and that 1 2 ln T + iφ T is a function of 1 2 s + iω, so that the Cauchy-Riemann equations give dφ T (ω)/dω = −d ln T (s)/ds at s = 0 (and similarly for the reflection coefficient). The relation between stored energy (or charge) and phase delay time is well known in different contexts: as Friedel's identity [16] in the context of screening of charge around impurities, or as Jauch formula relating phase delay time to the local density of states in scattering theory [17]; see also Ref. [18] for a related discussion and a list of relevant references. In this case the relation between stored energy in radiative transfer and the mesoscopic density of states is controlled by the transport velocity v E , which, for the moment, appears just phenomenologically in the EQRT. The phase delay time, in turn, can easily be related to the first moment of the scattered intensity, i.e. ∞ 0 dtT (t)t = T (ω)dφ T (ω)/dω , and similarly for reflection, which explains its interpretation as average delay time.
Let us next consider the stationary energy flow putting s = 0 and agree to have unit incident flux. Upon integrating Eq. (2) over angles we obtain the following equation for the energy density: This identifies the source in radiative transfer as S L (τ ) = v −1 E 1 0 dµI(0, µ) exp(−τ /µ) in terms of the incident radiation, and similarly for S R on the other side of the slab. Note that that is, the integral over optical depth of the energy source is determined by the incident flux density f + = F + /S. A useful identity can be obtained by defining K(τ ) ≡ 1 −1 dµµ 2 I(τ, µ). We easily see that ∂ τ K = −f . Since the total flux is conserved, K(τ ) = −f τ + const. We can check that this imposes the radiation I ∞ − 3 2 f (τ − µ) far from the boundary, with This relation connects the incident radiation on the surface, the reflected intensity and the one in the interior of the sample. For a conservative half space, f = 0, and the radiation pattern reaches the isotropic intensity I ∞ . The isotropic incident radiation I(0, µ > 0) = 2f in has incident flux density f in = F in /S, and Eq. (7) immediately gives us the (expected) result that I ∞ = 2F in /S, and hence the constant energy density w = 4F in /Sv E away from the boundaries. Previous work [19] showed that the delay time is essentially determined by a typical geometric length scale of the medium and largely independent from scattering details of bulk and surface, provided the incident radiation is isotropic. Pierrat et al. [20] confirmed this observation but emphasized the persisting role of transport velocity. Following this previous work, we can write for the total channel-averaged delay time t = T (t)t + R(t)t for an optically thick medium of arbitrary geometry with volume V and boundary surface S, For isotropic radiation incident on a slab we just derived that w = 4F in /Sv E away from the boundaries. This result looks "universal" and likely to be valid in more general geometries. Hence, by neglecting the surface layer, t = W/F in = 4V /Sv E which confirms the universal value ν = 1 reported in Refs. [19,20]. We recall that the energy transport velocity so far only appears phenomenologically in the dynamics of the EQRT, without any microscopic interpretation in terms of scattering properties. We can now provide a more microscopic definition. Equation (4) states that t = W/F in . Combining this with Eq. (8) yields for the energy transport velocity where, given isotropic incident radiation, w = W/V is the volume-averaged stored energy density in the medium and f in the incident flux density assumed constant across the boundary surface S. Because all quantities on the right-hand side are well defined and measurable, this equation can actually serve as the definition of the energy transport velocity, even in the localized regime. This definition is much in the spirit of the "energy velocity" as defined by Loudon [21]. It is also clear that via w the energy velocity becomes connected to the density of states in the medium. In the next section, we will find that, even in the localized regime with scale-dependent diffusion, W scales like the volume V so that ν is scaleinvariant though not always equal to 1 if the incident radiation is not isotropic.

III. EXAMPLES OF DELAY TIME CALCULATIONS
The most remarkable aspect of Eq. (4) is that total delay time-an intrinsic dynamic quantity-is related to the stored energy that can be found from the stationary EQRT (i.e. putting s = 0). In this section we will calculate total delay under different conditions. In the next section we will address delay in reflection and transmission separately, and will see that the dynamics then explicitly comes in.

A. Delay for isotropic incidence
For isotropic incident radiation equal on both sides of the slab we expect f = 0 and K(τ ) to be constant throughout the slab. Hence I is constant and isotropic everywhere. The main mathematical reason for this almost trivial result is that the orientational and depth dependencies of the specific intensity are strongly connected by EQRT. Absence of the first implies absence of the second.
Isotropic radiation I everywhere leads to constant energy density w = 2I/v E . In its turn, The result t = 2L/v E holds for the slab geometry but it is clear that the argument of constant energy density is valid for any geometry with a total surrounding surface S tot (= 2S for a slab) and a volume V , and even if the waves are localized by disorder. Hence ν = 1 in Eq. (8) and we recover the somewhat counterintuitive result by Blanco and Fournier [19]. If the number density n of the scatterers is small enough, we expect W = w 0 V + N W S , where W S is the total energy stored inside each scatterer, given the constant energy density w 0 outside. If v p is the phase velocity in the medium, we can write [11] which is the microscopic result, and scale independent. If the scatterer density is high enough for the waves to be localized, this result no longer applies. Nevertheless, W still scales with the volume V so that the definition (9) for v E does not reveal any scale dependence. A hand waving argument could have given W ∼ ξ 3 (with ξ the localization length) rather than W ∼ L 3 . The dependence of energy density on depth w(z) will be discussed in the next section and will help to explain why this argument is wrong.

B. Delay for normal incidence
For the incident wave normal to the surface z = 0 of a thick slab one finds the stationary reflection coefficient 3H(1, µ), which carries a unit flux density f − (0) = 1 0 dµµR(µ). At the same time, 1 0 dµµ 2 R(µ) = τ 0 , with τ 0 the extrapolation length in units of the mean free path [22]. Equation (7) then tells us that I ∞ = If we neglect the small energy contained in the boundary layers, we can impose W (B) = 0 so that the total flux density is f = (1+τ 0 )/B. From the Friedel identity, the total, channel-averaged de-lay time is Recall that dz = dτ (τ ) and that (τ ) = (B − τ ) if we assume that boundary conditions are identical on both sides of the slab. Thus, This result is the exact outcome of radiative transfer theory for normal incidence on a slab without internal reflections (τ 0 = 0.7104 . . .), for an arbitrary (symmetric) profile (z). We find here ν = 3(1 + τ 0 )/4 ≈ 1.28 in Eq. (8).

C. Delay in the diffusion approximation
In the diffusion approximation (DA) we replace EQRT (1) by the following diffusion equation: Here, G(z, z , s, q) is the Fourier-Laplace transformation of G(z, z , t, R) which stands for the energy density at time t, depth z and transverse distance R, given a source at t = 0, z , and R = 0. The solution for the energy density given an arbitrary incident radiation I(0, µ > 0, s, q) (where the q-dependence determines the transverse profile of the incident beam) is w(z, s, q) = dz S G(z, z S , s, q)S(z S , s, q) where we recall from the previous section that the source is S(z, s, where we omitted s and q as explicit arguments of G for brevity. For the total delay time we just need s = 0. We shall use simplified boundary conditions G = 0 at τ = B + τ 0 and τ = −τ 0 with B = L 0 dz / (z) the total optical thickness. This allows us to avoid mixed boundary conditions but generates small errors of order τ 0 /B in the energy density. For s = 0 and q = 0 the eigenfunctions are Φ n (τ ) = 2/B * sin[q n (τ + τ 0 )] with eigenvalues q n = πn/B * (we shall write B * = B + 2τ 0 and τ * = τ + τ 0 ), and the solution of Eq. (14) is The stored energy can be calculated as . We shall ignore the front factor that drops out in the delay time. Upon inserting the eigenfunction expansion (15) for G and assuming τ S B, we obtain And, finally, This result resembles closely Eq. (11) obtained from the radiative transfer theory. Note that we have assumed nothing yet about (z). The last simplification can be made if we assume that (z) = (L − z) which is true in the self-consistent theory of localization [4][5][6][7]. This implies τ (L − z) = B − τ (z) and Hence, to leading order, Thus DA yields ν = 3(τ S +τ 0 )/4 in Eq. (8). This is scaleindependent but depends on exact source depth and extrapolation length z 0 which, in turn, is known to depend on internal reflections on the sample surfaces. Without internal reflection at the boundaries τ 0 = 2 3 . We see that DA correctly reproduces the two cases considered within EQRT: the normal incidence has τ S = 1 and the isotropic incidence has τ S = 2 3 . For an incidence from direction µ 0 we easily find τ S = µ 0 .

D. Delay time for a sphere
In the following we will obtain the total delay time for a 3D sphere of radius R with equal radiation incident on all points of the outer surface that still may depend on µ. We expect the total delay time-integrated over all outgoing points on the surface-to be independent of angle, so we adopt a spherically symmetric source F in δ(r − r S )/4πr 2 . We shall put w(R + z 0 ) = 0 as a boundary condition in the diffusion equation Clearly, the outgoing flux is normalized, since −4πR 2 × 1 3 v E (R)dw/dr r=R,s=0 = F in . The Friedel identity applies and the delay time is equal to We can substitute dX = −dr/ 1 The solution of this equation satisfying w(X = B) = 0 and w(X → ∞) < ∞ to have finite energy in the center of the sphere is The delay time equals As follows from Eq. (26), for 0 < r < r S the energy density w(r) is constant and proportional to The second integral is a surface contribution that is smaller than the first one by a factor R/z S . Thus, for F in ∼ R 2 we find W ∼ R 3 . This "normal" scaling implies that This result is similar to what we have found for the slab. For normal radiation incident from the far field, and without internal reflection, the front factor equals 1 + 2/3 = 5/3, or equivalently ν = 1.25 in Eq. (8). For isotropic incident radiation the front factor equals 4/3, and we recover the universal value ν = 1.

IV. QUASI-1D TRANSPORT
In this section we investigate the delay time in a quasi-1D disordered wave guide to see what remains of the universal relation (8) if the measurement is done only in transmission or reflection. This question is extremely relevant for experiments where measuring both transmission and reflection may be problematic. We know that, in principle, in a quasi-1D wave guide with N transverse modes all waves are localized with a localization length ξ ∼ N B [24]. Here B is the transport mean free path in the absence of Anderson localization effects. In the self-consistent theory of localization, the quasi-1D geometry is described by the diffusion equation (14) with q = 0, supplemented by a self-consistent equation for the position-dependent transport mean free path [5,23]: which also depends on s. The Green's function of Eq. (14) is given by Eq. (15). For τ B, i.e. for a semiinfinite wave guide, the sum in Eq. (15) is essentially an integral and G(τ, τ ) = 3(τ + τ 0 ). We then easily find (z, The dynamics can be included using standard perturbation theory and treating s (τ, s)/v E as a small perturbation in Eq. (14). As in the previous section, we account for boundary conditions by introducing an extrapolation length z 0 , B * = B + 2τ 0 and τ * = τ + τ 0 . Since (τ, s) is suppressed by localization and since s is supposed to be a small hydrodynamic frequency, the perturbation can be argued to be"small". In first-order perturbation theory the eigenvalues of the diffusion equation change by and the eigenfunctions by Hence, We easily find that W nm is intrinsically a parameter of the dynamics, and it is of no relevance when s = 0. We can map the denominator in the sum of Eq. (32) onto the familiar diffuse It is tempting to conclude that each diffusion mode has now achieved its own diffusion constant D n (s = 0) = 1 We could proceed by saying that the transport velocity is affected by the dynamic kernel W nn / B ∼ ξ/B < 1 according to v (n) E ∼ v E /W nn and would thus be enhanced by the scaledependent diffusion. In this logic, the transport velocity would be associated with the energy density W ∼ (τ )W which is the conserved quantity featuring in the diffusion equation (14). Such a definition of v E is based on long-time tails of energy density in contrast to Eq. (9) that relies on the average delay time and the genuine energy density W which scales normally as W ∼ V even in the presence of scale-dependent diffusion. Apparently, the transport velocity is no longer uniquely defined when localization effects come into play. A deeper analysis is clearly necessary here but it is beyond the scope of this work.

A. Delay time in reflection
We are interested in calculating the weighted delay time in reflection: R(t)t = R(ω)dφ/dω = −dR(s)/ds at s = 0, with R(s) being the Laplace transform of the average reflection coefficient R(t). In the diffuse regime, it is easy to see that R(t)t ∝ L/v E up to a factor related to boundary effects, which provides a unique opportunity to measure the transport velocity directly. For normal diffusion, the energy density decays essentially algebraically between the two boundaries of the medium, and the total energy is thus proportional to L. In the following we ignore the s-dependence of imposed by localization effects, but we shall perform numerical calculations to make a comparison.
We consider a perfect point source δ(τ − τ S ) in Eq. (14). The average total reflection coefficient (integrated over all angles −1 < µ < 0) is It can easily be seen that with x = (τ S + τ 0 )/B * the transmission, exponentially small in the localized regime, and for a source near z = 0. We have two contributions to the delay time R(t)t = −∂ s R(s = 0). The first comes from the modi-fied eigenvalues: The second contribution to the delay time stems from the modified eigenfunctions, which generate Upon interchanging n and m in the first line we rewrite Eq. (38) as The approximation assumes τ S , τ 0 B. We can interchange x and 1−x and write (1−x) 2 where the delay in transmission T (t)t = 3ξτ * S (1 − L/B B )/v E is obtained in the next section. Anticipating this result gives us where all corrections of order z 0 /L have been ignored. This expression approaches (τ S + τ 0 )L/v E in the diffuse regime B ξ, and converges to 3 2 (τ S + τ 0 )L/v E deep in the localized regime. Contrary to the total delay, the delay time in reflection varies, though little, upon going from the diffuse into the localized regime. We illustrate this in Fig. 1(a) where results of different approaches to the calculation of R(t)t are compared. We see, in particular, that our perturbational result (41) is not exact and corresponds to a solution assuming (z, s) = (z, 0). Its dependence on the strength of localization effects quantified by the ratio ξ/L is, however, similar to the one exhibited by the exact solution of Eqs. (13) and (29).
A naive argument would suggest that in the localized regime, a wave penetrates only a distance ξ into the medium. This would lead to a much shorter delay time of order ξ/v E in reflection. This argument is apparently wrong by a factor L/ξ, because the energy does not decay as exp(−z/ξ). The lack of a penetration depth ξ in FIG. 2. Energy density, in units of F + /vE (with F + the incident flux) inside a disordered wave guide for a wave incident at z = 0, as a function of position z for four different values of the ratio of localization length ξ to the length L of the wave guide. We used z0 = 0 and zS = B for this figure. Independent of ξ/L, the total energy W inside the wave guide is always equal to 3(τS + τ0)LF + /2vE = 3LF + /2vE. the energy density w(z) is seen when we translate the (1 − τ /B)−profile of the energy density back to the z variable. We observe that w(z), in the localized regime, is actually flat on both sides of the medium with a steep descent in a region of size ξ around z = L/2. This is illustrated in Fig. 2 where we show profiles w(z) for several values of ξ/L. We observe that w(z) evolves from a linear decay in the diffuse regime ξ/L → ∞ to a steplike function deep in the localized regime ξ/L 1. For a given incident flux density F + , the energy density is always the same in the middle of the slab, whatever ξ/L. As counterintuitive as this can appear, the naive assumption of initial decay of w(z) over a region of size ξ is not confirmed by the z-dependent description of localization, since near the boundaries the waves are not localized. However, for a source in the middle of the wave guide, we find w(τ ) ∼ |τ − 1 2 B| − 1 2 B. The energy density as a function of z is then i.e. w(z) decays exponentially on both sides from the source and the waves are localized deep in the sample. The delay time, on both sides, then varies as t R,T = L 2 /D B for L ξ, and grows exponentially as exp(L/2ξ) deep in the localized regime.

B. Delay time in transmission
The transmission is given by and we easily obtain where we recall that x = τ * S /B * . Similar to Eq. (13) we find the first contribution to delay time in transmission, The second contribution is Again upon interchanging n and m in the first term, we obtain Hence, the mean weighted delay time is transmission is For ξ B B (diffuse regime), T (t)t = (τ S + τ 0 )L/2v E and thus the un-weighted delay time is t T = L 2 /6D B .
Upon entering the localized regime (ξ < B B ), this value saturates exponentially towards the L-independent value T (t)t = 3(τ S + τ 0 )ξ/v E . The un-weighted delay time is now equal to t T = 3ξB/v E . Naively, we could have expected t T = L 2 /6D with a reduced "scale-dependent" diffusion constant D = D B L/B B . This argument turns out to be wrong by a factor of order L/ξ. The dependence of T (t)t on ξ/L is illustrated in Fig. 1(b). Similarly to the case of reflection, we observe deviations of Eq. (48) from the exact numerical calculation, that also takes into account the dependence of on the dynamical parameter s. Quite remarkably, even though our results for both R(t)t and T (t)t are only approximate, their sum t = R(t)t + T (t)t is equal to the total delay time (12) exactly.

V. 3D SLAB
In this section we use the perturbation theory to study the delay time in a 3D slab. We also calculate other properties characteristic for 3D media: coherent backscattering, transverse spreading of a focused incident beam, and time-reversal focusing in the localized regime. As in section II, we consider a slab of disordered medium confined between planes z = 0 and z = L.

A. Delay time
The results obtained above for a quasi-1D system can be used-mutatis mutandis-for the 3D slab geometry, provided that we integrate over the transverse surface (hence q = 0). Physically this corresponds to measuring the delay times upon integration over the entire boundary surfaces. It can easily be seen that the total weighted delay in reflection, for example, is R(t)t = −dR(s = 0, q = 0)/ds. However, the diffusion coefficient is no longer determined by the return probability τ (1 − τ /B) formula valid in a quasi-1D wave guide. For q = 0, Eq. (13) remains valid and τ can be defined similarly in terms of (z). We define and find Let us consider a profile (z) = B /(1 + z/ξ c ) for 0 < z < L/2 (and mirrored on the other side of z = L/2), typically valid at the 3D mobility edge [4]. Since B dτ = dz(1 + z/ξ c ), we find with τ ξc = ξ c / B and B = L/ B + (L/2) 2 / B ξ c . For L ξ c we obtain R(t)t = (τ S + τ 0 )L/v E . For L ξ c , we have R(t)t = 3(τ S + τ 0 ) 23 60 (L 2 /4ξ c v E ) × 2ξ c /L. The correlation length ξ c drops out and R(t)t → 23 20 L/v E scales with the sample size L. The front factor 1.15 is smaller than the factor 1.5 obtained for the localized regime, and slightly exceeds the value 1 in the diffuse regime.
In transmission we proceed in a similar way and obtain This yields in the localized regime T (t)t → 7 20 (τ 0 + τ s )L/v E , smaller than but of the same order of magnitude as the weighted delay in reflection. Note that the Friedel identity (22) is obeyed (7/20 + 23/20 = 3/2). Again, the s−dependence affects both delays, but not their sum. In the diffuse regime, the ratio of weighted delays in reflection and transmission is 2 : 1, at the mobility edge this is close to 3 : 1. In the localized regime the ratio further grows up to L/ξ.

B. Transverse diffusion
In Ref. [12], dynamic transverse diffusion was used as a probe for Anderson localization. In the following we will treat 1 3 (τ, s) 2 q 2 as a small perturbation in Eq. (14) and study stationary properties. The Green function of the diffusion equation (14) is most conveniently written as G(τ, τ , q, s) = n Φ n (τ, q, s)Φ n (τ , q, s) where we have introduced Note that the Green function (52) looks as if anisotropic diffusion processes different for each mode were at work. The modification of eigenfunctions due to the perturba- Given a stationary source of waves of small transverse size at {x S = y S = 0, z S ∼ B }, the stationary energy at depth z and transverse distance R is ρ(z, R) = G(z, z = B , R, s = 0), with transverse Fourier transform ρ(z, q). The mean transverse energy spread at a depth z can be quantified by R(z) 2 = ρ(z, R)R 2 / ρ(z, R) = −(1/q)∂ q [q∂ q ρ(z, q)]/ρ(z, q) at q = 0 (we use R = {x, y}). Hence we need to expand to order q 2 . This expansion is similar to the expansion in s used to find the delay time, and we can copy the result (48) mutatis mutandis, Since ρ(τ, q = 0) = 3τ * S (1 − τ * /B * ) we find Near τ = B (in transmission) the second term is negligible, and we obtain In the diffuse regime = B and τ = z/ B so that R 2 (L) = 2L 2 /3. In the localized regime, (z) can be approximated by its profile in quasi-1D: B / = 1 + (τ /τ ξ )(1 − τ /B), so that These findings agree with previous results [23,25]. At a given optical depth τ B we can take the limit of the half-space (L → ∞) to see that and, in particular, R 2 (0) = 4ξz 0 in reflection. This result is surprisingly simple: the mean-square size of the transverse region in which the wave energy is concentrated grows linearly with the depth into the medium. For normal diffusion we also find a linear growth R 2 (z) = 4 3 L(z + z 0 )(1 − z/2L) with, however, a much larger slope that even diverges for a half-space since the transverse profile then becomes algebraic.

C. Coherent backscattering and time reversal
In the diffusion approximation and for normal incidence, the stationary CBS profile is approximately given by [26] C(Q) = G(τ = 1, τ S , s = 0, q = Q) with Q = k + k that vanishes at exact backscattering and G given by Eq. (52). For τ S ≈ 1, this yields the familiar formula for CBS of a normally incident plane wave. For τ S somewhere inside the slab, this expression actually describes the ensemble-averaged time-reversed profile by a perfect pointlike time-reversal machine at an optical depth τ S [27,28]. Using the perturbational approach of the previous section, we find the emerging specific intensity to be δC(Q) = − 6Q 2 B 2 π 3 (1 + τ 0 ) n,m V nm nm sin(q n τ * s ). (61) Recalling the definition of V nm , we simplify this to The background is given by sin q n (1 + τ 0 ) sin q n τ * S q 2 n (63) so that the normalized CBS profile becomes (τ S > 1) This result can be discussed in various limits. In the weak-disorder limit (z) = B and B = L/ B . Hence, the rounding of the CBS cone is typically −Q 2 τ S B L, i.e. it is rounded due to finite-size effects. For a half-space (L → ∞) this result is of little interest since the line profile is known to turn into the familiar cusp −|Q|τ S B , which is beyond the present perturbation theory as Q 2 has been assumed to be the leading small parameter. However, in the localized regime the limit of L → ∞ can be taken since the integral in Eq. (65) converges and we obtain C(Q) = 1 − 4Q 2 ξ(z S + z 0 ).
For CBS of a plane wave at normal incidence z S = B and we recover the rounding proportional to −Q 2 B ξ as predicted previously [4].
Finally, for a time-reversal experiment with a timereversal machine at depth z S we see, quite surprisingly, that the angular size of the focal spot δθ narrows down with (the genuine) depth according to δθ ∝ 1/k √ ξz s in contrast to δθ ∝ 1/kz s in the diffuse regime. We could have expected δθ ∼ 1/kξ, and arguably an impossibility to time-reverse well as z S > ξ, but this turns out to be a wrong expectation. Note, however, that we do not expect auto-focusing (i.e. focusing in the absence of ensemble averaging) to occur in the localized regime, because strong and long-range correlations should prevent self-averaging of a signal with even a relatively large bandwidth. A full discussion of this issue is beyond the scope of the present work.

VI. CONCLUSIONS
We have presented very general arguments based on the phenomenological equation of radiative transfer to derive simple expressions for the average delay time of a wave in a disordered medium. Our reasonings apply in all regimes of wave scattering, including the regime of strong (Anderson) localization. Specific examples of delay time calculations are provided for different geometries (a slab or a sphere) and different incident waves (an isotropic source or an incident plane wave). Detailed considerations of wave dynamics allowed us to suggest definitions for the energy transport velocity in the localized regime and to demonstrate that a unique definition for the latter may be difficult to achieve. In addition, we develop a novel perturbational approach to radiative transfer of localized waves in quasi-1D and 3D disordered media. This has enabled us to calculate the delay time measured separately in transmission or reflection which may be important to design experiments. We also apply our perturbation theory to study how well-known mesoscopic phenomena such as the transverse spreading of a wave packet, coherent backscattering, and time-reversal are affected by Anderson localization effects. A future study may be devoted to calculation of time-dependent quantities, such as the time-dependent transmission and reflection coefficients, in the framework of our perturbational approach.