Longitudinal modes in diffusion and localization of light
LLongitudinal Modes in Diffusion and Localization of Light
B.A. van Tiggelen ∗ and S.E. Skipetrov † Univ. Grenoble Alpes, CNRS, LPMMC, 38000 Grenoble, France (Dated: December 22, 2020)In this work we include the elastic scattering of longitudinal electromagnetic waves in transporttheory using a medium filled with point-like, electric dipoles. The interference between longitudinaland transverse waves creates two new channels among which one allows energy transport. Thispicture is worked out by extending the independent scattering framework of radiative transfer toinclude binary dipole-dipole interactions. We calculate the diffusion constant of light in the newtransport channel and investigate the role of longitudinal waves in other aspects of light diffusion byconsidering the density of states, equipartition, and Lorentz local field. In the strongly scatteringregime, the different transport mechanisms couple and impose a minimum conductivity of electro-magnetic waves, thereby preventing Anderson localization of light in the medium. We extend theself-consistent theory of localization and compare the predictions to extensive numerical simulations.
Keywords: light scattering, radiative transfer, Anderson localization
I. INTRODUCTION
The traditional and widely used picture of elastic mul-tiple scattering of light is one of a plane wave that expo-nentially extincts on the length scale of a mean free pathwhile propagating from one particle to another, with anelectric field orthogonal to the direction of propagation,and with subsequent scattering into a different directionin space, with exactly the same frequency. It is wellknown that transversality in real space ( r · E ( r ) = 0)is only valid in the far field of the scatterers, at distancesmuch larger than the wavelength. In the near field of adielectric object, the electric field achieves a “dipolar”structure, with a component directed along the prop-agation direction, while still being divergence-free, i.e. ∇ · E ( r ) = 0. In many approaches of multiple light scat-tering, these longitudinal modes are widely appreciated,yet considered “virtual”, in the sense that they do notcarry a Poynting vector so that they cannot transportelectromagnetic energy themselves, though they can me-diate the propagation of other waves, such as mechanical[2] or matter [3] waves.However, in inhomogeneous media, the dielectric con-stant ε ( r ) of the matter varies in space, and Gauss’ equa-tion imposes ∇ · [ ε ( r ) E ] = 0. As a result, true longitudi-nal electric fields exist, with ∇ · E (cid:54) = 0 and a finite densityof states (DOS) in phase space, to which elastic scatter-ing could take place. Induced polarization charges pos-sess Coulomb energy, and also stock dipole-dipole energyamong different scatterers but have no Poynting vector,so how can they transport energy? In atomic physics,the well-known process of F¨orster coupling [1] facilitatesa non-radiative transport mechanism to exchange quan-tum states and to move Coulomb energy from one atomto another. Like spontaneous emission, this process isinherently inelastic and incoherent, and is de facto ex- ∗ [email protected] † [email protected] cluded in a picture where only elastic multiple scattering,including interferences, is allowed. Much in the spiritof F¨orster coupling, Ref. [4] added explicitly the quasi-static dipole-dipole coupling as a new channel in trans-port theory of electromagnetic waves. In this work wewill show that this transport channel naturally emergesfrom a rigorous electromagnetic transport theory. The fi-nite Poynting vector of this channel is shown to originatefrom the interference between longitudinal and transversemodes.The transverse picture of electromagnetic wavesemerges naturally in the so-called “independent scatter-ing approximation” (ISA) of diffuse transport. In this ap-proximation, the longitudinal waves are usually ignored,and only transverse, propagating states are counted, as-sociated with damped plane waves with wave numbersclose to the frequency shell p ≈ k = ω/c in phase space.A fundamental question is whether this picture is signifi-cantly altered, within and beyond the ISA, or if just quan-titative modifications occur. Longitudinal states have afinite density of states (DOLS), proportional to the imag-inary part of the (longitudinal) dielectric constant of theeffective medium. Being mainly confined to scatterers,they exist far from the frequency shell in phase space,typically at very large wave vectors p (cid:29) k . We will showthat, due to the dipole size that is much smaller thanthe optical wavelength, excitations with large wave num-bers can scatter and mode-convert to both transverse andlongitudinal states. As such they take fully part in thediffuse transport.The purpose of this work is to include longitudinalwaves into the transport theory of electromagnetic wavesat scales well beyond the mean free path, identify sin-gularities, thereby respecting the conservation of energy.Finally we will perform the first study about the rolethat longitudinal electromagnetic waves play in weak andstrong (Anderson) localization. Static electric dipole cou-pling was already identified as a possible source of de-localization of mechanical waves [2]. Recent numericalsimulations with electromagnetic wave scattering frompoint-like electric dipoles revealed the absence of a mo- a r X i v : . [ c ond - m a t . d i s - nn ] D ec bility edge [5], and are difficult to explain within the tra-ditional picture that only acknowledges the transversefield as a mechanism for diffuse transport. II. TRANSPORT OF ELECTROMAGNETICWAVES
In standard transport theory [6], the dispersion and ex-tinction of waves are described by a complex self-energy Σ ( k, p ), associated with the effective medium. For elec-tromagnetic waves this is a second-rank tensor, depend-ing on frequency and wave vector p . Scattering betweentwo states in phase space is described by the four-rankscattering vertex U pp (cid:48) ( k ). Energy conservation in mul-tiple scattering is guaranteed by the Ward identity [6, 7], − Im Σ ( k + i(cid:15), p ) = (cid:88) p (cid:48) U pp (cid:48) ( k ) · − Im G ( k + i(cid:15), p (cid:48) ) (1)with the notation Im A ≡ ( A − A ∗ ) / i where A ∗ denotesthe Hermitian conjugate of a 3 × A ( A ∗ ij = ¯ A ji ).The left hand side stands for the extinction of an elec-tromagnetic excitation at wave vector p , the right handside puts this equal to the elastic scattering of the sameexcitation from p towards all other accessible states p (cid:48) inthe phase space. The “spectral tensor” − Im G ( k + i(cid:15), p (cid:48) )is positive (as (cid:15) ↓
0, for positive frequencies) and deter-mines the availability of microstates at the wave vector p (cid:48) , given the frequency ω = kc that is conserved inelastic scattering. For convenience we will drop explicitreference to (cid:15) and assume its presence in k + i(cid:15) implic-itly. Both Σ ( k, p ) and U pp (cid:48) ( k ) will be discussed in moredetail below. A. Dyson Green’s function
In Fourier space the Dyson Green’s tensor of an elec-tromagnetic “quasiexcitation” with frequency ω = kc and wave vector p of the effective medium is given by [6] G ( k, p ) = (cid:2) k − p ∆ p − Σ ( k, p ) (cid:3) − = ˆpˆp k − Σ L ( k, p ) + ∆ p k − p − Σ T ( k, p ) (2)split up into a longitudinal and a transverse part, with Σ ( k, p ) = Σ L ( k, p ) ˆpˆp + Σ T ( k, p ) ∆ p , with ∆ p = − ˆpˆp the projection tensor to transverse states. In transporttheory, the tensor G ( k, p ) ⊗ G ∗ ( k, p (cid:48) ) is the buildingblock of multiple scattering, and it is important to un-derstand G ( k, p ) in great detail on all scales.The longitudinal part of G ( k, p ) is associated with lo-cal Coulomb interactions between induced charges insidescatterers, often referred to as “non-radiative, static”,dipole-dipole coupling at a distance. The transverse partdescribes propagating waves. In the following we inves-tigate both components in real space for small and large distances. We demonstrate that at small distances, thelongitudinal part of the Dyson Green’s function domi-nates very generally and takes the form of dipole-dipolecoupling with the usual Lorentz contact term [28], andsurprisingly, is seen not to be static. At large distances,only transverse excitations contribute and G ( k, r ) is un-der very general conditions equal to an exponentiallysmall, propagating excitation with a polarization trans-verse to the direction of propagation r . This implies that G ( k, r ) contains the familiar near and far fields of elec-tromagnetism, without the need to add the first by hand[4]. At large distances the traditional picture, describedearlier, emerges.In real space, the Green’s tensor G ( k, r ) is the Fouriertransform of Eq. (2) and describes the propagation ofelectromagnetic waves over a distance r in the effectivemedium. The near-field component is “non-radiative” inthe sense that a longitudinal field E (cid:107) k induces no mag-netic field as k B ∼ k × E = 0. Alone, it carries thereforeno Poynting vector. However, we will show later in thiswork that the interference of longitudinal and transversecomponents in the tensor product G ⊗ G ∗ does carry aPoynting vector and facilitates a new channel to trans-port energy.With K L ( p ) ≡ k − Σ L ( k, p ) the square of a complexlongitudinal wave vector, one obtains in real space, G L ( k, r ) = (cid:88) p ˆpˆp K L ( p ) exp( i p · r )= δ ( r )3 K L ( ∞ ) + 1 − ˆrˆr πK L ( ∞ ) r + D ( r ) (3)where we have split off the singularity of the integral atlarge wave numbers, leaving the rest term D ( r ) as a con-tribution to the traceless dipole-dipole coupling describedby the second term. Since D ( r ) is, by construction, theFourier transform of a function that decays to zero forlarge p , it is free from a Dirac distribution, and evennon-singular as r →
0. We will show this explicitly insection II C for the recurrent scattering from two dipoles.As a result, the first two terms in Eq. (3) dominate onsmall scales. The first, subtle Lorentz contact term is agenuine Dirac distribution and vanishes for r (cid:54) = 0, butfor r = 0 makes a genuine contribution to DOS. Sincethe transverse field G T ( r ) ∼ /r for kr < G ( k, r → → G ,L ( K L ( ∞ ) , r ) (4)This takes the same form as the familiar dipole-dipoleregime of the bare Green’s function G ( r ), with how-ever the wave number k = ω/c in vacuum replaced bya complex-valued and frequency-dependent wave-vector K L ( ∞ ). For finite-size dielectric scatterers one may ar-gue that at small scales described by p → ∞ the effectivemedium is homogeneous and K L ( ∞ ) must be some real-valued wave number. For atomic atomic dipolar scatter-ers however, we will see that the complex value of K L ( p )extends up to infinity. The complex value of K L ( ∞ ) indi-cates that the dipole-dipole coupling, dominating in thenear field, is not static but depends on frequency andcontributes to the DOS. In section II D we will calcu-late K L ( ∞ ) numerically in all orders of the density for amodel of randomly positioned electric dipoles.At long distances kr → ∞ , small wave numbers prevailin Eq. (3) so that G L ( k, r → ∞ ) = 1 − ˆrˆr πK L (0) r (5)with K L ( p ) now evaluated at p = 0. If this propa-gator would not be compensated, the far field wouldcontain an algebraically small longitudinal term whichwould severely affect the random-walk picture of trans-verse electromagnetic wave transport. However, it iscompensated very generally by a part of the transversepropagator G T ( k, p ). For kr (cid:29) G T ( k, r ) = (cid:88) p ∆ p K T ( p ) − p exp( i p · r )= 12 π (cid:0) − ∇ + ∇∇ (cid:1) ir (cid:90) Γ dp e ipr p K T ( p ) − p + (cid:0) − ∇ + ∇∇ (cid:1) πK T (0) r (6)Here Γ denotes the line ( −∞ , + ∞ ) that avoids the origin p = 0 via a small contour in the upper complex p -plane,and which generates the last term. In the far field, sincenecessarily K T (0) = K L (0), the last term of Eq. (6) can-cels exactly against the longitudinal far field in Eq. (5).The Green’s function G ( k, r ) as a whole is therefore de-termined by the denominator of the first term and G ( k, r → ∞ ) = ∆ r π ir (cid:90) ∞−∞ dp p e ipr K T ( p ) − p (7)This indicates that the electric field is asymptoticallydominated by transverse modes and also transverse tothe direction of propagation r . If K T ( p ) has an analyt-ical extension at least over a small sheet Im p < K (cid:48)(cid:48) T inthe upper complex p -plane, G ( k, r ) will decay at least asexp( − K (cid:48)(cid:48) T r ) /r . Different “effective medium” approachesexist to calculate G ( k, r ) for various models [7]. The eas-iest method is to assume the presence of a simple pole K T ( p ) = k T + i/ (cid:96) , in which case normal exponentialbehavior emerges with the decay length equal to (twice)the elastic scattering mean free path (cid:96) .We conclude that the Green’s tensor of the ef-fective medium has a true longitudinal component( ∂ i G ij ( k, r ) (cid:54) = 0) that affects wave propagation at smallscales r < /k . In the far field, the electric field is alwaystransverse to propagation (ˆ r i G ij ( k, r ) = 0). Decay is ex-ponential under broad conditions with a decay length (cid:96) .This implies that radiative transfer should still be com-patible with a random walk with step length (cid:96) , thoughwith possibly new mechanisms for energy transport inthe near field provided by the presence of longitudinalfields, that can become dominant when k(cid:96) ≈
1. This idea will be worked out concretely in the next sectionsfor an ensemble of randomly distributed dipolar electricscatterers (“dipoles” for short).
B. Independent electric dipole scattering
In the independent scattering approximation (ISA) ap-plied to point-like electric dipole scatterers with num-ber density n and T -matrix t ( k ), Σ ISA ( k, p ) = nt ( k ).In this work we assume each dipole to be impenetra-ble for light outside, and to have only longitudinal ex-citations in its vicinity, at scales much smaller than thewavelength. This conveniently labels material energy aslongitudinal states that take part in the scattering pro-cess. By definition, the T -operator of a general polar-izable scatterer perturbs wave propagation in free spaceaccording to G ( k ) = G ( k ) + G ( k ) · T ( k ) · G ( k ). If weset T ( k ) = | r d (cid:105) t ( k ) (cid:104) r d | to describe an a electric dipole atposition r d , and impose (cid:104) r | G ( k ) | r d (cid:105) = 0 for any r (cid:54) = r d for it to be “impenetrable”, then it follows that t ( k ) = − (cid:104) r d | G ( k ) | r d (cid:105) = − (cid:34)(cid:88) p (cid:18) ˆpˆp k + ∆ p k − p + i (cid:19)(cid:35) − (8)This model can be refined to acknowledge finite penetra-tion of light into the dipoles [8], but the present choicehighlights the role of longitudinal waves and is arguablythe best description of elastic scattering from an atomwithout going into the details of atomic physics. Boththe longitudinal and the transverse integral diverge, thefirst essentially due to the Lorentz contact term. We willregularize the first as (cid:80) p ˆpˆp = 1 / u and the transversepart as (cid:80) p ∆ p /p = 1 / π Γ. It follows that t ( k ) = − π Γ k k − k − ik Γ (9)Both Γ (with dimension of length) and u (a volume)are genuine properties of the dipole, independent of fre-quency or polarization of the light. In particular k =2 π Γ /u determines the resonant frequency of the dipole.For k = k longitudinal and transverse singularities, op-posite in sign, cancel each other.For small k , the static polarizability α (0) is relatedto the t -matrix as t = − α (0) k [6], and we can iden-tify α (0) = 3 u . This relation can be understood fromclassical electrodynamics. We recall the Lorentz rela-tion E (0) = E − P for the homogeneous electric fieldinside the dipole, assumed spherical. Since we have im-posed E (0) = 0, the polarization density must equal 3times the local electric field E . The dipole moment isthus u P ≡ α (0) E = 3 u E with u the volume of thedipole, and hence α (0) = 3 u . The line width in fre-quency near the resonance is related to Γ according to γ = k c Γ = α (0) k c / π , a known relation for the ra-diative decay rate of a semi-classical two-level atom inthe electric-dipole approximation [9]. We can identifythe quality factor Q = ω /γ = 6 π/α (0) k . Near theresonance, we can thus write t ( k = ω/c ) = − πk γ/ ω − ω − iγ/ t -matrix satisfies the optical theorem, − Im t = (cid:88) p (cid:48) | t ( k ) | · ∆ p πδ ( k − p ) = | t ( k ) | k π This expression is consistent with Eq. (1), worked out lin-early in the dipole density n on both sides, with U ISA pp (cid:48) = n | t ( k ) | the ISA collision operator and Σ ISA ( p ) = nt ( k ).For its relative simplicity, many exact numerical simula-tions have been carried out with media filled randomlywith electric dipoles [5, 10–12], and many theoreticaltreatments exist already [13–15], not only because onecan go far without making further approximations butalso because they constitute a good and complete modelfor multiple scattering of light from simple atoms. Wenotice that the t -matrix of a single dipole is independentof both polarization, p and p (cid:48) . As a result, a single dipolecan scatter microstates with arbitrary state of polariza-tion, and with arbitrary p towards arbitrarily large p (cid:48) . C. Extinction involving two electric dipoles
The extinction caused by recurrent scattering from twodipoles was discussed in Ref. [15] for scalar waves, in Ref.[16] for low-energy electrons, and in Refs. [13, 17, 19]for electromagnetic waves. The last two works mainlyfocused on diffusion of transverse light, but used the fullGreen’s tensor (2) to describe recurrent scattering. InRef. [19] correlations between dipoles were included andcompared successfully to numerical simulations. Despitethe singular Green’s tensor G ( k, r ) in the near field, nonew divergencies were encountered provided the wholeseries of recurrent scattering is summed. In the followingsection we will explicitly include the longitudinal field inthe transport. To that end, we need to understand thebehavior of the self-energy tensor Σ ( k, p ) at large p . Theself-energy involving one or two different dipoles is givenby [13] Σ ( k, p ) = nt + n (cid:90) d r t G ( r ) − t G ( r )+ n (cid:90) d r t G ( r ) − t G ( r ) e i p · r + O (cid:0) n log n (cid:1) (11)We have dropped the explicit reference to k = ω/c in t ( k ) and in G ( k + i(cid:15), r ). The first term is the ISA,the second term involves recurrent loops between twodipoles. They are both independent of p and necessar-ily isotropic tensors. We will show in Sec. II D that, inour model, loop diagrams of arbitrary order rigorously −1−0.5 0 0.5 1 1.5 2 2.5 0 1 2 3 4 5 6 7 B oo m e r ang s e l f − ene r g y (r ea l ) momentum TransverseLongitudinal −0.8−0.6−0.4−0.2 0 0.2 0.4 0.6 0 1 2 3 4 5 6 7 B oo m e r ang s e l f − ene r g y ( i m ag i na r y ) momentum TransverseLongitudinal
Figure 1. Real (top, on resonance) and imaginary (bot-tom, for a detuning δ = ( ω − ω ) /γ = 1 .
5) part of thewave-number dependent “boomerang” self-energy associatedwith two dipoles. The wave number p (“momentum”) is ex-pressed in units of k = ω/c in free space, the self-energyis expressed in units of (4 πn/k ) × k . The transverseself-energy Σ T ( k, p ) converges asymptotically to zero (dashedline) for all detunings, meaning that the Lorentz local fieldterm Σ LL ( k, p ) = − n t , part of the boomerang diagramsbut independent of p , is canceled. The longitudinal self-energy Σ L ( k, p ), on the other hand, converges asymptoticallyto − n t (dashed line), as expressed by Eq. (13). The struc-ture visible near p = 3 k is caused by the triple round tripof light between two dipoles in the boomerang diagrams de-scribed by Eq. (11). determine the energy stored in longitudinal modes, andexploit this notion numerically. The third term, sum-ming up the so-called boomerang diagrams Σ B (see Fig.1) provides the first p -dependent contribution and causesΣ T ( p ) (cid:54) = Σ L ( p ). Higher orders in number density involve3 different dipoles or more. The boomerang diagramsgenerate a subtle contribution via the Lorentz contactterm δ ( r ) / k in G ( r ) [14, 18], which gives rise to thewell-known Lorenz-Lorentz correction − n t / k to boththe longitudinal and transverse dielectric functions, andthat is independent of p . Nevertheless, as p → ∞ , thisterm is compensated, again subtly, in the transverse self-energy and reappears as a purely longitudinal self-energy.This can be seen by subtracting the transverse photonfield G ,T ( r ) in free space, derived in Eq. (6) for the ef-fective medium, which is free from the Lorentz contactterm. The boomerangs become, Σ B ( p ) = n t (cid:90) d r (cid:20) G ( r )1 − t G ( r ) − G ,T ( r ) (cid:21) exp( i p · r )+ n t (cid:90) d r [ G ,T ( r ) − G ( r )] exp( i p · r ) (12)In the first term, the Lorentz contact term at r = 0 nolonger contributes and the integral vanishes for large p .The second integral is just equal to (minus) the longitu-dinal Green’s tensor G L, ( k, p ) in Fourier space. Hencewe find the somewhat surprising relation for infinite p ,lim p →∞ Σ ( p ) = Σ ISA + Σ
Loop − n t k ˆpˆp (13)Figure 1 illustrates by numerical integration of Eq. (11)that for large wave vectors, the Lorentz contact termis canceled in the transverse (boomerang) self-energy.It converges always to zero, whereas the longitudi-nal boomerang self-energy converges asymptotically toΣ L ( p ) = − n t /k . Neither one of them converges to − n t / k , associated with the Lorentz contact term.The asymptotic limit established in Eq. (13) is impor-tant since it demonstrates that K L ( ∞ ) (cid:54) = K T ( ∞ ), thefirst introduced earlier in Eq. (3) describing the dynamicdipole-dipole coupling in the near field.It is instructive to calculate the longitudinal Green’sfunction (3) associated with the self-energy in Eq. (11).Only the boomerang diagrams Σ B depend on wave num-ber p . Hence, up to order n , G L ( k, r ) = (cid:88) p ˆ p ˆ p k − Σ L ( p ) exp( i p · r )= − ∇∇ · (cid:88) p (cid:20) k − Σ + 1 k Σ B ( p ) (cid:21) exp( i p · r ) p with Σ = Σ ISA + Σ
Loop . Upon inserting the boomerangdiagrams and using − ∇∇ (1 / πr ) = δ ( r ) / − r ˆ r ) / πr = k G ,L ( r ), one obtains, G L ( k, r ) = G ,L (( k − Σ ) / , r )+ n k (cid:90) d r (cid:48) G ,L ( r − r (cid:48) ) · t G ( r (cid:48) )1 − t G ( r (cid:48) )This determines the longitudinal Green’s tensor at all dis-tances, and also depends on frequency for all distances.The first term stands for ordinary dipole-dipole couplingof the type 1 /r with a modified prefactor from the ef-fective medium that arises because we consider the elec-tromagnetic Green’s tensor and not the potential energyof the dipoles. The second term really changes the prop-agator from r = 0 to r (cid:48) , because a dipole can be situatedat r = 0, that first couples via a high-order dipole in-teraction to a dipole at r (cid:48) (a single coupling is alreadycounted in the effective medium) before finally arrivingat r . In the following we show 1) that this coupling fullydisappears at large distance (contrary to Ref. [4]) and2) that for small distances we recover the dipole-dipole coupling found earlier in Eq. (3), with the complex wavenumber K L ( ∞ ).For kr (cid:29)
1, we can take G ,L ( r ) out of the inte-gral, and recognize the remainder as the boomerang self-energy at p = 0. Hence, G L ( k, r → ∞ ) = G ,L ( K L (0) , r ) (14)This result agrees with Eq. (5) and was seen to cancelagainst a similar term in the transverse part of the DysonGreen’s function. For kr (cid:28) G L ( k, r →
0) = G ,L (( k − Σ ) / , r ) − n t k (cid:90) d r (cid:48) G ,L ( r − r (cid:48) ) · G ( r (cid:48) )+ n k (cid:90) d r (cid:48) G ,L ( r − r (cid:48) ) · t G ( r (cid:48) )1 − t G ( r (cid:48) )The last term is regular and r = 0 can be inserted. Thesecond term is equal to − ( n t /k ) G ,L ( r ) and adds upto the first term. Since by Eq. (13) we have K L ( ∞ ) = k − Σ + n t /k , G L ( k, r →
0) = G ,L ( K L ( ∞ ) , r ) + D (0)This agrees with Eq. (3) and attributes a finite complex,frequency-dependent value to D ( r = 0), D (0) = n t k (cid:90) d r (cid:48) G ,L ( r (cid:48) ) · G ( r (cid:48) )1 − t G ( r (cid:48) ) (15)We note that D (0) is negligible compared to the dipolarcoupling G L ∼ /r . D. Density of states
The total electromagnetic spectral density at frequency ω = kc in a polarizable medium is defined by N tot ( k ) = | k | c TR δ (cid:0) k − H (cid:1) with H = ε ( r ) − / ( p − pp ) ε ( r ) − / the Helmholtz oper-ator and TR the trace in the Hilbert space spanned by alleigenfunctions, including a strongly degenerate longitu-dinal eigenspace with eigenvalue 0. Written in this way,the spectral density is defined (and equal) for positiveand negative frequencies and normalized to the dimen-sion of the Hilbert space, (cid:90) ∞−∞ dω N tot ( k ) = TRindependent of ε ( r ), and formally infinite. We can workout the trace in real space as N tot ( k ) = (cid:90) d r | k | c (cid:104) r | Tr δ (cid:0) k − H (cid:1) | r (cid:105) with Tr the trace over 3 polarizations only, and identifythe integrand as the local density of states, N ( k, r ) = − kc π ImTr G H ( k + i(cid:15), r , r )with G H = [( k + i(cid:15) ) − H )] − . After ensemble-averagingit becomes independent of r , and we can express it interms of the Dyson Green’s function (2), (cid:104) N ( k, r ) (cid:105) = (cid:28) − kc π × ImTr (cid:104) r | ε / ( r ) · G ( k + i(cid:15) ) · ε / ( r ) | r (cid:105) (cid:29) = − kc π (cid:88) p Im Tr p ∆ p k · G ( k + i(cid:15), p ) (16)Both lines in this expression count, by construction, allstates but, quite surprisingly, the second line projects onthe transverse states only with however a large weight onlarge wave numbers p (cid:29) k . The reason is that the firstline counts electrical energy, including the longitudinalmodes, whereas the second line counts magnetic energy,that has only transverse modes. Equation (16) statesthat the density of states can be calculated from eitherthe magnetic or electrical energy, provided the latter in-cludes also the longitudinal states.For our model of electric dipoles we expect that theDOS is the sum of transverse traveling waves and stockedlongitudinal waves. To show this we go back to thefirst line of Eq. (16). For ε ( r ) = 1 + δε ( r ), we identify V = − δε ( r ) k as the interaction operator in the Bornseries of light scattering [6]. Before doing the config-urational average, we can consider M dipoles in a finitevolume V (see also Appendix A). Rigorous scattering the-ory imposes the operator identity V · G ( k ) = T · G ( k ).Hence N ( k, r ) = − kc π ImTr G ( k + i(cid:15), r , r )+ kc π ImTr (cid:104) r | T k · G ( k + i(cid:15) ) | r (cid:105) This equation is still exact and depends on the position r . Since the polarizability density δε ( r ) has disappearedexplicitly we can consider the special case of scatteringfrom identical, impenetrable electric dipoles, associatedwith a dielectric susceptibility δε ( r ) → ∞ , and describedby Eq. (8). For M such dipoles, T ( k ) = M (cid:88) mm (cid:48) T mm (cid:48) ( k ) | r m (cid:105)(cid:104) r m (cid:48) | (17)with, for m, m (cid:48) fixed, the 3 × T mm (cid:48) ( k ). To have G ( r m , r ) = 0 inside all dipoles at r m and for arbitrary r outside imposes that T mm (cid:48) ( k ) be given by the inverse ofthe 3 M × M matrix − G ( k, r m , r (cid:48) m ). It easily followsthat (cid:104) r | T · G ( k + i(cid:15) ) | r (cid:105) = − M (cid:88) m =1 δ ( r − r m ) = − n ( r ) −2−1 0 1 2 3 4 5 6 0 0.5 1 1.5 2 2.5 3 D O S (r) kr U P (r)U Q (r) −1.5−1−0.5 0 0.5 1 1.5 2 2.5 0 0.5 1 1.5 2 2.5 3 D O S (r) kr U P (r)U Q (r) Figure 2. The contribution of dipole-dipole coupling to theDOS as a function of distance between the dipoles. The vol-ume integral of the functions shown produces the second termof Eq. (22). U P ( r ) is associated with the electric field perpen-dicular to r , U Q ( r ) with the electric field directed along r .Top: δ = ( ω − ω ) /γ = − . δ = 0 . U P ( r ) has a subradiant peak only for positivedetuning whereas U Q ( r ) only for negative detunings. Since this is purely real-valued, it cancels in the ex-pression above for N ( k, r ). Upon averaging and letting M, V → ∞ at constant number density, the remainingterm yields (cid:104) N ( k ) (cid:105) = − kc π (cid:88) p ImTr G ( k + i(cid:15), p ) (18)in terms of the Dyson Green’s function (2). This is rec-ognized as (cid:104)| E ( r ) | (cid:105) , proportional to the energy density (cid:104) E ( r ) (cid:105) / π , averaged over disorder and cycles, and hav-ing both longitudinal N L ( k ) and transverse N T ( k ) parts.We emphasize that Eq. (18) only applies for our modelthat excludes any light inside the scatterer. As a resultno stored energy density (cid:104) E ( r ) · P ( r ) (cid:105) exists as e.g. inMie scattering [6]. In this model, the stocked dipole-dipole energy is entirely described by longitudinal (elec-tric) waves.We can insert the Dyson Green’s function obtained in −3−2.5−2−1.5−1−0.5 0 0.5 1 1.5−1.5 −1 −0.5 0 0.5 1 1.5 D O S Detuning [line width]
DOS(P)DOS(Q)DOS
Figure 3. The contribution of dipole-dipole coupling to theDOS (in units of N × Q × (4 πn/k ) ) as a function of de-tuning δ = ( ω − ω ) /γ . The dashed lines show the separatecontributions of the modes with electric field perpendicular( P ) and parallel ( Q ) to the separation vector r . Sec. II A into Eq. (18), (cid:104) N ( k ) (cid:105) = − kπc Im (cid:34) K L ( ∞ ) (cid:88) p + (cid:88) p (cid:18) K L ( p ) − K L ( ∞ ) (cid:19) +2 (cid:88) p K T ( p ) − p (cid:35) Only the first term, stemming from the singular Lorentzcavity and entirely governed by longitudinal excitations,diverges and is regularized using (cid:80) p = 3 /α (0) = k Q / π consistent with section II B. The second term,proportional to D ( r = 0) in Eq. (3), is non-zero but afactor Q smaller and shall be neglected. Finally the lastterm is just the density of states of transverse waves. Weshall assume the existence of a well-defined complex pole K T = k T + i/ (cid:96) . This gives (cid:104) N ( k ) (cid:105) = k π c (cid:20) − Q Im k K L ( ∞ ) + k T (cid:21) (19)The ratio of longitudinal and transverse LDOS is thus (cid:104) N L ( k ) (cid:105)(cid:104) N T ( k ) (cid:105) = − Q k k T Im 1 K L ( ∞ ) (20)In view of the factor Q this can be a large number,proportional to the density of the dipoles. For low densityis K L ( ∞ ) ≈ K T ≈ k + i/ (cid:96) so that (cid:104) N L (cid:105) / (cid:104) N T (cid:105) = Q /k(cid:96) .This ratio will be discussed in the next sectionA rigorous expression can be derived for DOS withoutrelying on the existence of a complex pole by taking intoaccount the p -dependence of the self-energy Σ T ( k, p ) as-sociated with the scattering from two electric dipoles. A direct expansion in dipole density yields (cid:104) N ( k ) (cid:105) = − kπc ImTr (cid:88) p (cid:104) G ( k, p )+ G ( k, p ) · Σ ( k, p ) · G ( k, p )+ G ( k, p ) · Σ ( k, p ) · G ( k, p ) · Σ ( k, p ) · G ( k, p ) (cid:105) + O ( n )Several singular longitudinal terms, stemming from theLorentz cavity can be seen to cancel. The first termdescribes the free electromagnetic field and the longi-tudinal field drops out trivially. The longitudinal com-ponent of the second term contains a singular Lorentzcavity ( nt + Σ Loop − n t /k ) (cid:80) p G L ( p ) stemming fromEq. (13). Similarly, the third term generates a singu-lar longitudinal contribution n t (cid:80) p G L ( p ) that cancelsexactly against the local field − n t /k generated by theprevious term. We can work out the wave number in-tegral in the expression for (cid:104) N ( k ) (cid:105) exactly by insertingEq. (11), and use the cyclic property of the trace, (cid:104) N ( k ) (cid:105) = k π c + kπc ImTr (cid:20) nt ∂∂k G ( k, n t (cid:90) d r G ( r ) − t G ( r ) · ∂∂k G ( k, n t (cid:90) d r G ( r ) − t G ( r ) · ∂∂k G ( k, r )+ n t (cid:90) d r G ( k, r ) · ∂∂k G ( k, − r ) (cid:21) We have transformed the integral over wave vectors p ofthe last term G · Σ · G · Σ · G in (cid:104) N ( k ) (cid:105) to real space.Using again the relation /t ( k ) = − G ( k, r = 0) this canbe rearranged to (cid:104) N ( k ) (cid:105) − N ( k ) = − n π ddk Im log t − n πc ddk ImTr (cid:90) d r log (cid:2) − t G ( r ) (cid:3) (21)with N = k / π c the LDOS of transverse waves infree space. The appearance of a full frequency derivativein the DOS is a manifestation of Friedel’s theorem [20].The second term is recognized as the dipole-dipole energyexpressed as the “return trip operator”, widely used inthe theory of Casimir energy in matter [21], and involvesloop paths only. The integral is well-defined at both r = 0and r → ∞ . Since the dominating frequency dependencecomes from dt/dk ≈ − Q t / π , (cid:104) ∆ N ( k ) (cid:105) N ( k ) = − Q k ImTr (cid:20) nt + n (cid:90) d r t G ( r )1 − t G ( r ) (cid:21) = − Q k Im (Σ
ISA + Σ
Loop ) (22)This expression suggests that in general the modificationof DOS is dominated by ISA + loop diagrams, describing (cid:16) (cid:22) (cid:16) (cid:21) (cid:16) (cid:20) (cid:19) (cid:20) (cid:21) (cid:22)(cid:19)(cid:17)(cid:19)(cid:19)(cid:17)(cid:21)(cid:19)(cid:17)(cid:23)(cid:19)(cid:17)(cid:25)(cid:19)(cid:17)(cid:27)(cid:20)(cid:17)(cid:19) (cid:16) (cid:22) (cid:16) (cid:21) (cid:16) (cid:20) (cid:19) (cid:20) (cid:21) (cid:22) (cid:16) (cid:20)(cid:17)(cid:24) (cid:16) (cid:20)(cid:17)(cid:19) (cid:16) (cid:19)(cid:17)(cid:24)(cid:19)(cid:17)(cid:19)(cid:19)(cid:17)(cid:24)
Figure 4. Top: Numerical simulation of the averaged imagi-nary part of the diagonal elements of the T -matrix as a func-tion of detuning, for several densities n of the dipoles. ViaEqs. (19) and (23) this quantity determines the longitudinaldensity of states LDOS. Bottom: Same but with the ISAapproximation subtracted and normalized by 4 πn/k . Thedashed line shows the second-order in density term in Eq. (22)as also shown in Fig. 3. longitudinal excitations, even if mediated by transverse,propagating waves. This, in turn, implies that the com-plex longitudinal wave number K L ( ∞ ) is governed byloop diagrams only. In Appendix A we demonstrate thatthis statement holds rigorously. More precisely, if we re-call the T -matrix (17) of M electric dipoles randomlydistributed in a volume V , then1 K L ( ∞ ) = 1 k + 13 k n (cid:104) Tr T mm ( k ) (cid:105) (23)for M → ∞ at constant M/V = n . The first two termsin the density expansion clearly coincide with Eq. (22).All higher order terms are rigorously loop diagrams andthe ensemble-average of the diagonal element T mm overall other M − (cid:104) ∆ N (cid:105) /N hasa Lorentzian profile with a large peak height inverselyproportional to Q , to be associated with the excita-tion of a single dipole. The second term of Eq. (22)becomes important when 4 πn/k ≈ G ( k, r ) = − exp( ikr ) / (4 πr )[ P ( kr ) ∆ r + Q ( kr )ˆ r ˆ r ] for r (cid:54) = 0 [6] the integrand can be split up into two interac-tions U Q ( r ) and U P ( r ), that govern the near-field cou-pling of two dipoles in real space. Both are shown inFig. 2. The total dipole-dipole coupling, shown in Fig.3, is negative around the resonance. Because local fieldsingularities cancel in the DOS, the solid curve in Fig. 3is the same as found in Ref. [17].In Fig. 4 we have calculated numerically the diagonalelements of the 3 M × M matrix T mm for M = 10 dipoles homogeneously distributed in a sphere at den-sity n = M/V , thereby averaging over all 3 M diagonalelements as well as over 10 independent random config-urations of the dipoles. This calculation confirms that K L ( ∞ ) is a genuine complex quantity and in general dif-ferent from the complex wave vector K T = k T + i/ (cid:96) as-sociated with the transverse modes (see also Figs. 11 and12 in Appendix B). For low dipole densities, the calcula-tion agrees accurately with the analytical calculation ofthe loops between two dipoles in Eq. (22). The line pro-file of the longitudinal DOS broadens significantly wellbeyond the single-dipole line profile as the dipole densityincreases. Nevertheless, the total surface underneath re-mains constant. This is to be expected since each dipolecontributes exactly 3 microstates to the DOS and thisnumber cannot be affected by dependent scattering (seeAppendix A). E. Equipartition between longitudinal andtransverse waves
If we acknowledge the finite density of longitudi-nal states (DOLS) proportional to − Im G L ( k, p ) ≈− n Im t ( k ) /k , the right hand side of Eq. (1) allows thescattering towards longitudinal states with arbitrary p (cid:48) ,even in independent single scattering. The scattering isindependent of initial wave number p , but of course re-quires the existence of this initial state, governed by thespectral function − Im G L,T ( k, p ). For transverse wavesthe rate is proportional to the occupation number n T ( p )of the transverse states with wave number p and the totalavailable number of longitudinal states, ρ T ( p ) /τ ISA T → L ∼ n T ( p )( − )Im G T ( p ) × (cid:88) p (cid:48) U ISA · ˆp (cid:48) ˆp (cid:48) ( − )Im G L ( p (cid:48) )= − Q n | t | Im t πk × n T ( p )( − )Im G T ( p ) (24)Since the integral diverges at large p (cid:48) , we have used thesame regularization as the one employed earlier for the T -matrix of one dipole. The resulting scattering rate ispositive, converting the transverse state with wave vector p to available longitudinal state with mostly large wavevector p (cid:48) .The matrix element U (2) pp (cid:48) involving recurrent scatteringfrom two dipoles [17] can mode-convert any initial stateto transverse states, with a rate proportional to n , as isthe process in Eq. (24). A close look identifies only one event part of U (2) pp (cid:48) that gives rise to a singular scatteringrate (see Fig. 5). The irreducible ladder diagrams, aspart of U pp (cid:48) without external lines, add up to U (LAD) pp (cid:48) = n | t | (cid:90) d r (cid:20) G − t G (cid:18) G − t G (cid:19) ∗ − G G ∗ (cid:21) (25)Our tensor notation is ( AB ) ij | kl = A ik B lj , equivalentto ( AB ) · S = A · S · B . This event scatters again indepen-dent of p and p (cid:48) . The second term must be subtractedsince it stands for a reducible event that is not part ofthe collision operator. However, the substraction createsa diverging contribution at r = 0 in the integral due tothe singular longitudinal field. To repair this in a wayconsistent with previous sections, we extract the trans-verse photon propagator G ,T ( r ), and write U (LAD) pp (cid:48) = n | t | (cid:90) d r (cid:20) G − t G (cid:18) G − t G (cid:19) ∗ − G ,T G ∗ ,T (cid:105) + n | t | (cid:88) p (cid:48)(cid:48) (cid:2) G ,T ( p (cid:48)(cid:48) ) G ∗ ,T ( p (cid:48)(cid:48) ) − G ( p (cid:48)(cid:48) ) G ∗ ( p (cid:48)(cid:48) ) (cid:3) We have used Parseval’s identity to convert the secondintegral into an integral over wave vectors. The first termof U ( L ) pp (cid:48) now converges at small r , the second term canbe dealt with as before, giving U (LAD ,s ) pp (cid:48) = − Q n | t | πk S + Q n | t | πk (cid:18) − S (cid:19) (26)with S ≡ (cid:104) ˆpˆpˆpˆp (cid:105) the fully symmetric four-rank tensor.The first term of this collision operator can convert bothtransverse and longitudinal waves to available transversewaves with wave number p (cid:48) . The rate is proportional to ρ T ( p ) /τ LAD T → L ∼ n T ( p )( − )Im G T ( p ) (cid:88) p (cid:48) U (LAD ,s ) pp (cid:48) · ∆ p (cid:48) ( − )Im G T ( p (cid:48) )= − Q n | t | (6 π ) × n T ( p )( − )Im G T ( p ) (27)Since t satisfies the optical theorem, the two secular scat-tering rates (24) and (27) compensate each other on theleft hand side of the Ward identity (1), consistent with atotal scattering rate for transverse waves that was seenin Eq. (11) not to suffer from any singularity.Nevertheless, the two scattering events produce each adifferent polarization and thus affect the balance of lon-gitudinal and transverse energy. In the Bethe-Salpeterequation—to be introduced in Eq. (33) of the nextsection—the scattering rate for longitudinal states to Figure 5. Diagrammatic presentation of the ladder se-ries (without external lines) involving two different electricdipoles. Dashed lines connect identical dipoles, solid linesdenote the Green’s tensor G ( r ), crosses denotes transitionmatrix t ( k ), bottom line denotes Hermitian conjugation. Thefirst diagram on the left is reducible (a simple product) andis not part of the collision operator U pp (cid:48) . convert back to a transverse mode with wave number p must have the same matrix element as in Eq. (24) butnow depends on the total number of occupied longitudi-nal states n L ( p (cid:48) ), and the number of available transversestates, ρ L ( p ) /τ ISA L → T ∼ ( − )Im G T ( p ) × (cid:88) p (cid:48) U ISA · ˆp (cid:48) ˆp (cid:48) n L ( p (cid:48) )( − )Im G L ( p (cid:48) ) (28)This conversion, again arbitrarily fast for longitudinalexcitations with arbitrarily large p (cid:48) , does not cancel ingeneral the reverse process of Eq. (24) unless n T ( p ) = n L ( p ) = n for all p . This implies that the occupied phasespace is proportional to the number of available statesand equipartition of energy has been reached. These sin-gular terms therefore lead to a very fast equipartition be-tween longitudinal and transverse modes in phase space,and once established, stay constant and cancel in the suc-ceeding transport dynamics. In this regime, the ratio oflongitudinal and transverse energy densities is (cid:104) ρ L (cid:105)(cid:104) ρ T (cid:105) = N L ( k ) N T ( k ) = Q k(cid:96) (29)The rate for LT-mode conversions associated with thesingular event can be compared to the rate that governsthe equipartition among wave vectors of traveling trans-verse waves, and is in the same units proportional to U ISA × k/ π = k/(cid:96) . Hence τ − T L τ − S = 13 Q k(cid:96) (30)If k(cid:96)/Q (cid:28) (cid:104) E L (cid:105) (cid:29) (cid:104) E T (cid:105) and 1 /τ T L (cid:29) /τ S . This means that longitudinal states dominate inenergy, and equilibrate in phase space as fast as the trans-verse waves. In fact, when k(cid:96)/Q <
1, the intermediatescattering to a longitudinal wave becomes more efficientfor transverse waves to equilibrate among themselvesthan accomplished by ISA single scattering. The atomicquality factor Q is large, and experiments [22] and nu-merical simulations [5, 10, 12] exist where Q (cid:29) k(cid:96) .0 III. KUBO FORMALISM
In this section we use the rigorous Kubo formalismfor the DC conductivity, adapted from electron conduc-tion [20] to scalar classical waves [23] and electromag-netic waves [24]. We investigate how photon diffusion isaffected by the existence of longitudinal waves.Before averaging over the disorder, the electric field atfrequency ω = kc is given formally by the operator iden-tity E ( k ) = G ( k ) ⊗ s ( k ), with s ( k ) a source ( ⊗ stands forthe matrix product in full Hilbert space, whereas · standsfor matrix product in 3 × φ ij = (cid:104) E i ¯ E j (cid:105) ofthe electric field at two different frequencies and for twodifferent positions. We can formally relate it to the sourcecorrelation function S according to φ ij = R ij | kl ⊗ S kl ,which introduces the reducible four-rank vertex R . Itsatisfies the Bethe-Salpeter equation, R = GG † + GG † ⊗ U ⊗ R (31)This equation identifies the irreducible vertex U as thescattering operator, and GG † as the transport betweenscattering events. (We use † for Hermitian conjugatein full Hilbert space as opposed to ∗ for Hermitian con-jugate of a 3 × AB · S = A · S · B , orequivalently ( AB ) ij | kl = A ik B lj , with the matrix B dis-played as the bottom line of a Feynman diagram, prop-agating backwards in time. Similarly in Hilbert space,( GG † ) αβ | κγ = G ακ G † γβ = G ακ ¯ G βγ . After averaging,translational symmetry can be exploited so that the ver-tex in Fourier space (Fig. 6) can be written as R pp (cid:48) ( q ),with p (cid:48) and p interpreted as incident and outgoing wavenumbers, and q conjugate to distance between source andobserver. Thus, the electromagnetic “Wigner function”takes the form φ ij ( p , q ) ≡ (cid:104) E i ( p + ) ¯ E j ( p − ) (cid:105) = (cid:88) p (cid:48) R pp (cid:48) ; ij | kl ( q ) S kl ( p (cid:48) , q )(32)and R pp (cid:48) ( q ) = G ( p + ) G ∗ ( p − ) δ pp (cid:48) + G ( p + ) G ∗ ( p − ) · (cid:88) p (cid:48)(cid:48) U pp (cid:48)(cid:48) ( q ) · R p (cid:48)(cid:48) p (cid:48) ( q ) (33)with p ± = p ± q / δ pp (cid:48) ≡ (2 π ) δ ( p − p (cid:48) ).The two terms describe direct propagation with extinc-tion of the mode p and scattering from p (cid:48) towards p , re-spectively. One important property is reciprocity [24, 25].Since without external magnetic fields, the (unaver-aged) Green’s function satisfies (cid:104) p , i | G ( k + i | k, p (cid:48) (cid:105) = (cid:104)− p (cid:48) , k | G ( k + i | i, − p (cid:105) we easily check that R ij | kl, pp (cid:48) ( q ) = R kl | ij, − p (cid:48) − p ( − q ) (34) Figure 6. The diagrammatic convention associated with thereducible vertex R ij | kl, pp (cid:48) ( q ), with external lines. Top linedenotes retarded Green’s function G ( k + i G † ( k + i
0) = G ( k − i
0) is the advanced Green’s functionand travel in the opposite direction. The polarization labelsare ij on the left hand side (“observer”) and kl on the righthand side (“source”).The sum of incoming and outgoing wavenumbers is conserved. A second property follows from complex conjugation,equivalent to switching bottom and top lines of the dia-gram, R ij | kl, pp (cid:48) ( q ) = ¯ R ji | lk, pp (cid:48) ( − q ) (35)If Eq. (1) is satisfied, R is known to exhibit long-range dif-fusion ( q → R ij | kl, pp (cid:48) ( q ) = d ij ( p , q ) d kl ( p (cid:48) , q ) πN ( k ) D ( k ) q (36)where N ( k ) is the DOS given in Eq. (16) and the eigen-function associated with long-range diffusion is writtenas d ( p , q ) = − Im G ( p ) − i J ( p , q ) + O ( q ) (37)The first term − Im G ( p ) ≡ − [ G ( p ) − G ∗ ( p )] / i isproportional to the spectral function and implies per-fect equipartition of the electromagnetic energy in phasespace. The second term is linear in q and describes asmall perturbation due to gradients of Φ ij ( r ) in real spacethat trigger diffuse energy flow. For it to be small for allmomenta p imposes constraints to be discussed later. Be-cause d ( p , q ) describes an electric field correlation func-tion, it must satisfy d ij ( p , q ) = ¯ d ji ( p , − q ), consistentwith Eqs. (35) and (34). Thus, J ( p , q ) = − J ∗ ( p , − q )and, being linear in q by construction, we conclude thatthe tensor J ( p , q ) is Hermitian.Following common treatments in radiative trans-fer, many microscopic approaches interpret the expan-sion (37) as one in the angular anisotropy of scatteredradiation with wave numbers in equipartition and im-posed near the frequency shell, as described by the firstterm. If we ignore electromagnetic polarization and with-out any kind of explicit anisotropy in space, the only1possible choice of this expansion is, d ( p , q ) = − Im G ( p ) [1 − iJ ( p ) p · q + · · · ] (38)For diffusion of cold atoms, J ( p , q ) was obtained by solv-ing numerically the Bethe-Salpeter equation [26]. Alter-natively, the unknown function J ( p ) can be chosen suchthat the first angular moment (cid:80) p pd ( p , q ) matches thedivergence − i q · K of the energy current density. Thisleads to J ( p ) = 1 /p [7] and makes the vector K the onlyunknown. This choice conveniently circumvents diver-gencies that occur in rigorous theory for large p . For vec-tor waves, J ( p , q ) is a tensor containing longitudinal andtransverse components, and even their interferences. Wewill derive it more rigourously from the Bethe-Salpeterequation (33) whose orders linear in q give us J ( p , q ) = J D ( p , q )+ G ( p ) G ∗ ( p ) · (cid:88) p (cid:48) U pp (cid:48) · δ q G ( p (cid:48) , q )+ G ( p ) G ∗ ( p ) · (cid:88) p (cid:48) U pp (cid:48) · J ( p (cid:48) , q ) (39)The first term is often referred to as the Drude con-tribution to diffusion and depends only on the effectivemedium properties. It reads J D ( p , q ) = G ( p ) · L ( p , q ) · G ∗ ( p ) − δ q Re G ( p , q ) (40)in terms of the bilinear Hermitian tensor L ij ( p , q ) = 2( p · q ) δ ij − p i q j − q i p j and the notation is δ q Re G ( p , q ) =( q · ∂ p )Re G ( p ). The second and third terms in Eq. (39)are genuine contributions from scattering. They vanishonly for isotropic events in U pp (cid:48) but not in general. It isstraightforward to demonstrate that the (cycle-averaged)Poynting vector K = c Re ( E × ¯ B ) / π is related to thecorrelation function of the electric field according to K n ( k, q ) = c πk (cid:88) p (cid:18) p n δ ik − p k δ in − p i δ kn (cid:19) φ ki ( p , q )+ c πk (cid:88) p q k
12 ( φ kn ( p , q ) − φ nk ( p , q )) (41)In the absence of external magnetic fields, φ ki ( p , ) = φ ik ( p , ) so that the second term vanishes in linear orderof q . Upon inserting Eq. (37), the first term involves onlythe diffusion current tensor J . Some manipulations leadto iq n K n = 14 πN ( k ) (cid:88) p L ik ( p , q ) J ki ( p , q ) × D ( k ) q c πk (cid:88) p (cid:48) − Im G lj ( p (cid:48) ) · S lj ( p (cid:48) ) (42)The factor that has been split off on the right hand sidecan be identified as the (cycle-)averaged energy density ρ ( r ) = (cid:104)| E | + | B | (cid:105) / π released by the source and dif-fusing out. This can be established by noting that the tensor J , being odd in p , does not contribute to theenergy density. Since the first term in Eq. (37) obeysequipartition, (cid:104)| E ( r ) | (cid:105) = (cid:104)| B ( r ) | (cid:105) , the energy densityis (cid:104) ρ ( q ) (cid:105) = 1 πN ( k ) kc Tr (cid:88) p p k ∆ p · − Im G ( k, p ) × D ( k ) q c πk Tr (cid:88) p (cid:48) − Im G ( p (cid:48) ) · S ( p (cid:48) ) (43)If we recall that the electromagnetic DOS is given byEq. (16), the first factor in Eq. (43) equals one. In realspace Eq. (42) thus becomes ∇ · K = − D ∇ ρ ( r ) with πN ( k ) D ( k ) = 14 Tr (cid:88) p L ( p , ˆ q ) · J ( p , ˆ q ) (44)This is the Kubo formula for the electromagnetic diffu-sion constant. Since D is a scalar in this work, the righthand side does depend on the direction of q . The lefthand side can be identified as the electromagnetic DCconductivity σ ( k ) (here in units of 1 / m) expressed as the(Einstein) product of DOS and diffusion constant. Withthis definition, that we prefer in view of the presence of πN D in Eq. (36), the ”electromagnetic conductance” ofa slab with surface A and length L takes the form ofa Landauer formula (cid:104) (cid:80) ab T ab (cid:105) = 4 σA/L [27]. In termsof the energy density ρ ( q ) the electric field correlationfunction is φ ij ( p , q ) = d ij ( p , q ) πN ( k ) × πkc ρ ( q ) (45) A. Diffusion current tensor
The diffusion current tensor J ( p , q ) must be a parity-even, Hermitian tensor, linear in the gradient vector q .For our problem, with no explicit anisotropy present, thisleaves us with the following general form J ( p , q ) = J ( p )( p · q )∆ p + J ( p )( p · q )ˆ p ˆ p + J ( p )( pq + qp ) + J ( p ) i ( pq − qp ) (46)with four real-valued functions J i ( p ) to be determined. Afifth term i(cid:15) ijk q k is in principle allowed but is excludedfor scattering that respects parity symmetry. Alterna-tively, we could have defined the mode J ( p ) in terms ofthe tensor pq + qp − p · ˆ q )ˆ p ˆ p in which case all 4 modeswould be mutually orthogonal.The four functions can be associated with four differentaspects in diffuse transport. By restricting only to thefirst, the transport problem reduces to the common ap-proximation made in Eq. (38). The modes J , J and J are clearly genuine vector effects, absent in a scalar the-ory. However, only J and J carry a Poynting vector,with J associated with the transport of transverse wavesin the far field, and J associated with a novel process2that involves the interference of longitudinal and trans-verse waves. By restricting to the purely transverse term J ( p ), transport theory almost reduces to a scalar theory.The term J describes how the longitudinal energy den-sity | E L ( p ) | achieves an anisotropy in phase space dueto the spatial gradient of energy, but without inducingan energy current. The presence of J is more subtle andcan be associated with the imaginary part of the com-plex Poynting vector, discussed for instance in Ref. [28].Let us call Im K = c Im ( E × ¯ B ) / π . We readily find,similar to the derivation of its real part in Eq. (41), thatin terms of the field correlation function φ ik ( p , q ),Im K n ( k, q ) = − ic πk (cid:88) p (cid:18) q n δ ik − q k δ in − q i δ kn (cid:19) φ ki − ic πk (cid:88) p p k ( φ ki − φ ik ) (47)The first term is independent of J and can be evaluatedwithout any approximation. The integral over wave num-bers is proportional to the total DOS N ( k ) and cancelsthis same factor in the denominator of Eq. (45). As aresult it is completely independent of the presence of thedipoles. The second term requires anti-symmetry in thediffusion tensor J ij , described only by J ( p ). We obtainIm K ( k, q ) = 23 (cid:32) c k + 1 πN ( k ) (cid:88) p p J ( p ) (cid:33) ( − i q ) ρ ( q )(48)Like the real part, the “current density” Im K is pro-portional to minus the gradient in energy density, withhowever a very small “fictitious” diffusion constant D I = c /k associated with Im K , and a correction from J calculated in the next section.Even if J and J do not carry current themselves,they cannot be ignored because the Bethe-Salpeter equa-tion (39) couples in principle all J i through scattering.From Eq. (39) we can identify four different contribu-tions to J ( p , q ), written as J ( p , q ) = J D + J δ Σ + J δG + J S (49)In this expression, the Drude diffusion current in Eq. (40)has been further split up into the first two terms above.The first is given by, J D ( p , q ) = G ( p ) · L ( p , q ) · G ∗ ( p ) − G ( p ) · L ( p , q ) · G ( p )(50)The second term is generated by the explicit dependenceof the self-energy on wave number, J δ Σ ( p , q ) = − Re G ( p ) · ( q · ∂ p ) Σ ( p ) · G ( p ) (51)with the convention that Re A = ( A + A ∗ ) /
2. The fi-nal two terms J δG and J S are defined as the two lastscattering terms involving U pp (cid:48) in Eq. (39).The mode J implies a new mechanism of long-rangediffusion mixing near and far fields. One peculiarity is the direction of the Poynting vector associated with thediffuse mode expressed by Eq. (37). The mode J of thepure transverse field generates a Poynting vector whosecomponent along the gradient varies as cos θ in phasespace, with θ the angle between wave vector p and gra-dient vector q , and is thus largest along the gradientvector. For the mode J this component varies as sin θ ,which is largest orthogonal to the gradient vector.In the following subsections III A 1–III A 3 we discussthese 4 contributions to J separately, and show that thescattering from two electric dipoles generates all fourchannels in Eq. (46). The results are summarized in sub-section III A 5 and in Table (I).
1. Drude current tensor
The Drude current tensor J D ( p , q ) defined Eq. (50) isindependent of the collision operator U pp (cid:48) and is there-fore the easiest to calculate. We will split J D ( p , q ) fur-ther up into a pure transverse part and an interferenceterm and write J D ( p , q ) = J DT T ( p , q ) + J DT L ( p , q ) (52)The first part stems from purely transverse propagationand contributes only to the J -channel in Eq. (46). Thesecond term is produced by a mixture of longitudinal andtransverse propagation and contributes to the channels J , J and J . Since p · L ( p , q ) · p = 0, the Drude currenttensor features no purely longitudinal mode J DLL ( p , q ),proportional to | G L ( p ) | .The transverse Green’s function G T ( p ) is given by thesecond term in Eq. (2). It follows J DT T ( p , q ) = 2( p · q ) ∆ p (cid:0) | G T ( p ) | − Re G T ( p ) (cid:1) = 4( p · q ) ∆ p Im G T ( p ) (53)This function is heavily peaked near the frequency shellof the effective medium. We can ignore any p -dependencein Σ T ( p ) and approximate it by Σ T ( p = k ). For G T =( K T − p ) − and K T = k e + i/ (cid:96) a complex wave vectorindependent of p , we can use, (cid:80) p p Im G T ( p ) (cid:80) p − Im G T ( p ) = k e (cid:96) In terms of the density of transverse states (DOTS) thisproduces the classical Drude diffusion constant in the J -channel, D D ( k ) = 13 × (cid:18) c k e k N T ( k ) N ( k ) (cid:19) × (cid:96) ( k ) ≡ v E (cid:96) (54)It is customary to write k e /k = c /v p in terms of thephase velocity v p . The ratio N T /N = N T / ( N L + N T )is a factor that can be very small near the resonance ω . We recall that for our electric dipole scatterers allstored energy resides in the longitudinal field. In the3Drude approximation for the transverse field we recoverthe familiar picture of light diffusion, with the extinctionlength as the mean free path, and v E as energy transportvelocity [6].The perturbation expansion in q is valid for the trans-verse waves as long as 2 pq | Im G T ( p ) | < | Im G T ( p ) | .This is most stringent near the frequency shell p = k e where the spectral function | Im G T ( p ) | is maximal andnot stringent at all for large momenta. This gives q < | Im Σ T ( k e ) | / k e = 1 / (cid:96) . This could have been an intu-itive estimate.The diffusion tensor J DT L is given by J DT L ( p , q ) = 2Im G T Im G L (2ˆ p ˆ p ( p · q ) − pq − qp )+ i Im [ ¯ G L G T ] ( pq − qp ) (55)with contributions to the channels J , J and J inEq. (46). We focus first on the Poynting vector for whichonly the channel J ( p ) is relevant. Inserting the first terminto Eq. (44) gives πN ( k ) D D ( k ) = (cid:88) p Im G T ( p )Im G L ( p ) [ p − ( p · ˆ q ) ](56)This diffusion is clearly determined by the overlap oftransverse and longitudinal modes in phase space. Be-cause the longitudinal spectral function is essentially in-dependent of p , this overlap is significant and the integraleven diverges as (cid:80) p /p . We can extract and regularizeit as earlier by Q k / π , πN ( k ) D D ( k ) = Q π (ImΣ ISA ) k (57)This singular term is proportional to the square of thedensity of the dipoles and will later be seen to cancel. Asa result, the leading current tensor is, J DT L = 2 πk(cid:96) δ ( k − p ) 1 p (2ˆ p ˆ p ( p · q ) − pq − qp )+ i Im [ ¯ G L G T ] ( pq − qp ) (58)The use of the Dirac distribution implies here implicitlythat a typical Kubo integral of the kind (cid:80) p pJ ( p , q ) con-verges for large p , with no need for regularization. Thisexpression will turn out to be leading for J and domi-nating for J . The Drude diffusion constant associatedwith the mixture of transverse and longitudinal waves isthus given by D D = 13 πN ( k ) (cid:88) p πk(cid:96) δ ( k − p ) ≈ v E k (cid:96) (59)This diffusion constant can be considered as the ISA ofelectromagnetic diffusion in the J -channel. Its value is positive and, apart from the universal pre-factor v E indiffusion, depends linearly on the density of the dipoles.In Ref. [4] one finds a correction induced by dipole-dipolecoupling that can be written as ∆ D = v E F ( δ ) /k (cid:96) , with the function F varying over the resonance. Like thediffusion found in Eq. (59), it is positive and proportionalto v E , but unlike Eq. (59) it scales as n . The interfer-ence of longitudinal and transverse waves is excluded inRef. [4] which explains why this leading term (59) is notfound.For the hydrodynamic expansion made in Eq. (37)to hold for the transport channel J , we must have | Im G T ( p ) | > | J DT L ( p ) | pq/
2, or equivalently, pq < / | Im G L | ≈ | k / | Im Σ | . Since transverse waves alreadyimpose q < /(cid:96) we conclude that p < k (cid:96) . This becomesrestrictive once k(cid:96) approaches unity.It is straightforward to obtain the Drude approxima-tion for the fictitious diffusion constant in Eq. (48) asso-ciated with Im K , and which was seen to be governed by J . Since J = Im ( ¯ G L G T ) we can write1 πN Im (cid:88) p p ( ¯ G L G T ) = 1 πN Im (cid:88) p (cid:20) p ¯ G L − p + z ¯ z G T (cid:21) ≈ − c k N L + N T N L + N T where we used the expression (18) of the DOS split upin its longitudinal part N L and its transverse part N T .With v E = c N T / ( N L + N T ) we find from Eq. (48), D DI = 13 v E k − (60)The singular longitudinal DOS, proportional to Q , can-cels. In the Drude approximation the fictitious diffusionconstant D I of the mode J is a factor k(cid:96) larger thanthe diffusion constant D of the channel J , and a fac-tor k(cid:96) smaller than the transverse diffusion of mode J .This suggests that they all become of the same order near k(cid:96) = 1.
2. Self-energy dependent on wave number
Any dependence on p of the self-energy contributesto the diffusion current via the term J δ Σ derived inEq. (51). For electric dipoles such dependence on wavenumber comes in via the boomerang diagrams discussedin Eq. (11) with the subtle local field correction at largemomenta derived in Eq. (13), on which we shall focus. Ifwe insert this term into Eq. (51) we find an interferenceterm between longitudinal and transverse propagation, J δ Σ ( p , q ) = − Re (cid:18) n t k G T ( p ) G L ( p ) (cid:19) p × (2ˆ p ˆ p (ˆ p · q ) − ˆ pq − q ˆ p )Its contribution to the Poynting vector in the J -channeldiverges again as (cid:80) p /p . Restricting to large wavevectors, πN ( k ) D δ Σ2 ( k ) = Q π Re (cid:0) n t (cid:1) k (61)4The remainder of Σ ( p ) in Eq. (11) provides contribu-tions to J , J and J , and is proportional to n once thedivergency has been removed. Some formula manipula-tion gives the following closed expression for the diffusionconstant caused by the dependence on wave number ofthe self-energy of two electric dipoles, πN ( k ) D δ Σ ( k ) = 14 n Re Tr (cid:90) d r ( r · ˆ q ) (cid:18) t G − t G − t G · G ,T (cid:19) (62)This expression is free from any singularity but is beyondthe scope of this work, being a factor 1 /k(cid:96) smaller thanwhat was found in Eq. (59) for the J -channel, and evena factor 1 / ( k(cid:96) ) smaller than the leading contribution inthe J -channel.
3. Scattering diffusion current tensor
The scattering diffusion current tensor J δG ( p , q ) isgiven by the second term in Eq. (39). It vanishes forany isotropic scattering in U pp (cid:48) , among which (here)single scattering. Among the different scattering eventsgenerated by two electric dipoles, only the most-crosseddiagrams and the forward-crossed diagrams induce ananisotropy in scattering. They are given by U MC pp (cid:48) = n | t | (cid:90) d r t G − t G (cid:18) t G − t G (cid:19) ∗ e i ( p + p (cid:48) ) · r (63)and U FC pp (cid:48) = n | t | (cid:90) d r (cid:20) − t G (cid:18) − t G (cid:19) ∗ − (cid:21) e i ( p − p (cid:48) ) · r (64)The most-crossed diagrams generate a contribution to J δG ( p , q ) of the type J leading to a diffusion constantfree from any singularity at large p , and of the same orderas was found in Eq. (62). We will ignore them for thesame reason and focus on the forward-crossed diagrams.We can write (cid:88) p (cid:48) U FC pp (cid:48) · δ q G ( p (cid:48) ) = n | t | (cid:90) d r ( i q · r ) (cid:20) − t G · Re G · (cid:18) − t G (cid:19) ∗ − Re G (cid:21) e i p · r This integral is regular for all p , but does not decay fastenough with p to prevent singularities in the channels J and J . To see this, the factor between brackets is writtenas F = F ( r ) + F ( r )ˆ r ˆ r . The space integral above canbe done to get, J δG ( p , q ) = n | t | G ( p ) G ( p ) ∗ · (ˆ p · q ) f ( p ) + ˆ p ˆ p f ( p ) + (ˆ pq + q ˆ p ) f ( p ) with 3 known functions related to F i ( r ). The first termwith f ( p ) is part of J , and constitutes a high-ordercorrection to transverse diffusion, of no interest here. Thelongitudinal term with f produces no Poynting vector.We concentrate on the term with f , given by pf ( p ) = − (cid:90) d r F ( r ) j ( pr )This integral is finite for all p but does not decay with p since lim p →∞ pf ( p ) = lim p →∞ − p (cid:90) d y F ( y/p ) j ( y )= − πk (cid:90) d y j ( y ) y = − k We have used that for small r , F ( r ) = − δ ( r ) / k − (1 − r ˆ r ) / πk r . The local contact term does not contribute.The diffusion constant in the J -channel is given by, πN ( k ) D δG ( k ) = − n | t | (cid:88) p pf ( p )Re G ,L ( p ) ¯ G ,T ( p )(65)This equation thus suffers from a divergence. Upon split-ting it off and regularizing (cid:80) p /p = Q k/ π we find πN ( k ) D δG ( k ) = − Q n | t | πk + O ( n ) (66)This is the third diverging term that will cancel againstthe two already found earlier. The term proportional to f ( p ) also produces a contribution to the J -channel, J δG ( p ) = − n | t | k f ( p ) p Im G T ( p ) (67)This function decays rapidly as 1 /p for large p . It is eas-ily checked that the integral (cid:80) p pJ ( p ) is not singular atlarge p and produces a correction of order n in Eq. (47)that will not be further discussed.
4. Weak localization
The last term in Eq. (39), defined as J S ( p , q ), mixes inprinciple all four transport mechanisms J i . For our modelof electric dipoles, the ISA makes no contribution sinceisotropic, but the diagrams (63) and (64) do. The leadingorder is obtained by inserting on the right hand the Drudeexpression for the transverse diffusion current tensor J T T found in Eq. (53). Since this current is strongly peakednear p = k we can approximate J T T ( p , q ) = 2 π(cid:96) (ˆ p · q ) ∆ p δ ( k − p ) so that, J S ( p , q ) = k(cid:96) π G ( p ) G ( p ) ∗ · (cid:90) d ˆ k π U pk · ∆ k (ˆ k · q ) (68)Only the angle-dependent scattering U MC and U F C sur-vive this integral. For convenience we can summarize5Eqs. (63) and (64) by U pp (cid:48) = (cid:90) d r (cid:104) U MC ( r ) e i ( p + p (cid:48) ) · r + U F C ( r ) e i ( p − p (cid:48) ) · r (cid:105) (69)The angular integral over ˆ k can be performed to get J S ( p , q ) = k(cid:96) πi G ( p ) G ( p ) ∗ · (cid:90) d r e i p · r · (cid:2) U MC ( r ) − U F C ( r ) (cid:3) · (cid:20) j ( kr ) kr (ˆ rq + q ˆ r ) − (ˆ r · q ) (cid:18) j ( kr ) − j ( kr ) kr + j ( kr )ˆ r ˆ r (cid:19)(cid:21) The integrand of this equation for J S ( p , q ) involvesthe difference U = U MC − U F C between most-crossedand forward-crossed diagrams. They both contain sub-radiant poles (where t G ≈ J S ( p , q ) = J S ( p )(ˆ p · q ) + J S ( p )ˆ p ˆ p (ˆ p · q )+ J S ( p )(ˆ pq + q ˆ p ) + J S ( p ) i (ˆ pq − q ˆ p )We will show that the mode J S exhibits the standardweak localization correction, of relative order 1 /k(cid:96) andnegative in diffusion constant. Also the mode J S is sub-ject to a weak localization correction, of order 1 / ( k(cid:96) ) and positive , showing that not all modes are affected sim-ilarly by interference.We first focus on J S . Contrary to U F C , U MC asso-ciated with two dipoles induces a singular angular de-pendence of the kind 1 / | p + p (cid:48) | , and therefore dom-inates J S . The space integral is dominated by large r so we insert U MC = (6 π/(cid:96) ) C ( r ) C ( r ) ∗ with C ≈− ∆ r (exp( ikr ) / πr ). The angular integral over ˆ r can bedone. The end result is written as J MC ( p ) = − k(cid:96) | G T ( p ) | × (cid:90) ∞ dr (cid:18) j ( kr ) − j ( kr ) (cid:19) (cid:18) j ( pr ) − j ( pr ) (cid:19) = − π
20 1 (cid:96) | G T ( p ) | (cid:18) k p + 97 k p (cid:19) The last equality holds only for p ≥ k ; J MC ( p ) decaysrapidly as 1 /p and has most of it weight near p = k .The weak localization correction can be obtained fromEq. (44),∆ D W L = 14 πN ( k ) Tr (cid:88) p L ( p , ˆ q ) · ∆ p (ˆ p · ˆ q ) × − π k δ ( k − p ) = − v E k π
280 (70)or equivalently ∆ D W L /D T = − . /k(cid:96) . The numeri-cal factor is actually larger than the leading one ( π/ . J channel. If we extrapo-late to small values for k(cid:96) , we conclude that the diffusionin the J -channel compensates the first weak localizationcorrection in the J -channel for k(cid:96) < . J S is more complicated. It is instructivesplit the Green’s function up into G ( r ) ∼ P ( r ) ∆ r + Q ( r )ˆ r ˆ r , and to express the tensor U ( r ) = U MC − U F C as, U ( r ) = U T T ( r ) ∆ r ∆ r + U LL ˆ r ˆ r ˆ r ˆ r +Re U T L (ˆ r ˆ r∆ r + ∆ r ˆ r ˆ r ) + i Im U T L ( ∆ r ˆ r ˆ r − ˆ r ˆ r∆ r )This corresponds to 4 different scattering events involv-ing two dipoles at distance r with the electric field vectoreither along or perpendicular to r , as well as their in-terferences. With the angular integral of ˆ r performedanalytically, they give each the following contribution to J S , J T T ( p ) = 2 (cid:96)k Re G T ( p ) (cid:90) ∞ dr r U T T ( r ) × (cid:18) j ( kr ) − j ( kr ) kr (cid:19) j ( pr ) prJ LL ( p ) = − (cid:96)k Re G T ( p ) (cid:90) dr ∞ r U LL ( r ) j ( kr ) kr j ( pr ) prJ T L ( p ) = 2 (cid:96)k Re G T ( p ) (cid:90) ∞ dr r Re U T L ( r ) × (cid:18) j ( pr ) − j ( pr ) pr (cid:19) j ( kr ) kr and finally J T L ( p ) = 2 (cid:96)k Im G T ( p ) (cid:90) ∞ dr r Im U T L ( r ) j ( pr ) j ( kr ) kr The weak localization correction is found by∆ D ( k ) = −
13 1 πN ( k ) (cid:88) p p ( J T T + J LL + J T L + J T L )(71)To perform the integral over the wave vector p we usethat (cid:88) p p ( k + i(cid:15) ) − p j ( pr ) = − ik π h (1)1 ( kr )and (cid:88) p p ( k + i(cid:15) ) − p j ( pr ) pr = − ik πr (cid:20) h (1)2 ( kr ) + 3 i ( kr ) (cid:21) Consequently, only the radial integrals (cid:82) dr remain tobe done numerically. The weak-localization correction∆ D ( k ) is proportional to the density of the electricdipoles.Figure 7 shows the total diffusion constant D D + ∆ D in the J channel around the resonance frequency, as6 −2−1.5−1−0.5 0 0.5 1 1.5 2 2.5−4 −3 −2 −1 0 1 2 3 4 D / D [ / ( k l ) ] Detuning [line width]
WL(TT)WL(LL)WL(TL)WL(TL2)Drude + WL
Figure 7. The diffusion constant in the J channel (solidline), being the sum of the Drude approximation D D plusthe weak-localization correction ∆ D ( k ) in Eq. (71) from twodipoles, as a function of the detuning δ = ( ω − ω ) /γ . It isnormalized by the diffusion constant in the J channel. The4 weak localization corrections discussed in this section areshown separately as dashed lines. well as the contributions stemming from the 4 individualterms in Eq. (71). The weak localization correction ∆ D is dominated by the purely transverse and longitudinalchannels T T and LL between the two dipoles, who havecompeting signs. For negative detuning the purely lon-gitudinal mode LL dominates, for positive detunings the T T channel dominates and is more than twice as largeas the Drude contribution D D . We note that the weaklocalization correction to the diffusion constant D of thechannel J ( p , q ) is of the same order as the Drude ap-proximation in Eq. (59). The sum of the 4 weak localiza-tion terms and the Drude approximation is strictly pos-itive. At fixed density, positive detuning has the largestdiffusion constants in the J channel. Note that the ra-tio D /D is of same order 1 / ( k(cid:96) ) , but of opposite signcompared to the standard (Cooperon) weak localizationcorrection − / ( k(cid:96) ) . This will be discussed more in de-tail in the next section, for which it will turn out usefulto define a function F ( δ ) = ( k(cid:96) ) D /D .
5. Summary of previous subsections
We have identified four mechanisms in the transportof electromagnetic waves, expressed by the diffusion cur-rent tensor (46). The results have been summarized inTable I. The mechanism described by J ( p ) is the familiarpicture of transverse wave diffusion near the shell p ≈ k and results in the diffusion constant (54). It is inverselyproportional to the density of the electric dipoles andcontains an energy velocity that can be small since theimpenetrable electric dipole scatterers contain temporar-ily stored, longitudinal energy. The mechanism associ-ated with J is caused by interference of longitudinal and Drude J D J δ Σ J δ G WL J S J + (cid:96) /k (cid:96) /k (cid:96) − . /k (54) (62) NC (70) J +1 /k (cid:96) − Q ( δ ) / k (cid:96) − Q ( δ ) / k (cid:96) + F ( δ ) /k (cid:96) + Q ( δ ) /k (cid:96) (59) (61) (66) (71) Fig. 7 J +1 /k /k (cid:96) /k (cid:96) /k (cid:96) (60) (62) (67) NCTable I. Contributions to the transport mean free path forvarious transport channels J i and the 3 different diagram-matic classes identified in the Bethe-Salpeter equation (49).When an explicit sign is found, it is indicated. Most valuesdepend also explicitly on detuning, not indicated if not cal-culated. Numbers refer to the corresponding equations. Thechannel J does not contribute to transport mean free path, J only contributes to the transport mean free path associatedwith the imaginary part of the Poynting vector. NC standsfor “not calculated”, WL for “weak localization”. The termsproportional to Q ( δ ) are regularized singularities dependingon detuning δ that cancel in the transport mean free path. transverse fields, necessary condition to carry a Poyntingvector. The leading term (59), linear in the dipole den-sity, comes from the Drude approximation. Upon consid-ering all scattering events involving two dipoles, we havebeen able to identify three singular terms. After regular-ization, they are expressed by Eqs. (57), (61) and (66)and proportional to the large quality factor Q and thedensity squared. They add up to πN ( k )∆ D ( k ) = Q πk n (cid:18) (Im t ) + 12 Re t − | t | (cid:19) = 0 (72)This explicit cancelation in the J -channel is very im-portant and not entirely obvious since the 3 terms stemfrom entirely different parts in transport theory (Drudediffusion, Lorentz local field and enhanced forward scat-tering). Without cancelation they would have given anelectromagnetic conductivity Q /k(cid:96) , and not small at allwith respect to the traditional transverse conductivity, oforder k (cid:96) since Q is large for an atomic oscillator. Theircancelation also supports the general renormalizability ofelectromagnetic transport theory with point-like dipoles.It is highly plausible that this cancelation happens in allorders of perturbation theory, but this is currently im-possible to prove in general. We will use this hypothesisin the next session.The J channel, in the leading order modified by theweak localization from 2 electric dipoles, exhibits a pos-itive diffusion constant, linear in the dipole density. Al-though usually small compared to standard transversediffusion, it must be realized that this diffusion stemsfrom a sofar unexplored mechanism for electromagneticwave diffusion, involving the interference of longitudinaland transverse waves. In this transport channel, thefirst weak localization correction induced by two elec-tric dipoles is actually of the same order as the Drude7value and again positive , showing that in the channel J interferences behave differently in comparison to the tra-ditional transverse channel. IV. RADIATIVE FORCE DENSITY
A well-known relation exists between diffuse flow andradiative forces. In radiative transfer, the energy flux isdriven by the spatial gradient of the total energy density,and automatically carries momentum. In the presence ofan induced polarization density P , the electromagneticforce density f is caused by the Lorentz force acting onthe induced Coulomb charge density ρ ( r ) = − ∇ · P andon the induced current density j = ∂ t P . Maxwell’s equa-tions allow the formulation of a momentum conservationlaw that is very generally valid. It takes the form (beforecycle averaging) [29, 30], ∂ t G + f = ∇ · T (73)with G = ( E × B ) / πc the electromagnetic momentumdensity, and T the momentum stress tensor, T = 14 π (cid:20) EE + BB − (cid:0) E + B (cid:1) + X (cid:21) (74)The tensor X is related to internal angular momentuminside the particle that we shall ignore here.In the regime of multiple scattering, and after cy-cle averaging, Eq. (37) expresses that (cid:104) E i ( r ) ¯ E j ( r ) (cid:105) = (cid:104)| E ( r ) | (cid:105) δ ij , and idem for the magnetic field. Thestress-tensor T is thus diagonal on average, meaningthat the i th component of the electromagnetic momen-tum only flows in the direction i . For stationary flow, wethus obtain (cid:104) f ( r ) (cid:105) = − ∇ (cid:28) | E ( r ) | + | B ( r ) | π (cid:29) (75)For a medium filled with impenetrable electric dipoleswe have shown in Eq. (18) that | E | / π is the to-tal electric energy density having both longitudinaland transverse components, and equal to the mag-netic energy density. In the diffusion approximation,we write the averaged Poynting vector as (cid:104) K (cid:105) = − D ∇ (cid:10)(cid:2) | E ( r ) | + | B ( r ) | (cid:3) / π (cid:11) . This leads to a simplerelation (cid:104) f (cid:105) = 13 D (cid:104) K (cid:105) (76)between Poynting vector and radiative force density. Inthe ISA, D = v E (cid:96)/ (1 − (cid:104) cos θ (cid:105) ), and this reduces to thealmost intuitive expression (cid:104) f (cid:105) = nσ (1 − (cid:104) cos θ (cid:105) ) (cid:104) K (cid:105) /v E involving the product of particle density and pressurecross-section of one scatterer. The second factor ac-counts for transfer of momentum from the light to a singlescatterer, and of course for independent electric dipoles (cid:104) cos θ (cid:105) = 0. The factor 1 /v E is less intuitive in this model. Forone isolated scatterer this would clearly be 1 /c , sincefor a plane wave with arbitrary direction in vacuum, mo-mentum current density and energy current density (thePoynting vector) differ by a factor 1 /c . In a mediumfilled with resonant dipoles, stocked, longitudinal energycontributes to the momentum current density T but notto the energy current density K . Put otherwise, scatter-ing of a transverse state with wave number p ≈ k to alongitudinal mode with large wave number induces a sig-nificant recoil, but does not generate an energy current.For the medium filled with dipoles, the ratio f /K thusachieves a factor ( N L + N T ) /N T c ≈ /v E . V. NO ANDERSON LOCALIZATION OFLIGHT?
In the following we will make a first attempt to includethe 4 transport mechanisms, introduced in the previoussection, into the self-consistent transport theory for lo-calization of light. This theory is celebrated by somefor its surprisingly simple description of the transitionfrom long-range diffusion to localization. Others criti-cize the theory for its oversimplified nature, neglectingmany scattering events in the collision operator U pp (cid:48) ( q )introduced in Eq. (33). The self-consistent theory pre-dicts the Ioffe-Regel criterion k e (cid:96) ≈ U pp (cid:48) ( q ). By reciprocity, the most-crossed diagrams also contain a hydrodynamic pole, fea-turing the same diffusion constant. This immediatelyturns the calculation of D into a self-consistent problembecause the most-crossed diagrams, modify the diffusioncurrent J ( p , q ) in Eq. (39). We recall that in the caseof electromagnetic waves the diffusion current is a ten-sor with 4 independent parts. No Anderson localizationwas seen to occur in recent numerical simulations withelectric dipoles [5]. The intention of this section is todiscover what exactly breaks down in this theory whentaking into account longitudinal waves.We here summarize the various approximations made,which are basically equivalent to the ones made in previ-ous works, even if often adopted implicitly [7, 31, 32]. • The most-crossed diagrams, involving scatteringevents associated with many dipoles, are the onlyangle-dependent scattering events that influencethe diffusion current tensor J ( p , q ) when going be-yond the Drude approach. The existence of other8 diagrams is only acknowledged implicitly to guar-antee flux conservation. Weak localization effectscaused by low-order scattering events, such as thosedescribed by Eqs. (70) and (71) are not includedeither, although this could be done without dra-matic changes in the theory. The implicit exis-tence of other diagrams is also necessary to jus-tify the cancelation of UV-singularities in trans-port theory. For scattering events involving onlytwo electric dipoles, UV-divergencies were seen tocancel explicitly earlier in this work, but no generaldemonstration is known. • The diffuse regime of the most-crossed diagrams,only valid on spatial scales well beyond the meanfree path, is assumed to hold on scales up to themean free path. On this scale we may expect thediffusion kernel to be of the type D ( q ) q whichis disregarded in the standard version of the self-consistent theory. • The electromagnetic self-energies Σ
T/L ( k, p ) are as-sumed not to depend on p . In particular this meansthat Σ T ( k ) = Σ L ( k ). This is definitely an approxi-mation, even for point-like dipoles, that needs morestudy, but in general, such wave number depen-dence is not believed to be essential for Andersonlocalization.The contribution of the most-crossed diagrams to thescattering vertex U pp (cid:48) ( q ) can be obtained from the re-ducible vertex R pp (cid:48) ( q ) introduced in Eq. (31) by re-moving the four external Dyson propagators, and time-reversing the bottom line. This gives U MC pp (cid:48) ; ij | kl ( q ) = ˜ d il ( f + q , Q ) ˜ d kj ( − f + q , Q ) + O ( Q ) πN D Q (77)with the notation Q = p + p (cid:48) and f = ( p − p (cid:48) ) /
2. In thisexpression the tensor ˜ d ( p , Q ) is the diffuse eigenfunctiondefined in Eq. (37) stripped from the 4 external lines inFig. 6 (transforming − Im G + i J / − Im Σ + i j / j again a Hermitian bilinear form). This leads to˜ d ( ± f + q , Q ) = − Im Σ( ± f + q ) + j ( ± f + q , Q ).A generalized Ward identity,( q · ∂ p )Re Σ ( p ) = (cid:88) p (cid:48) U pp (cid:48) (0) · ( q · ∂ p (cid:48) )Re G ( p (cid:48) )+ (cid:88) p (cid:48) δ q U pp (cid:48) ( q ) · Im G ( p (cid:48) ) (78)can be used to eliminate the second term in Eq. (39),which then transforms into J ( p , q ) = J D ( p , q )+ G ( p ) · ( q · ∂ p )Re Σ ( p ) · G ∗ ( p ) − G ( p ) G ∗ ( p ) · (cid:88) p (cid:48) δ q U pp (cid:48) ( q ) · Im G ( p (cid:48) )+ G ( p ) G ∗ ( p ) · (cid:88) p (cid:48) U pp (cid:48) · J ( p (cid:48) , q ) (79) The first and second terms cannot depend on diffusionconstant. Because U MC depends on both diffusion con-stant D and the entire diffusion tensor J ( p , q ), the self-consistent theory would, in its most advanced version, bea non-linear integral equation for the second-rank tensor J . In the following we apply the approximations speci-fied above. The above hydrodynamic limit of U MC isassumed valid when | p + p (cid:48) | (cid:28) /(cid:96) . In the standard ap-proach of the self-consistent theory one focusses on itsdiffuse pole near p ≈ − p (cid:48) , and neglects all other de-pendence on p (cid:48) . Secondly, wave number dependence ofthe self-energy is ignored. In that case the self-consistentproblem simplifies to the following equation, J ( p , q ) ≈ J D ( p , q ) + G ( p ) G ( p ) ∗ · (cid:88) | Q | In most applications of the self-consistent theory forlocalization of light one ignores polarization and focusseson the transverse channel J ( p ) and, not unrelated, as-sumes this channel to be governed by excitations nearthe frequency shell p ≈ k e of the effective medium wheretheir DOS is largest. In this approximation weak local-ization of light becomes essentially equivalent to the oneof scalar waves. As a matter of fact, this approximationapplies to localization of elastic waves with all polariza-tions modes propagating with the same velocity every-where. We refer to the work of Zhang and Sheng [33]where the self-consistent theory for localization of scalarwaves is derived and discussed in great detail.9We will first neglect the weak localization found inEq. (70) associated with two dipoles and incorporate it inthe next section when dealing with the mode T L . Uponputting J ( p , q ) = J ( p )(ˆ p · q ) ∆ p into Eq. (80), and byassuming that Im Σ( p ) is independent of p , the explicitsolution is just J ( k, p ) = J D ( k, p ) (cid:104) σ c σ A ( k, p ) (cid:105) − (82)with the dimensionless function A ( k, p ) = | G T ( k, p ) | × Im Σ T ( k, p ), and a critical conductivity defined as σ c ≡ (cid:80) Q /Q = q m / π . Note that A ( k, p ) ≤ p . Equation (82) says that the amount ofweak localization varies in phase space, and is maximalat the frequency shell p = k e of the transverse waves,and small when p (cid:29) k . From Eq. (81) we obtain a closedequation for σ , σ ( k ) = 13 (cid:88) p p J D ( k, p )1 + ( σ c /σ ) A ( k, p ) (83)The Kubo formula attributes a large weight to large p , nevertheless the integral converges for all σ > 0. If σ > σ c large wave vectors p are not relevant in thedenominator since J D ( k, p ) = 4(Im G T ( k, p )) decaysrapidly with p . The integral is dominated by p nearthe frequency shell p ≈ k e so that σ ( k ) ≈ σ D ( k )1 + ( σ c /σ ) ⇒ σ ( k ) = σ D ( k ) (cid:18) − σ c σ D (cid:19) = σ D ( k ) (cid:18) − π k e (cid:96) ) (cid:19) (84)where q m = 1 /(cid:96) has been chosen. This result, whenextrapolated, locates the mobility edge at k e (cid:96) = 0 . σ < σ c however, the p -dependence of the denomi-nator shifts the integral over p to larger values for p . Atthe mobility edge σ = 0 and σ c = 43 (cid:88) p p | G T ( k, p ) | (85)This involves an integral that diverges as (cid:80) p /p . It canbe verified that this divergence is not an artifact of theapproximation made in Eq. (80). This kind of divergenceis absent in standard approaches of the self-consistenttheory [7, 31] because of the a priori assumptions for J ( p , q ) to be “strongly peaked near the frequency shell” p ≈ k e . This is clearly the case for J D ( p , q ), but notnecessarily true after solving Eq. (80). Since the mo-bility edge in scalar wave scattering by point dipoles isobserved in numerical simulations [34] near k e (cid:96) ≈ 1, thisdivergence must be an artifact and should be eliminated.We could subtract the singularity (cid:80) p /p by hand, as-suming it cancels against other terms that have been ig- nored, to get σ c = 43 (cid:88) p (cid:18) p | G T ( p | − p (cid:19) = k e (cid:96) π (cid:18) − k e (cid:96) ) (cid:19) (86)This would locate the mobility edge at k e (cid:96) = 0 . 866 withthe choice q m = 1 /(cid:96) . This is close to the extrapolatedvalue above, and we could argue that the extrapolation inEq. (84) is satisfactory up to the mobility edge and con-sistent with both previous theory [7, 33] and numericalsimulations [34]. It is nevertheless tempting to speculatethat this divergence highlights a true breakdown of theself-consistent theory and that a more rigorous regular so-lution may actually exhibit a critical exponent differentfrom one, the value predicted by the extrapolation (84). B. Inclusion of longitudinal modes In this section we give a simplified description of howthe self-consistent theory is extended when the other 3diffusion modes are included. Let us start with Eq. (46)and write the diffusion current tensor as J ij ( p , q ) = (cid:88) n =0 J n ( p ) χ nij ( p , q ) (87)Let us set U MCij ; kl = ( U/σ ) δ kj δ il with U = (Im Σ) σ c (withdimension 1 /m ) and σ = πN D the conductivity (withdimension 1 /m ). We can check that, G ni G ∗ jm δ kj δ il χ kl = | G T | χ nm G ni G ∗ jm δ kj δ il χ kl = | G L | χ nm G ni G ∗ jm δ kj δ il χ kl = R ( χ nm − χ nm ) + Iχ nm G ni G ∗ jm δ kj δ il χ kl = − Rχ nm + I ( χ nm − χ nm )where we abbreviated R ( p ) = Re G L ¯ G T and I ( p ) =Im G L ¯ G T . This gives the following self-consistent set ofequations J ( p ) = J D ( p ) − (cid:18) Uσ | G T ( p ) | + 0 . k e (cid:96) (cid:19) J ( p ) (cid:18) | G L | Uσ (cid:19) J ( p ) = J D ( p ) + 2 R ( p ) Uσ J ( p )+2 I ( p ) Uσ J ( p ) (cid:18) R ( p ) Uσ (cid:19) J ( p ) + I ( p ) Uσ J ( p ) = J D ( p ) I ( p ) Uσ J ( p ) + (cid:18) − R ( p ) Uσ (cid:19) J ( p ) = J D ( p ) (88)The equation for the transverse mode J discussed in theprevious section is not altered and decouples from the0others. We have added the weak-localization contribu-tion caused by 2 dipoles found in Eq. (70), since it is notcovered by the diffusion approximation, and assumed itenters just as a number in the equation for J ( p ). This isa clear oversimplification but has no huge consequencesfor what follows. The purely longitudinal diffusion cur-rent J is known once the others are known, but is notrelevant for Poynting vector and can likewise be ignored.The modes J and J however, couple and the solutionfor J is J ( p ) = J D ( p ) + ( U/σ ) C ( p )1 − U | G L ( p ) G T ( p ) | /σ (89)We recall from Eq. (55) that J D ( p ) = − G L Im G T < J D ( p ) = − I ( p ). Thus, with K = k e + i/ (cid:96) thecomplex pole of G T ( p ), the function C ( p ) is given by C ( p ) = I ( p ) − J D ( p ) R ( p )= (cid:18) k e (cid:96) (cid:19) | G ( p ) | + | K | | G ( p ) | | K | (90)which is strictly positive.Before calculating diffusion constant we first discussthese results. Since the wave number integral of J ( p )contributes to the diffusion constant via Eq. (81), its de-nominator cannot possess any non-integrable singularity.This implies that σ ( k ) > U | G L ( p ) G T ( p ) | (91)to be valid for all p . This inequality excludes de facto amobility edge. It is most stringent near the transversefrequency shell p ≈ k e (more precisely p = Re K = k e − / (cid:96) , positive as long as k e (cid:96) > / 2) where G T =1 / ( − i Im Σ). Furthermore, since we neglect p -dependencein self-energies we set | G L | = 1 / | K | and neglect the factthat near the transverse shell the complex wave num-bers of longitudinal and transverse modes are not neces-sarily equal. Recalling that U = (Im Σ) σ c and setting q m = q/(cid:96) , with q of order unity, the minimal possibleelectromagnetic conductivity is given by σ ( k ) > c ( k e (cid:96) ) σ D ( k ) (92)with c ( x ) = (3 q/πx )( x + ) − for k e (cid:96) > / 2. Equiv-alently, if the transport mean free path (cid:96) ∗ is defined asusual via σ = k e (cid:96) ∗ / π [27], then k e (cid:96) ∗ > qπ k e (cid:96) ) + (93)for k e (cid:96) > / 2. For k e (cid:96) < / p = 0 and we find k e (cid:96) ∗ > qπ ( k e (cid:96) ) [( k e (cid:96) ) + ] (94)The very existence of this minimum conductivity forvector waves is determined by scattering properties of t r an s po r t m ean f r ee pa t h k l * Ioffe−Regel parameter kl full solutionmininum thresholdconventional solutionfictitious Figure 8. The self-consistent solution for the electromagnetictransport mean free path (cid:96) ∗ defined by σ = k e (cid:96) ∗ / π . Shownare the values for k e (cid:96) ∗ for the full solution in this section,the conventional picture described by Eq. (81) with only thetransverse mode J considered, with a mobility edge predictedaround k e (cid:96) ≈ 1, the lower threshold imposed by the existenceof the diffusion modes J and J , as well as k e (cid:96) ∗ associatedwith the fictitious conductivity and J . We used a cut-off q m = 1 /(cid:96) . longitudinal and transverse waves near the frequencyshell and not by large wave numbers p that are subjectto uncertain regularization procedures. It neverthelessrelies on our choice for q , and the approximation that K L ( p ) = K T ( p ) = K . The above lower bound becomesstringent for k e (cid:96) ≈ p < k , so that setting K L ( p ) = K T ( p ) = K may not bea bad approximation, knowing that for p (cid:28) k it is valid(see for instance the p -dependent self-energies in Fig. 1).If q = 1, we find for k e (cid:96) = 1, k e (cid:96) ∗ > . 76, and uponentering the evanescent regime k e (cid:96) = 1 / k e (cid:96) ∗ > . k e (cid:96) = 0 . 35 the maximum value is 2 . σ ≡ σ/σ D = ˆ σ + ˆ σ . Since the mobility no longervanishes, the transverse diffusion mode J , which decou-ples from the others, can be given the same treatmentas done in Eq. (84), with the denominator removed andtaken outside at its maximum value. This gives the firstequation for the conductivity of the transverse channel,ˆ σ = 11 + c ( k e (cid:96) ) / ˆ σ + 0 . /k e (cid:96) (95)with c ( x ) = 3 q/πx . We can apply the same proce-dure for the diffusion current J . However, as was seenin Eq. (59) to be the case for the Drude component, theremaining integral suffers from a divergence at large p ,again of the kind (85). The regularization proposed inEq. (86) is not satisfactory here since it changes sign at k e (cid:96) = 0 . 86 and would produce a negative Drude con-ductivity in the T L -channel, arguably not physical. In1section III A 1 we found that for p > k e (cid:96) the diffusiontheory in the J channel breaks down so that the presenttheory is not valid for too large p . We therefore proposea regularization (cid:88) p p | G T ( p ) | → | K | (cid:88) p | G T ( p ) | = | K | (cid:96) π with K = k e + i/ (cid:96) the transverse complex wave num-ber. In real space is G T ( r ) = − exp( iKr ) / πr and thisregularization comes down to (cid:90) d r (cid:12)(cid:12)(cid:12)(cid:12) ∇ (cid:18) exp( iKr ) − πr (cid:19)(cid:12)(cid:12)(cid:12)(cid:12) → | K | (cid:90) d r (cid:12)(cid:12)(cid:12)(cid:12) exp( iKr ) − πr (cid:12)(cid:12)(cid:12)(cid:12) meaning that the regularization only considers the farfield when taking the spatial derivative. In particularthis leads to the Drude diffusion constant of channel J , D D = − πN ( k ) (cid:88) p p J D ( p ) = 4 π v E | K | (cid:96) (cid:88) p p | G T ( p ) | → v E | K | (cid:96) This is a satisfactory, positive extrapolation of the result v E /k (cid:96) obtained in Eq. (81) for low density, and wherethe divergence was seen to cancel explicitly. If we adoptthis regularization, we find in the J -channel,ˆ σ = F ( δ ) c ( k e (cid:96) )1 − ( c ( k e (cid:96) ) / ˆ σ ) (cid:18) − c ( k e (cid:96) )ˆ σ (cid:19) (96)with c ( x ) defined earlier, c ( x ) = ( x + 1 / − and c ( x ) = (3 q/ π )(9 / x / x + 1 / − . We recall that F ( δ ) is the function that describes the explicit depen-dence on detuning of the diffusion constant in the channel J , shown in Fig. 7.Equations (95) and (96) lead to a cubic equation forˆ σ that can be solved analytically. The resulting formulais quite lengthy and we do not present it here. The so-lution for k(cid:96) ∗ = ˆ σ × k(cid:96) is shown in Fig. 8. We haveput F ( δ ) = 1, its role will be discussed later, in whichcase the self-consistent theory has only one parameter,the product k e (cid:96) , as in the scalar case. According toEqs. (95) and (96) the traditional weak localization cor-rection δσ = − c in the transverse channel is partiallycompensated by the positive conductivity δσ = c ofthe J channel, and even exactly when q ≈ 1. This ex-plains why k e (cid:96) ∗ is well in excess of the traditional pre-diction (84), for values as small as k e (cid:96) = 1 . 8, and closeto the Drude value k e (cid:96) of the transverse channel. Theterm containing c > J mode as 1 / ( k e (cid:96) ) but the coupling to J described by c reverses this trend. Around the region k e (cid:96) ≈ J , as described by Eq. (48). The self-consistent solution (cid:19) (cid:20) (cid:21) (cid:22) (cid:23) (cid:24)(cid:19)(cid:17)(cid:19)(cid:19)(cid:17)(cid:24)(cid:20)(cid:17)(cid:19)(cid:20)(cid:17)(cid:24)(cid:21)(cid:17)(cid:19)(cid:21)(cid:17)(cid:24)(cid:22)(cid:17)(cid:19)(cid:22)(cid:17)(cid:24) Figure 9. The ratio of transport and scattering mean freepaths (cid:96) ∗ /(cid:96) as a function of k e (cid:96) compared to the self-consistenttheory for 4 πn/k = 3 . 77 and q = 0 . n fordetunings δ = ( ω − ω ) /γ ∈ [ − , . 5] from the resonance. Thedashed line is the lower bound for (cid:96) ∗ /(cid:96) described by Eqs. (93)and (94), again with q = 0 . is given by J ( p ) = − Im G L ( p ) G T ( p )1 − U | G L ( p ) G T ( p | /σ × (cid:20) Uσ Re G L ( p ) G T ( p ) (cid:21) (97)and upon inserting this into Eq. (48), the same proce-dure as above provides an expression for the “fictitious”conductivityˆ σ I = 1 k e (cid:96) + 11 − ( c ( k e (cid:96) ) / ˆ σ ) (cid:18) d ( k e (cid:96) )ˆ σ − d ( k e (cid:96) )ˆ σ (cid:19) (98)with the functions d ( x ) = (3 q/π ) x − ( x +1 / − and d ( x ) = (3 q/ π ) x − ( x + 1 / − . The transport meanfree path associated with the fictitious diffusion is alsoshown in Fig. 8. For k e (cid:96) ∼ q = 0, Eqs. (95) and (96) simplify to the sum ofthe diffusion constants associated with one or two dipolesin the channels J and J , without any cross-talk,ˆ σ = (cid:96) ∗ (cid:96) = 11 + 0 . /k e (cid:96) + F ( δ )( k e (cid:96) ) + (99)For k e (cid:96) < J channel startsdominating. If we ignore the explicit dependence on δ byputting F ( δ ) = 1, this equation yields ˆ σ < k e (cid:96) > . 73 and below this value starts increasing monotonically.In the same limit of q = 0, we have ˆ σ I = 1 /k e (cid:96) .2 C. Comparison with numerical simulations In Fig. 9 we compare the predictions of the self-consistent theory for electromagnetic waves developedabove to numerical simulations in which we simulate themultiple scattering of light by an ensemble of dipolar res-onant point scatterers. The results of the simulationsallow us to estimate k e , (cid:96) and (cid:96) ∗ . Both the details ofthe simulations and the way in which we interpret theirresults are detailed in Appendix B. We repeat calcula-tions for several atomic number densities n and detun-ings δ = ( ω − ω ) /γ ; the resulting ratios (cid:96) ∗ /(cid:96) are pre-sented in Fig. 9 by circles of different colors as functionsof the Ioffe-Regel parameter k e (cid:96) . The numerical resultsare bounded from below by Eqs. (93) and (94) for theminimum conductivity (dashed line). Equations (93) and(94) impose a sharp rise of the ratio (cid:96) ∗ /(cid:96) at small valuesof k e (cid:96) where one would normally have expected a mobil-ity edge. This rise is well reproduced by the numericalresults.A striking feature of the numerical results is the cleartendency of data to group together along two different“branches”. A careful inspection of Fig. 9 shows thatthe lower branch is composed of data corresponding to δ < δ > 0. This means that -apart from the absence of a localization transition - thereis no one-parameter dependence with k e (cid:96) either. Thedouble-branch structure actually follows from the explicitdependence of the J channel on detuning δ , describedby the factor F ( δ ), which is larger for positive detunings(see Fig. 7). Figure 9 shows the prediction of the self-consistent theory for the ratio (cid:96) ∗ /(cid:96) with the inclusion ofthe function F ( δ ) and with the dimensionless parameter k e (cid:96) calculated from the averaged incident field (see Ap-pendix B) for one fixed dipole density 4 πn/k = 3 . 77 andfor various detunings. Predictions for k e (cid:96) correspond-ing to other densities are not shown since they all ex-hibit the same overall appearance. Despite the fact that,strictly speaking, (cid:96) ∗ /(cid:96) is a function of two independentparameters ( k e (cid:96) and δ or equivalently k e (cid:96) and 4 πn/k ,we see that all results for the quite wide explored den-sity range 4 πn/k = 0 . πn/k = 3 . 77. Theagreement between numerical and analytical results isnot perfect but we believe that it can be further improvedby distinguishing explicitly between transverse and lon-gitudinal complex wave numbers (see Sec. II D), whichare known to be different (see Appendix A, and Figs. 11and 12). This can be done in future work. VI. CONCLUSIONS AND OUTLOOK In this work we have included longitudinal excitationsinto a transport theory for electromagnetic waves prop-agating in a medium with randomly distributed dipo- lar electric scatterers (dipoles). We identify four dif-fuse modes, triggered by the gradient in electromagneticenergy, among which two carry a Poynting vector andcontribute to the diffusion constant. We have developedthis theory by extending the independent scattering ap-proximation (the elementary scattering unit is a singledipole) to include rigourously recurrent scattering fromtwo dipoles. This has led to the following results. 1) Lon-gitudinal and transverse waves of the effective mediumare characterized by different complex wave numbers K L and K T , respectively, and dominate near and far field inscattering. 2) The interference between longitudinal andtransverse waves creates a new diffuse transport chan-nel with a diffusion constant proportional to the numberdensity of dipoles, to be compared to the usual diffusionconstant that is inversely proportional to this density.3) Divergent terms appear at large wave numbers in thediffusion constant, in the longitudinal density of statesand in the collision operator. Many of them cancel, inparticular for the electromagnetic Kubo diffusion con-stant all divergent terms cancel. We postulate that thiscancelation holds in all orders of perturbation theory. 4)When extending the self-consistent theory of localization,with all its usual assumptions, to include the four diffusemodes, we find a minimum conductivity that prevents theonset of Anderson localization of light, as also observedin numerical simulations [5]. 5) The predictions of thedeveloped self-consistent theory are surprisingly close tothe results of independent numerical simulations, includ-ing the explicit dependence of the new transport channelon frequency detuning from the dipolar resonance. Thesefindings demonstrate that, due to the presence of longitu-dinal, non-propagating waves, (weak) localization of lightis fundamentally different from what was believed so far.Early stages of this work were supported by collabora-tions with Yvan Castin, Ad Lagendijk, Nicolas Cherroretand Dominique Delande. We thank Denis Basko for use-ful discussions. Appendix A: Longitudinal and transverse DOS of Melectric dipoles In this Appendix we pose the problem of light scatter-ing from M point electric dipolar scatterers (“dipoles” forshort) in a volume V and derive the DOS for both longi-tudinal and transverse excitations in the thermodynamiclimit M, V → ∞ at constant density M/V = n .The real-space Green’s function G ( r , r (cid:48) ) of M point-like dipoles at positions { r m } , m = 1 , . . . , M , is definedin terms of their collective T -matrix as G ( r , r (cid:48) ) = G ( r − r (cid:48) )+ (cid:88) mm (cid:48) G ( r − r m ) · T mm (cid:48) · G ( r m (cid:48) − r (cid:48) )(A1)If we impose that all dipoles be impenetrable, we must3have G ( r n , r (cid:48) ) = 0 for all r (cid:48) outside the dipoles. Thus,0 = (cid:88) m (cid:48) (cid:34) δ nm (cid:48) + (cid:88) m G ( r n − r m ) · T mm (cid:48) (cid:35) · G ( r m (cid:48) − r (cid:48) )For this to be true for all r (cid:48) , the 3 M × M matrix betweensquare brackets must vanish, { T mm (cid:48) } = − ( { G ( r m − r m (cid:48) ) } ) − (A2)We can split off the singular diagonal elements anduse the fact that the t -matrix of one single dipole is t = − G − (0) and here proportional to the 3 × { T mm (cid:48) } = t ( δ mm (cid:48) − t { G ( r m − r m (cid:48) (cid:54) = 0) } ) − (A3)This matrix is regular as long as the dipoles do not over-lap. For any point source s ( r (cid:48) ) located at r (cid:48) the electricfield anywhere in the medium is given by E ( r ) = G ( r , r (cid:48) ) · s ( r (cid:48) ), and the incident field is E ( r ) = G ( r , r (cid:48) ) · s ( r (cid:48) ). Ifthe source is located in the far field of the M dipoles,and the origin r = 0 is chosen inside the scatteringmedium, we have | r m | (cid:28) | r (cid:48) | , and we can approxi-mate the incident field inside the medium as E ( r ) =[ − exp( ikr ) / πr ] ∆ r (cid:48) · s ( r (cid:48) ) exp( − ik ˆ r (cid:48) · r ) and equal to atransverse plane wave with wave vector k = − k ˆ r (cid:48) . Itfollows that E ( r ) = E ( k , r ) + (cid:88) mm (cid:48) G ( r − r m ) · T mm (cid:48) · E ( k , r m (cid:48) )(A4)The fields { E ( r n ) } vanish because the T -matrix has ear-lier been designed to do so. For r = r n we can extract thesingular term G ( r n , r n ) = − /t and define the “macro-scopic field” ˜E ( r n ) in the vicinity of the dipole n as theone scattered from all others, ˜E ( r n ) = E ( k , r n )+ (cid:88) m (cid:54) = n ; m (cid:48) G ( r n , r m ) · T mm (cid:48) · E ( k , r m (cid:48) )= E ( k , r n ) + (cid:88) m (cid:54) = n t G ( r n , r m ) · ˜E ( r m )= 1 t (cid:88) m T nm · E ( r m ) (A5)The omission of the diagonal term m = n gives a goodimpression of the electric field inside the medium andthe solution of Eq. (A5) is equivalent to the calculationof the matrix T mm (cid:48) as is apparent from the last identity.However, it misses completely the singular field scatteredby the dipole at r n . At a small distance x from dipole n the relation between the fields E ( r n ) and ˜E ( r n ) is E ( r n + x ) = ˜E ( r n + x ) + t G ( x ) · ˜E ( r n ) (A6)Especially the longitudinal part is strongly singular as x → /x but carries a finite energy flux andposes less a problem.To illustrate this consider first one electric dipole lo-cated at r = 0, and for which only one diagonal termexists. Equation (A1) reduces to G ( k, r , r ) = G ( k, 0) + t G ( k, r ) . According to the analysis that has led toEq. (18), the local density of states (per unit volume,here per interval dk = dω/c ) at position r is given by N ( k, r ) = − ( k/π )Im Tr G ( k, r , r ). Formally, N ( k, 0) =0, but for r (cid:54) = 0 we can identify longitudinal states closeto the dipole, and transverse states far away. The totalextra number of states due to the presence of the dipoleis d Π( k ) = dk (cid:90) d r [ N ( k, r ) − N ( k, r )]= − kdkπ Im t ( k ) Tr (cid:88) p G ( k, p )= − kdkπ Im t ( k ) (cid:88) p (cid:20) k + 2[( k + i + ) − p ] (cid:21) This clearly separates into a strongly diverging longitu-dinal and a regular transverse component. Regularizingthe first to Q k / π , with Q = k c /γ the quality fac-tor, as proposed in section II B, we find that for k withina line width of k d Π( k ) = dk π (cid:20) − Q Im t ( k ) + 12 Re t ( k ) (cid:21) (A7)The first term is missed by ignoring divergencies, andthus difficult to capture by a numerical simulation. Itlargely dominates near resonance and is Lorentzian as isthe cross-section. The second term describes the mod-ification of transverse energy density and can be inter-preted as the change in local refractive index due to thepresence of the dipole. Upon integrating over the entireresonance, using that (cid:90) ∞−∞ dk t ( k ) = − π iQ we find that only the first term survives and giving a totalnumber of extra states per dipole is (cid:82) d Π( k ) = 3, equalto the number of degrees of freedom associated with theoptical polarization.The analysis above can be straightforwardly gener-alized to M dipoles. This yields expressions such asEq. (23) for the longitudinal complex wave number K L ( ∞ ) and associated with the DOS of longitudinalstates, and a similar one for transverse waves. From theensemble-averaged Dyson Green’s function G ( r , r ) in theunbounded medium, given in Eq. (3), upon splitting offterms, singular as x → 0, we obtainTr (cid:104) G (0) (cid:105) = δ (0) K L ( ∞ ) + Tr D (0)+ (cid:88) p p K T ( p ) K T ( p ) − p − (cid:88) p p (A8)4The Lorentz contact term of the effective mediumemerges as Tr G L (0) = δ (0) /K L ( ∞ ) that can be reg-ularized as before. Since D (0) is a finite longitudinalcontribution, it will be neglected. The transverse diver-gence described by the last term is real-valued and playsno role for DOS and is also independent of dipole den-sity. The integral over transverse wave numbers can bedefined as − iK T / π , which would be the value if K T ( p )were independent on p . Before ensemble averaging, theGreen’s function satisfies Eq. (A1). Let us first focus onthe trace of the diagonal terms that average to (cid:42) Tr M (cid:88) m =1 G ( r − r m ) · T mm · G ( r m − r ) (cid:43) =Tr M (cid:88) m =1 V (cid:90) d rG ( r − r m ) · (cid:104) T mm (cid:105) · G ( r m − r )In the thermodynamic limit, the average (cid:104) T mm (cid:105) shouldbe independent of the dipole m and its position r m .Hence the integral over r can be converted to Fourierspace to become n Tr (cid:88) p G ( p ) · (cid:104) T mm (cid:105) = n (cid:20) δ (0)3 k + i πk (cid:21) Tr (cid:104) T mm (cid:105) The contact term that occurs in this expression must beidentified with the one in Eq. (A8), so that1 K L ( ∞ ) = 1 k + n k Tr (cid:104) T mm ( k ) (cid:105) (A9)valid as M, V → ∞ , at constant density n = M/V . Sincewe expect (cid:104) T mm ( k ) (cid:105) ∝ the trace compensates the fac-tor 3 in the denominator.The off-diagonal elements m (cid:54) = m (cid:48) in Eq. (A1) arenegligible for the longitudinal states [and identified as D (0)], but not for the transverse waves, (cid:90) d r G ,T ( r m (cid:48) − r ) · G ,T ( r − r m )= (cid:88) p G ,T ( p ) exp[ i p · ( r m − r m (cid:48) )] ≡ πk H T ( k r )with H T ( y ) = 12 (cid:26) ie iy ∆ y − ddy (cid:18) e iy iy + e iy − y (cid:19) (1 − ˆyˆy ) (cid:27) which is regular ( H T (0) = i/ n πk Tr (cid:88) m (cid:48) (cid:54) = m (cid:104) T mm (cid:48) · H T ( k r mm (cid:48) ) (cid:105) Again, we suppose that the sum over m (cid:48) not to depend on m in the thermodynamic limit. Comparing to Eq. (A8) this gives the following expression for the complex trans-verse wave number K T = k − n k Tr (cid:104) T mm (cid:105) + i n k Tr (cid:88) m (cid:48) (cid:54) = m (cid:104) T mm (cid:48) · H T ( k r mm (cid:48) ) (cid:105) (A10)If we neglect any recurrent scattering (r.s.) from twoor more dipoles, we have T mm = t and T mm (cid:48) = t G ( r mm (cid:48) ). Converting the sum over m (cid:48) into the in-tegral n (cid:82) d r Tr G ( r ) · H ( r ) = in/ k gives1 K L ( ∞ ) = 1 k + ntk + r . s .K T = (cid:18) k − nt k − n t k + r . s . (cid:19) = k − nt + r . s . and we recover the ISA approximation. In particular, theoff-diagonal elements in Eq. (A10) are not negligible.We emphasize that the complex wave numbers K T and K L relate to transverse and longitudinal DOS and shouldnot be interpreted as effective medium parameters of elec-tromagnetic excitations. Appendix B: Numerical simulation of scattering andtransport mean free paths We consider a sample having the shape of a cylinderof radius R and thickness L parallel to the z axis of thereference frame and confined between the planes z = 0and z = L (see the inset of Fig. 10). The sample is madeof M point-like resonant scatterers described by Eq. (10)with a resonant frequency ω = k c and a decay rate γ (cid:28) ω of the excited state. The scatterers are locatedat random positions r j , j = 1 , . . . , M , inside the sample.The scatterer number density is n = M/V with V = πR L being the volume of the sample. In the following,we set k L = 10 and k R = 30, which implies M =2827–14137 for n/k = 0 . n/k = 0 . 02 atwhich we set k L = 30, k R = 60 and M = 6786.It is convenient to introduce dimensionless quantitiesand to neglect the frequency dependence of the Green’stensor over the bandwidth of interest that is assumed tobe much less than ω though can exceed γ considerably.Thus we put G ( k ) = G ( k ). Given an incident wave E ( k, r ) and using Eqs. (A4) and (A5), we obtain equa-tions for the electric field E ( k, r ) at any point in space: E ( k, r ) = E ( k, r ) + t ( k ) M (cid:88) j =1 G ( k , r − r j ) · ˜E ( k, r j )(B1)with G ( k , r ) = − e ik r πr [ P ( k r ) ∆ r + Q ( k r ) ˆrˆr ] (B2)5 (cid:19) (cid:21) (cid:23) (cid:25) (cid:27) (cid:20)(cid:19) (cid:16) (cid:19)(cid:17)(cid:19)(cid:24)(cid:19)(cid:17)(cid:19)(cid:19)(cid:19)(cid:17)(cid:19)(cid:24)(cid:19)(cid:17)(cid:20)(cid:19)(cid:19)(cid:17)(cid:20)(cid:24)(cid:19)(cid:17)(cid:21)(cid:19)(cid:19)(cid:17)(cid:21)(cid:24) Figure 10. Average, on-resonance, electric field in the mostdense sample. Red and blue lines show the real and imaginaryparts of the field, respectively. The black solid line is the fit ofEq. (B8) to the numerical data; the dashed line is the resultfor the imaginary part obtained from the fit to the real part. the free-space Green’s tensor for r (cid:54) = 0. Here P ( x ) =1 − /ix − /x and Q ( x ) = 2 /ix + 2 /x .The magnetic field B ( k, r ) can be found by applying B ( k, r ) = ∇ × E ( k, r ) / ( ik ) to Eq. (B1): B ( k, r ) = B ( k, r ) − it ( k ) M (cid:88) j =1 G ( B )0 ( k , r − r j ) · ˜E ( k, r j )(B3)where G ( B )0 ( k , r ) = − k π ( (cid:15) · ˆr ) e ik r [1 − P ( k r )] (B4)with the transverse, antisymmetric matrix( (cid:15) · ˆr ) = z/r − y/r − z/r x/ry/r − x/r (B5)and r = x ˆx + y ˆ y + z ˆ z . Note that G ( B )0 ( k , r ) only di-verges as 1 /r for small r , whereas G ( k , r ) divergesas δ ( r ) + 1 /r . In addition, the angular integral of G ( B )0 ( k , r ) vanishes. For these reasons, B ( k, r ) suffersfrom less fluctuations and is better suitable for numeri-cal studies.Equations (B1) and (B3) can be written for the fields ˜E ( k, r m ) and ˜B ( k, r m ) at each scatterer, excluding thefields scattered by themselves: ˜E ( k, r m ) = E ( k, r m )+ t ( k ) M (cid:88) j (cid:54) = m G ( k , r m − r j ) · ˜E ( k, r j ) (B6) ˜B ( k, r m ) = B ( k, r m ) − it ( k ) M (cid:88) j (cid:54) = m G ( B )0 ( k , r m − r j ) · ˜E ( k, r j ) (B7) (cid:212)(cid:212)(cid:212)(cid:212)(cid:212)(cid:212)(cid:212)(cid:212)(cid:212)(cid:212)(cid:212)(cid:212)(cid:212)(cid:212)(cid:212)(cid:212)(cid:212)(cid:212)(cid:212)(cid:212)(cid:212)(cid:212)(cid:212)(cid:212)(cid:212)(cid:212)(cid:212)(cid:212)(cid:212)(cid:212)(cid:212)(cid:212)(cid:212)(cid:212)(cid:212)(cid:212)(cid:212)(cid:212)(cid:212)(cid:212)(cid:212)(cid:212)(cid:212)(cid:212)(cid:212)(cid:212)(cid:212)(cid:212)(cid:212) (cid:16) (cid:22) (cid:16) (cid:21) (cid:16) (cid:20) (cid:19) (cid:20) (cid:21) (cid:22) (cid:16) (cid:20)(cid:19)(cid:20)(cid:21)(cid:22) Figure 11. The effective wave number obtained from the fitto the coherent field (see Fig. 10), the real part of K L ( ∞ )calculated using Eq. (A9), and the ISA result. We solve the system of equations (B6) for ˜E ( k, r m )( m = 1 , . . . M ) assuming that the sample is illuminatedby an incident linearly polarized plane wave: E ( k, r ) =ˆ x exp( ik z ). Magnetic fields on the scatterers and elec-tric and magnetic fields everywhere in space can be thenfound from Eq. (B7) and Eqs. (B1), (B3), respectively. 1. Scattering mean free path and effective wavenumber of coherent wave To determine the scattering mean free path (cid:96) , the so-lution of Eqs. (B6) is averaged over many (up to 10 ) in-dependent configurations of scatterers inside the sample,over slices of width ∆ z = L/ 100 along the z axis, and overthe central part of the cylinder with radius R = R − L (see the inset of Fig. 10). The average field (cid:104) E ( k, z ) (cid:105) ob-tained in this way should mimic the average field in aslab of infinite lateral extent ( R → ∞ ). A typical resultobtained from these calculations is illustrated in Fig. 10.To determine the scattering mean free path (cid:96) and theeffective wave number k e of the transverse waves, we fitthe results for the real part of (cid:104) E x ( k, z ) (cid:105) to the expressionRe (cid:104) E x ( k, z ) (cid:105) = A cos( k e z + φ ) exp (cid:16) − z (cid:96) (cid:17) (B8)where (cid:96) , k e , φ and A are free fit parameters (see Fig.10). In order to reduce the influence of boundary effects,we ignore the data corresponding to z < L/ 10 and z > (9 / L in the fits. The resulting effective wave numberand scattering mean free path are shown in Figs. 11 and12 (red line) as functions of detuning δ = ( ω − ω ) /γ for n/k = 0 . 2. Diffuse field We compute the average energy density of light insidethe sample ˜ ρ ( k, z ) by averaging the square of the macro-6 (cid:194) (cid:194) (cid:194) (cid:194) (cid:194) (cid:194) (cid:194) (cid:194) (cid:194) (cid:194) (cid:194) (cid:194) (cid:194)(cid:212)(cid:212)(cid:212)(cid:212)(cid:212)(cid:212)(cid:212)(cid:212)(cid:212)(cid:212)(cid:212)(cid:212)(cid:212)(cid:212)(cid:212)(cid:212)(cid:212)(cid:212)(cid:212)(cid:212)(cid:212)(cid:212)(cid:212)(cid:212)(cid:212)(cid:212)(cid:212)(cid:212)(cid:212)(cid:212)(cid:212)(cid:212)(cid:212)(cid:212)(cid:212)(cid:212)(cid:212)(cid:212)(cid:212)(cid:212)(cid:212)(cid:212)(cid:212)(cid:212)(cid:212)(cid:212)(cid:212)(cid:212)(cid:212) (cid:16) (cid:22) (cid:16) (cid:21) (cid:16) (cid:20) (cid:19) (cid:20) (cid:21) (cid:22)(cid:19)(cid:17)(cid:19)(cid:19)(cid:17)(cid:24)(cid:20)(cid:17)(cid:19)(cid:20)(cid:17)(cid:24)(cid:21)(cid:17)(cid:19)(cid:21)(cid:17)(cid:24)(cid:22)(cid:17)(cid:19) Figure 12. Scattering and transport mean free paths forthe most dense sample. We also show the scattering length1 / K L ( ∞ ) associated with longitudinal waves. (cid:19) (cid:21) (cid:23) (cid:25) (cid:27) (cid:20)(cid:19)(cid:19)(cid:20)(cid:21)(cid:22)(cid:23)(cid:24)(cid:25) Figure 13. The average energy density, consisting of a co-herent and a diffuse part. The solid straight line is a linearfit to the diffuse energy density for k z = 4–9. scopic magnetic field ˜B ( k, r m ) on the scatterers. Thediffuse energy density is obtained by subtracting the co-herent intensity:˜ ρ ( k, z ) = c π (cid:104)| ˜B ( k, r m ) | (cid:105) (B9)˜ ρ dif ( k, z ) = ˜ ρ ( k, z ) − c π |(cid:104) ˜B ( k, r m ) (cid:105)| (B10)where, as previously, the averaging (cid:104) . . . (cid:105) is done over scat-terer configurations as well as over the central part of thecylindrical sample. Because of equipartition, electric andmagnetic energies should be equal on average. Neverthe-less, because we calculate ˜ B , and not B , the magneticenergy density still misses the singular stored energy in-side the dipole.A typical profile of energy density inside the sampleis shown in Fig. 13. We fit ˜ ρ dif ( k, z ) by a linear func-tion for k z between 4 and 9 to determine its gradient d ˜ ρ dif ( k, z ) /dz .The average energy flux is given by the average Poynt- (cid:16) (cid:21) (cid:19) (cid:21) (cid:23) (cid:25)(cid:19)(cid:20)(cid:21)(cid:22)(cid:23)(cid:24) Figure 14. Transport mean free paths for different densi-ties n . The dashed line shows the ISA results for the lowestdensity 4 πn/k = 0 . ing vector (cid:104) K ( k, z ) (cid:105) = c π Re (cid:104) E ( k, r ) × ¯B ( k, r ) (cid:105) (B11)We compute the z component of (cid:104) K (cid:105) outside scatterersusing Eqs. (B1)–(B4) and extend the calculation to thespace in front and behind the sample. Inside the sample,the calculation fails to average because of the large fluctu-ations stemming from the near fields of the M scatterers.However, the calculation converges very well when aver-aging the Poynting vector calculated outside the sample.However, only in a slab of infinite lateral extent ( R → ∞ )without lateral leakage, the energy flux (cid:104) K z (cid:105) would beindependent of z and equal in- and outside the sample.For finite R , this equality is only valid approximately.To correct for this, we perform a linear fit of (cid:104) K z ( k, z ) (cid:105) calculated at z/L ∈ [ − . , − . 1] in front of the sample,and z/L ∈ [1 . , . 3] just behind the sample, and use thefit to find the value of (cid:104) K z ( k, z ) (cid:105) at z/L = 0 . 65 as thebest estimate of (cid:104) K z (cid:105) for the slab of infinite lateral ex-tent. The point z/L = 0 . 65 is chosen in the middle ofthe depth range where ˜ ρ dif ( k, z ) is seen to exhibit a clearlinear decay (see Fig. 13).The transport mean free path (cid:96) ∗ is obtained by usingthe Fick’s law (cid:104) K z (cid:105) = − D ddz ˜ ρ dif ( k, z ) (B12)where D = ( c /v p ) (cid:96) ∗ / v p = c k /k e is the phase velocity. No energy velocityappears here since ˜ ρ dif ( k, z ) does not count the storedenergy. Expressing (cid:96) ∗ from this equation yields (cid:96) ∗ = − k k e × (cid:104) K z (cid:105) d ˜ ρ dif ( k, z ) /dz (B13)The results following from this equation are shown inFigs. 12 (blue line) and 14. The comparison of theseresults with the analytic theory is presented in Fig. 9and is discussed in the main text.7Figure 12 shows that the transport mean free path dif-fers significantly from any of the scattering mean freepaths, including the scattering length 1 / K L ( ∞ ) as-sociated with longitudinal waves. The transport mean free path is an asymmetric function of the detuning fromthe resonance, and it is larger for positive than for nega-tive detunings. [1] For a review see: G.D. Scholes, Annu. Rev. Phys. Chem. , 57 (2002).[2] L.S. Levitov, Phys. Rev. Lett. , 547 (1990); EPL ,83 (2007).[3] C. Menottia, M. Lewenstein T. Lahaye and T. Pfau, Dipolar interaction in ultra-cold atomic gases , in: AIPConference Proceedings , 332 (2008).[4] Th. M. Nieuwenhuizen, A. Burin, Yu. Kagan and G.V.Shlyapnikov, Phys. Lett. A , 360 (1994).[5] S.E. Skipetrov, I.M. Sokolov, Phys. Rev. Lett. ,023905 (2014).[6] A. Lagendijk and B.A. van Tiggelen, Phys. Rep. ,143–215 (1996).[7] P. Sheng, Introduction to Wave Scattering, Localization,and Mesoscopic Phenomena (Academic, 1995).[8] Th.M. Nieuwenhuizen, A. Lagendijk, B.A. van Tiggelen,and A. Tip Phys. Lett. A , 191–194 (1992).[9] L. Allen and J.H. Eberly, Optical Resonance and Two-Level Atoms (Dover, 1987).[10] M. Rusek, A. Orlowski, and J. Mostowski, Phys. Rev. E , 4892 (1997).[11] F.A. Pinheiro, M. Rusek, A. Orlowski, and B. A. vanTiggelen Phys. Rev. E , 026605 (2004).[12] O. Leseur, R. Pierrat, J.J. Saenz, and R. Carminati,Phys. Rev. A , 053827 (2014).[13] B.A. van Tiggelen, A. Lagendijk and A. Tip, J.Phys-Condensed Matter , 7653–7677 (1990).[14] O. Morice, Y. Castin, and J. Dalibard Phys. Rev. A ,3896 (1995).[15] B.A. van Tiggelen and A. Lagendijk, Phys. Rev. B ,16729 (1994).[16] F. Eckert, A. Buchleitner, and T. Wellens, J. Phys. A:Math. Theor. (cid:80) p /p discussed in our Eq. (9) shows up in Eq. (13) of thisreference.[17] N. Cherroret, D. Delande, B. A. van Tiggelen Phys. Rev.A , 012702 (2016).[18] A. Lagendijk, B. Nienhuis, B.A. van Tiggelen, P. de Vries,Phys. Rev. Lett. , 533 (1997).[19] C.C. Kwong, D. Wilkowski, D. Delande, and R. Pierrat, Phys. Rev. A , 043806 (2019).[20] G.D. Mahan, Many-Particle Physics , third edition(Kluwer/Plenum, 2000).[21] P.W. Miloni, The Quantum Vacuum (Academic, 1994),Chapter 7.[22] G. Labeyrie, R. Kaiser, and D. Delande, Applied PhysicsB - Laser and Optics , 1001-1008 (2005).[23] Yu. N. Barabanenkov and V.D. Ozrin, Phys. Lett. A ,38 (1991).[24] Yu. N. Barabanenkov and V.D. Ozrin, Phys. Lett. A ,116 (1995).[25] B.A. van Tiggelen and R. Maynard in: Wave Propagationin Complex Media , edited by G. Papanicolaou (Springer-Verlag, 1998), p. 247-271.[26] M. Piraud, L. Sanchez-Palencia, and B.A van TiggelenPhys. Rev. A , 063639 (2014).[27] The electromagnetic conductivity is defined by the lefthand side of Eq. (44), in analogy with electric conductiv-ity. For a slab with depth L the total average transmissionof one channel a is (cid:104) (cid:80) b T ab (cid:105) = 4 (cid:96) ∗ / L . The number ofchannels is equal to N = 2 × k e A/ π , with A the lat-eral surface, k e the effective medium wave number, andthe factor two stemming from the two polarizations. Thisgives (cid:104) (cid:80) b T ab (cid:105) = 4 πN ( k ) D ( k ) A/L = 4 σA/L , with den-sity of states N ( k ) = k e / π v E and diffusion constant D ( k ) = v E (cid:96) ∗ as discussed in the text.[28] J.D. Jackson, Classical Electrodynamics (Wiley, 1999).[29] R. Loudon, L. Allen, and D. F. Nelson Phys. Rev. E ,1071 (1997).[30] B.A. van Tiggelen, Phys. Rev. A , 053826 (2019).[31] D. Vollhardt and P. W¨olfe, in : Electronic Phase Transi-tions , edited by W. Hanke and Yu. V. Kopaev (Elsevier,1992).[32] E. Akkermans and G. Montambaux, Mesoscopic Physicsof Electrons and Photons (Cambridge, 2011).[33] Z.Q. Zhang and P. Sheng, in : Scattering and Localizationof Classical waves in Random Media , edited by P. Sheng(World Scientific, 1990).[34] S.E. Skipetrov, I.M. Sokolov, Phys. Rev. B , 064207(2018); S.E. Skipetrov Phys. Rev. B94