Diffusion of monochromatic classical waves
Sijmen Gerritsen and Gerrit E. W. BauerKavli Institute of Nanoscience, Delft University of Technology, Lorentzweg 1, 2628CJ Delft, The Netherlands
共Received 27 July 2005; published 26 January 2006兲
We study the diffusion of monochromatic classical waves in a disordered acoustic medium by scattering theory. In order to avoid artifacts associated with mathematical point scatterers, we model the randomness by small but finite insertions. We derive expressions for the configuration-averaged energy flux, energy density, and intensity for one-, two-, and three-dimensional 共3D兲 systems with an embedded monochromatic source using the ladder approximation to the Bethe-Salpeter equation. We study the transition from ballistic to diffusive wave propagation and obtain results for the frequency dependence of the medium properties such as mean free path and diffusion coefficient as a function of the scattering parameters. We discover characteristic differences of the diffusion in 2D as compared to the conventional 3D case, such as an explicit dependence of the energy flux on the mean free path and quite different expressions for the effective transport velocity.
DOI:10.1103/PhysRevE.73.016618 PACS number共s兲: 42.25.Dd, 05.60.⫺k, 43.35.⫹d
I. INTRODUCTION
The ongoing interest in the field of classical waves in complex media is caused by the importance of detection and imaging techniques based on wave propagation and scatter-ing. This ranges from electromagnetic waves in optical and near infrared tomography 关1兴 and microwave radars 关2兴 to acoustic waves in ultrasonics关3兴 and geophysics 关4兴. Com-plexity is often associated with inhomogeneities that cause scattering which considerably complicates most imaging processes. However, when used cleverly, the scattered field can also be used to improve imaging 关5兴. Although length scales共with respect to the wavelength兲 and the degree of the disorder may vary considerably from field to field, methods and results have been shown to be interchangeable without much difficulty关6兴. Recent topics of interest include local-ization of classical waves关7,8兴, the transition from ballistic to diffusive wave propagation关9兴, acoustic time-reversal im-aging关10兴, etc. Direct simulation by the exact solution of a well-known Helmholtz wave equation for a given realization of the medium is often the method of choice for given appli-cations. The drawbacks of the brute force computational ap-proach are the limited system size and statistics that can be achieved with given computer resources as well as the diffi-culty to distill general principles out from the plethora of output data. The need for simple models with transparent results therefore remains.
An analytic theory of wave propagation in disordered me-dia necessarily relies on simple model scatterers, for which point scatterers, i.e.,共regularized兲 ␦ functions in real space, are often chosen 关11,12兴. Unfortunately, the scattering re-sponse of a single point scatterer can become noncausal, a pathological behavior that can not be solved by a simple momentum cutoff关13兴. Especially for the study of the fre-quency dependence over a wider range it is therefore neces-sary to use more realistic model scatterers.
In this paper we study a simple but not unrealistic experi-ment for the determination of the scattering properties of scalar waves in a disordered bulk material. A signal is emit-ted by a source and detecemit-ted by a receiver, both embedded in the medium at sufficiently large distances from the
bound-aries. Ultimately, we are interested in the detector signal caused by a pulsed共broadband兲 signal emitted by the source. After a first arrival we then expect the so-called coda that arrives at later times due to multiple scattering at the random scatterers关14兴. However, combining both the effects of mul-tiple scattering and the full frequency dependence of the scattering processes renders an analytical treatment difficult without additional approximations, such as a complete ne-glect of the frequency dependence of the scattering ampli-tudes when fitting the diffusive halo. In order to understand how to justify certain approximations and eventually find better ones, we have carried out a study of the frequency dependence of the scattering properties of random media. We concentrate on the steady state in the presence of strictly monochromatic sources, which distinguishes the present work from related studies of the propagation of narrow band pulses关15,16兴. As main results we obtain the frequency de-pendence of macroscopic effective medium properties such as the mean free path and the diffusion constant that depend on the microscopic parameters of the random scatterers.
When the ratio between source-receiver distance and mean free path is small, wave propagation is predominantly ballistic. When this ratio is large, energy and intensity propa-gation is governed by the diffusion equation 关12,15兴. Both these regimes are well understood. However, many imaging applications operate on length scales where the mean free path and the source-receiver distance are comparable. This is especially the case in geophysics where mean free paths range from a few hundred meters up to tens of kilometers 关17兴. The behavior at this crossover regime between ballistic and diffuse wave共intensity兲 propagation is of considerable interest关9兴 and also the subject of the present study.
source-receiver distance increases. We did not find many the-oretical studies of wave propagation in two-dimensional ran-dom media in the literature关6,19兴, although several experi-ments on quasi-2D systems have been carried out关15,20兴. Another possible test for our 2D theory is comparison with numerical studies, which for very large systems are much cheaper than in the 3D case.
The remainder of this paper is organized as follows. In Secs. II–IV we start by defining our model system and the basic equations, addressing the scattering matrices of single scatterers and discussing the average amplitude propagators in the frequency domain. The intensity, energy flux, and en-ergy density are discussed in Sec. V. Results on the fre-quency dependence of the diffusion constant and its depen-dence on the model parameters are discussed in Sec. VI. Generally, the results for 1D systems are easily obtained, whereas our results for 3D systems agree with findings pre-viously reported by others. The mathematics in the 2D case is not trivial, however, and the derivations are summarized in the Appendix. We end with the conclusions.
II. DEFINITIONS AND BASIC EQUATIONS A. Microscopic equations
We describe the propagation of共scalar兲 acoustic waves in a microscopic model system. Specifically, we consider a 1D, 2D, or 3D acoustic medium with wave velocity c0 and a mass density 0. The medium contains n randomly distrib-uted scatterers per unit length, area, or volume and we treat the dilute limit in which the average distance between scat-terers is much larger than their radius a. The internal wave velocity of a scatterer is cintand, for simplicity, the difference in mass density with the surrounding medium is disregarded. The waves are emitted by a monochromatic point source oscillating at frequencypositioned at the origin. The wave amplitude , related to the pressure by p=t and to
local particle velocity by v= −0−1, then obeys the wave equation
关2− c−2共r兲
t
2兴
共r;t兲 = − Q0␦共d兲共r兲cos共t兲. 共1兲 The source term chosen here corresponds to a volume injec-tion term, with␦共d兲the Dirac delta function and d the dimen-sion. The source emits plane waves for 1D, cylindrical waves for 2D and spherical waves for 3D media. In all cases Q is in units of length per unit time. The wave velocity profile of the entire medium c共r兲 contains the information of the positions of the scatterers共in 1D r=x兲.
The Green function of the Helmholtz equation共1兲 in the real space and frequency domain reads
关2+ 0
2− V共r;兲兴G共r,r
⬘
;兲 =␦共d兲共r − r
⬘
兲, 共2兲 where0=/ c0, the length of the wave vector in the homo-geneous medium. V共r;兲 is the scattering or impurity poten-tial, a sum over all individual scattering potentialsV共r;兲 =02共1 −␥−2兲
兺
i=1 N
⌰共a − 兩r − ri兩兲. 共3兲
⌰ is the Heaviside step function, with ⌰共x兲=0 when x⬍0 and 1 otherwise. The velocity contrast is defined as ␥= cint/ c0so that the single scatterer potential is “attractive” when ␥⬍1 and “repulsive” when ␥⬎1. Equation 共3兲 de-scribes a spherical potential, however, the precise shape is not relevant when the scatterers are sufficiently small com-pared to the wavelength.
The amplitude of the wave field is related to the Green function
共r;t兲 = − Q0Re兵e−itG共r,r
⬘
= 0;兲其. 共4兲 The intensity I共r;t兲 is the square of this expression. Related physical properties are the energy fluxF共r;t兲 = − 1 0
t共r;t兲共r;t兲, 共5兲
and the energy density
W共r;t兲 = 1
20
关共r;t兲兴2+ c−2共r兲关t共r;t兲兴2, 共6兲
recognized as the sum of the potential and kinetic energy contributions, respectively. For a monochromatic source with frequency these observables contain a time independent contribution and a second term oscillating with frequency 2. We concentrate on the constant part by time averaging over one period. Expressed in terms of the Green function this yields I共r兲 =Q 2 0 2 2 兩G共r,r
⬘
= 0;兲兩 2, 共7兲 F共r兲 = −Q 2 0 2 Im兵G共r,r⬘
= 0;兲G *共r,r⬘
= 0;兲其, 共8兲 W共r兲 =Q 2 0 4冉
兩G共r,r⬘
= 0;兲兩 2+ 2 c2共r兲兩G共r,r⬘
= 0;兲兩 2冊
. 共9兲 B. Macroscopic equationsThe properties of the wave field depend, via the Green function, on the exact configuration of scatterers. However, in large systems, different realizations of the ensemble give similar responses 共ergodicity兲. The similarities in the re-sponse can be studied by calculating the configurational av-erage. This average is the connection between the micro-scopic description and the macromicro-scopic 共effective兲 medium properties.
t具I共r;t兲典 = D2具I共r;t兲典, 共10兲
where the brackets denote the configuration average and D is the diffusion constant. In spite of neglecting the frequency dependence of the diffusion constant in this case, this ap-proximation is known to work well in cases where the source-receiver distance is much larger than the mean free path and the incoming pulse is a narrowband signal关12,15兴. In the case of a narrowband pulse, one can also write a transport equation for a wave packet with some inner and outer frequencies. In this way the frequency dependence of
D can be derived关16兴.
In order to obtain the steady-state diffuse intensity of a monochromatic wave field, the diffusion equation共10兲 is in-sufficient. The energy density共and not the intensity兲 of the wave field is the conserved property. Equation共10兲 is there-fore only valid if the intensity is strictly proportional to the energy density. In general, the averaged energy transport is governed by Fick’s first law
具F共r兲典 = − D共兲具W共r兲典, 共11兲 accounting for the frequency dependence of the diffusion constant. In the steady-state problem and outside the mono-chromatic source the proper Laplace equation is
2具W
共r兲典 = 0. 共12兲
III. SCATTERING MATRICES
Here we discuss the properties of a single model scatterer in the system 关N=1 in Eq. 共3兲兴. The response of a system containing a monochromatic source共in the origin兲, a receiver 共at r兲 and a single “s-wave” scatterer 共at ri兲 can be expressed
in terms of Green functions of the homogeneous system 共V=0兲 关21兴:
G共r,r
⬘
= 0;兲 = G0共r;兲 + G0共兩r − ri兩;兲t0共兲G0共r;兲. 共13兲 This expression is valid in the far field limit 共r,riⰇ兲 andwhen scattering is isotropic 共Ⰷa兲, where is the wave-length.
The transition 共t-兲 matrix elements for s-wave scattering are related to the scattering matrix elements by
t0共兲 =
冦
2i0R共兲 共1D兲, 2i关S0共兲 − 1兴 共2D兲, 2i0−1关S0共兲 − 1兴 共3D兲.
共14兲 In 1D, the s-wave scattering condition corresponds to equivalence of t0for either reflection or transmission. R共兲 is the reflection coefficient at a step discontinuity, and can be obtained by imposing flux conservation across the scatterer boundary. This gives关22兴
R共兲 = e−i02aR0共1 − e
i04a/␥兲
1 − R02ei04a/␥ , 共15兲
where R0=共␥− 1兲/共␥+ 1兲. By imposing the same condition we can derive an expression for the scattering matrix element
of the s-wave channel S0关related to the scattering phase shift ␦0by S0= exp共i2␦0兲兴. In 2D 关23兴 S0共兲=−␥J0共0a/␥兲H1 共2兲共 0a兲 − J1共0a/␥兲H0共2兲共0a兲 ␥J0共0a/␥兲H1共1兲共0a兲 − J1共0a/␥兲H0共1兲共0a兲 . 共16兲 In 3D the Bessel 共Ji兲 and Hankel 共Hi
共j兲兲 functions are re-placed by the spherical Bessel 共ji兲 and Hankel 共hi
共j兲兲 func-tions. The scattering matrix element then simplifies to关22兴
S0共兲 = e−i20acot共0a/␥兲 + i␥
cot共0a/␥兲 − i␥. 共17兲
In this calculation of the scattering matrices, the differ-ence in mass density is disregarded. However, including this does not fundamentally alter the calculation 共the scattering matrix is still calculated from flux conservation兲. Further-more, the scattering matrices calculated here describe acous-tic wave scattering where only acousacous-tic modes are allowed inside the scatterers. When solid scatterers are considered extra mode conversions from acoustic waves to shear waves back to acoustic waves occur which considerably compli-cates the calculation. In principle, this calculation can be done 关24兴 and it is known that the extra mode conversions cause extra resonances in scattering properties of the scatter-ing object关25兴.
IV. THE CONFIGURATION-AVERAGED PROPAGATOR Now we switch to the case of multiple scattering at the proposed model scatterers. The wave propagator in a disor-dered medium after configuration averaging is dressed with a self-energy⌺. In reciprocal space it reads 关12兴
具G共k,k
⬘
;兲典 = 1 0 2 − k2−⌺共k;兲共2兲 d␦共d兲共k − k⬘
兲. 共18兲 When n, the density of scatterers, is low, interference be-tween multiply scattered waves by different sites may be disregarded. In this “single site approximation” the self-energy does not depend on k and it is simply given by关12兴⌺共兲 = nt0共兲. 共19兲
关G共r;兲=具G共r,r
⬘
= 0 ;兲典兴. In 1D the amplitude propagators are exponentially damped plane wavesG共兩x兩;兲 = 1
2ie共兲
eie共兲兩x兩, 共20兲
in 2D they are cylindrical:
G共r;兲 =
冦
− i 4H0 共1兲„ e共兲r… if⬎ 0, i 4H0 共2兲„− e共兲r… if⬍ 0, 共21兲and in the 3D spherical
G共r;兲 = − 1
4re
ie共兲r 共22兲
关6兴. In Eqs. 共20兲–共22兲 e is the “renormalized” effective
wave vector
e共兲 =
冑
02− nt0共兲 ⬅ sgn共兲r共兲 + i1 2lf共兲
. 共23兲 r共兲=兩Re兵e共兲其兩 and lf−1共兲=2兩Im兵e共兲其兩, the mean free
path. We retrieve the Green functions for the homogeneous systems共G0兲 by letting n or t0 go to zero. Properties of the averaged response to a pulsed signal can be studied by cal-culating the Fourier transform to the time domain, as was done in Refs.关13,26兴.
V. THE CONFIGURATION-AVERAGED INTENSITY END ENERGY
We derive here the configuration averaged intensity, en-ergy flux, and enen-ergy density in the frequency domain.
A. The Bethe-Salpeter equation Ensemble averaging the intensity共7兲 gives us
具I共r兲典 =Q20 2
2 ⌸共r;兲, 共24兲
where ⌸共r;兲=具兩G共r,r
⬘
= 0 ;兲兩2典 is the average of the squared Green function propagator. It is given by⌸共r;兲 = ⌸0共r;兲 +
冕
ddr1ddrddr3ddr4具G共r,r1;兲典 ⫻具G*共r,r2;兲典⌫共r1,r2,r3,r4;兲 ⫻ 具G共r3,r
⬘
= 0;兲G*共r4,r
⬘
= 0;兲典. 共25兲This is the Bethe-Salpeter equation in position space, where ⌸0 is the coherent intensity 共⌸0=兩具G典兩2兲 and ⌫ is the irre-ducible vertex function. The lowest order approximation that still accounts for multiple scattering is
⌫共r1,r2,r3,r4;兲
= n⌫共兲␦共d兲共r
1− r3兲␦共d兲共r1− r2兲␦共d兲共r3
− r4兲, 共26兲
reducing the Bethe-Salpeter equation to
⌸共r;兲 = ⌸0共r;兲 + n⌫共兲
冕
ddr1⌸0共兩r − r1兩;兲⌸共r1;兲. 共27兲 In reciprocal space this integral equation becomes a geomet-ric series that can be summed as⌸共k;兲 = ⌸0共k;兲
1 − n⌫共兲⌸0共k;兲. 共28兲 In order to be able to calculate the Fourier transform of ⌸共k;兲, an expression for ⌸0共k;兲 is needed. It is calculated as the Fourier transform of the coherent intensity and this results in 1D in ⌸0共k;兲 = 2lf 3 关共2rlf兲2+ 1兴关共klf兲2− 1兴 , 共29兲 in 2D in ⌸0共k;兲 = lf2 arcsin
冉
冑
共2rlf兲 2−共kl f兲2冑
1 +共2rlf兲2冊
冑
1 +共klf兲2冑
共2rlf兲2−共klf兲2 , 共30兲 and in 3D in关12兴 ⌸0共k;兲 = lf 4 arctan共klf兲 klf . 共31兲The calculation of the vertex function⌫ is discussed in the next subsection.
B. Energy conservation and the Ward identity It is well known that for a given approximation for the self-energy, the vertex correction cannot be freely chosen. Here we take advantage of the flux conservation constraint to obtain ⌫ without additional microscopic calculations. The energy flux from the monochromatic source 共on average兲 points outwards. In the steady-state case the following con-dition must hold for the averaged flux in direction n:
具n · F共r兲典 ⬀ 1
rd−1n · rˆ, 共32兲
where rˆ is the unit vector in the radial direction. In 1D this condition reads
具F共x兲典 ⬀ sgn共x兲. 共33兲
The microscopic expression for the average energy flux is 具n · F共r兲典 = −Q 2 0 2 Im兵具G共r,r
⬘
= 0;兲n · G *共r,r⬘
= 0;兲典其 = −Q 2 0 2 Im兵⌸´ n共r;兲其, 共34兲which defines the function ⌸´n. The vertex function is the same as for the intensity, so we can express⌸´nin reciprocal
⌸´n共k;兲 = ⌸´0 n 共k;兲 1 − n⌫共兲⌸0共k;兲 . 共35兲 ⌸´0
n共k;兲 is the coherent energy flux in direction n that is
given by the Fourier transform of ⌸´0
n共r;兲 = G共r;兲n · G*共r;兲. 共36兲 In 2D and 3D the averaged microscopic expression for the energy flux should match the macroscopic condition
Im兵⌸´n共r;兲其 = − C
rd−1n · rˆ, 共37兲
which in reciprocal space reads
Re兵⌸´n共k;兲其 = − 共n · kˆ兲2d−1C
k, 共38兲
where C is real and depends on frequency and the model parameters.⌸0共k;兲 is an even function of k. We know how ⌸´0
n共k;兲 depends on k, as the Fourier transform in 2D reads
⌸´0 n 共k;兲 = − 共n · kˆ兲2i
冕
0 ⬁ drJ1共kr兲rG共r;兲rG*共r;兲, 共39兲 and in 3D ⌸´0 n 共k;兲 = − 共n · kˆ兲4i冕
0 ⬁ drj1共kr兲r2G共r;兲rG*共r;兲. 共40兲 The Taylor series of⌸´0n共k;兲 around k=0 only contains odd terms. So, in the limit that k→0, condition 共38兲 can only be fulfilled by Eq.共35兲 whenn⌫共兲 = ⌸0−1共k = 0;兲. 共41兲
In 1D showing that condition共33兲 can only be fulfilled when Eq. 共41兲 is fulfilled as well is straightforward. The Ward identities are relations between self-energy and vertex cor-rections. We can identify Eq.共41兲 as the Ward identity for our problem. We now have all the ingredients to calculate the Fourier transform of Eqs.共28兲 and 共35兲 to calculate the av-eraged intensity and energy flux respectively.
C. Flux
Using the Taylor expansions in the limit k→0, we find an expression for C关from Eq. 共37兲兴 in 2D and 3D:
C =
冕
0 ⬁ dr Im兵G共r;兲rG*共r;兲其rdd−1 1 2⌸0 −1共k = 0;兲 k 2⌸ 0共k;兲兩k=0 . 共42兲The average flux in 1D is obtained by directly Fourier trans-forming Eq.共35兲: 具F共x兲典 = Q20兩兩 8 r r2+ 1/共2lf兲2 sgn共x兲. 共43兲 We show how to calculate C in 2D case in the Appendix. With the result, the projection of the average flux becomes
具n · F共r兲典 =
Q20兩兩 82
arctan共2rlf兲
r n · rˆ, 共44兲
while in 3D calculating C from Eq.共42兲 is straightforward and the projection of the average flux then reads
具n · F共r兲典 =
Q20兩兩 2
r
共4兲2r2n · rˆ. 共45兲 Letting lf→ ⬁共r→兩0兩兲 recovers the flux of a monochro-matic source in an unperturbed medium.
It is interesting to see that, in contrast to the 3D case, in 1D and 2D the average flux depends on both the mean free path and the real part of the effective wave vector. So the scattering mean free path limits the energy flux in 1D and 2D, but not in 3D. In a strongly scattering 2D medium, in which the wave energy is not共yet兲 localized, the dependence on the arctangent should be observable.
D. Intensity
The total average intensity is proportional to the propaga-tor⌸共r;兲, which can be obtained by calculating the Fourier transform ⌸共r;兲 =
冕
d d k 共2兲d eik·r ⌸0 −1共k;兲 − ⌸ 0 −1共k = 0;兲. 共46兲 ⌸0共k兲 is given by Eq. 共29兲 in 1D, Eq. 共30兲 in 2D, and Eq. 共31兲 in the 3D case. In 1D and 2D this integral diverges because in the steady state case with a monochromatic source, energy does not escape fast enough to infinity due to the scatterers. This is analogous to the fact that the Poisson equation共the diffusion equation in steady state with source term兲 for a line or planar source has no well-defined solution. The gradient of the intensity exists in all cases. In 1D it is constant and the derivative of⌸共x;兲 is given byx⌸共x;兲 = − sgn共x兲 1 4lf 1 r2+ 1/共2lf兲2 . 共47兲
The gradient of⌸ in 2D is expressed as an integral by ⌸共r;兲 = − rˆ
冕
0 ⬁ dk 2 k2J1共kr兲 ⌸0−1共k;兲 − ⌸0−1共k = 0;兲 = rˆf共r;兲, 共48兲 which defines a function f共r;兲, that represents the gradient in the rˆ direction. We split this up into a coherent共coh兲 and a “totally diffusive”共td兲 part and a crossover correction 共cr兲fcoh共r;兲 =r兩G共r;兲兩2= − 1 8Re兵共r+ i/2lf兲H1 共1兲„共 r + i/2lf兲r… ⫻ H0共2兲„共r− i/2lf兲r…其. 共50兲
In the appendix it is shown that
ftd共r;兲 = −arctan共2rlf兲 22 rlf 1 rg −1共2 rlf兲, 共51兲 with g共2rlf兲 = 1 − 1 共2rlf兲2 + 1 2rlfarctan共2rlf兲 . 共52兲 This part decays as 1 / r, much slower than the coherent and crossover contributions. It is the part that describes the inten-sity gradient when energy transport is completely governed by Fick’s first law, so we refer to this term as the “totally diffusive” part. When the total gradient is approximated by just the sum of the coherent and the totally diffusive contri-bution, the gradient first decays exponentially until the source-receiver distance is approximately two to three mean free paths and then the 1 / r decay is dominant. However, in this approximation it is neglected that close to the source the diffusive field is different compared to the field far away from the source. The third term of f, the crossover term, describes this difference. In the Appendix it is shown that
ftd共r;兲 + fcr共r;兲=−
冕
0 ⬁ dk 2J1共kr兲k 2⌸ sc共k;兲, 共53兲 where ⌸sc共k;兲 = ⌸0 −1共k = 0;兲⌸ 0共k;兲 ⌸0−1共k;兲 − ⌸0−1共k = 0;兲 . 共54兲We did not find an analytical expression for the integral共53兲 and thus we need to evaluate it numerically. The crossover term vanishes for r / lf→0 or r/lfⰇ1 and peaks at
r / lf⬇0.3. Only around this value of r/lf, is the gradient共in
absolute value兲 overestimated significantly 共up to 25%兲 when
we approximate it by just the sum of coherent and “totally diffusive” terms.
In 3D the Fourier transform共46兲 converges and the inten-sity is well defined. We rewrite
⌸共r;兲 = 1 162rlf
冉
lf re −r/lf+ 3 + e−r/lfh共r/l f兲冊
, 共55兲 where h共r/lf兲 =冕
0 ⬁ d冉
4共+ 1兲 2 关2共+ 1兲 − ln共1 + 2/兲兴2+2− 1冊
e−r/lf. 共56兲 We were also not able to solve Eq. 共56兲 analytically. In the 3D case, the intensity is a function of lf only 共it doesnot depend on r兲. Equation 共55兲 consists of three terms
共⌸=⌸coh+⌸td+⌸cr兲. The first term is proportional to the co-herent intensity, the second term is the algebraically decay-ing diffuse term共the only term that is not exponentially de-caying兲. We plot the sum of the first and second term multiplied by r关in units of 1/共162l
f兲兴 in Fig. 1 as a function
of r / lf 共on the right axis兲. The third term is again the
cross-over correction to the total intensity if we approximate⌸ by only the first two terms. A plot of the crossover correction divided by the sum of the first two terms is shown in Fig. 1 共left axis兲. The crossover term vanishes for r/lf→0 or r/lf
Ⰷ1 and peaks at r/lf⬇0.3. It can thus be concluded that the
intensity can very well be approximated by just the sum of coherently and totally diffusively propagated intensities, as the total intensity in 3D is never overestimated by more than 5% using this approximation.
To complete this discussion we show the final results for the gradient of the average intensity in 1D and 2D:
x具I共x兲典 = − Q202 2 sgn共x兲 1 4lf兩e兩2 , 共57兲
FIG. 1. r共⌸coh+⌸td兲 in units of 共162lf兲−1
共dotted line, right axis兲 and ⌸cr/共⌸coh+⌸td兲
共solid line, left axis兲 as a function of r/lf 共the
具I共r兲典 ⬇ rˆQ20 2
2 关fcoh共r;兲 + ftd共r;兲兴, 共58兲 where fcoh and ftdare given by Eqs. 共50兲 and 共51兲, respec-tively. The average intensity in 3D is approximated well by
具I共r兲典 ⬇Q 2 0 2 2 1 162rlf
冉
lf re −r/lf+ 3冊
. 共59兲We obtain these expressions from our first-principles calcu-lations that enable us to study not only the ballistic and dif-fusive limits, but also the crossover regime when r / lf⬇1.
From this we observe that we can approximate the average intensity well by only the coherent and diffusive contribu-tions. Furthermore, we saw that already at r / lf⬇0.3 the
dif-fusive intensity is higher than the coherent intensity. This does not mean that when a pulsed source is used we should see signs of the crossover to the diffusive regime at this point because the diffuse peak is much broader than the coherent peak so this crossover point is at larger values of r / lfas was
previously reported关9兴. Obviously, our present model system has been assumed to be boundless. In a finite slab geometry boundary scattering, which is beyond the scope of this study, would of course affect the results.
E. Energy density
To derive a first-principles expression for the diffusion constant from Fick’s law共11兲, we still have to calculate the average energy density given by
具W共r兲典 =Q 2 0 4 共具兩G共r,r
⬘
= 0;兲兩 2典 +具2c−2共r兲兩G共r,r⬘
= 0;兲兩2典兲. 共60兲 The first term is the average potential energy density and the second term corresponds to the kinetic energy.We start with the potential energy term in 2D and 3D. We define
具兩ⵜG共r,r
⬘
= 0;兲兩2典 = ⌸˝共r;兲. 共61兲 The Fourier transform of ⌸˝0共=兩具ⵜG典兩2兲 diverges, which means that we can not use the same procedure as we used for the intensity and the flux. According to the Bethe-Salpeter equation⌸˝共r;兲 = ⌸˝0共r;兲 + ⌸0
−1共k = 0;兲
⫻
冕
ddr1⌸˝0共r1;兲⌸共兩r − r1兩;兲. 共62兲 This integral diverges as well because of the strong singu-larities in⌸˝0共also when the gradient is calculated in the 2D case兲. When averaging, scatterers are effectively moved around the medium, and for every configuration, the contri-bution to the total average response is calculated. However, because of the stronger singularities in⌸˝0共as every scatterer becomes a new source of spherical waves兲 this is not pos-sible when the receiver position coincides with a scattererposition. The reason for this is the point receiver assumption and the far field scattering approximation. We can circum-vent this problem by omitting a small volume/area around r1 with radius of approximately one wavelength. This slightly modifies the probability distribution function form “com-pletely random” to “nonoverlapping”共with the receiver兲 in order to avoid the divergencies. We then find that⌸˝ is given by:
⌸˝共r;兲 = 兩e兩2⌸共r;兲. 共63兲
In principle, our original expression for ⌸ should now be multiplied by a factor exp共−r0/ lf兲, where r0 is the radius of omission so as long as the mean free path is longer than a few wavelengths omitting this small volume does not influ-ence the results. Furthermore, even if scattering is strong and the mean free path is of the order of the wavelength, this factor is not important.
The second term of Eq.共60兲, the kinetic energy, can be split:
具2c−2共r兲兩G共r,r
⬘
= 0;兲兩2典=02⌸共r;兲 − 具V共r;兲兩G共r,r
⬘
= 0;兲兩2典. 共64兲 Now the condition that the scatterer position cannot coincide with the receiver position ensures that the second term van-ishes, due to the step function in the potential 共3兲. We can thus just disregard this term.In 1D proving that
⌸˝共兩x兩;兲 = 兩e兩2⌸共兩x兩;兲, 共65兲
always holds is straightforward. We have to impose the con-dition the receiver cannot coincide with a scatterer to ensure that
具2c−2共x兲兩G共x,x
⬘
= 0;兲兩2典 = 02⌸共兩x兩;兲. 共66兲 Only under the restrictions mentioned here, can the averaged energy density in 1D, 2D, and 3D be expressed as being proportional to the intensity
具W共r兲典 = 1 20
共兩e兩2+0 2兲具I
共r兲典 共67兲
and this thus means that only the gradient of the energy density is well defined in the 1D and 2D cases.
VI. THE DIFFUSION CONSTANT
Using the Bethe-Salpeter equation with the Ward identity we find expressions for the average energy flux 共43兲–共45兲, the共gradient of兲 the average intensity 共57兲–共59兲. The average energy density is just proportional to the average intensity 共67兲. When r/lfⰇ1 we expect Eq. 共11兲 to hold and, as the
D共兲 =1
dceff共兲lf共兲, 共68兲
where in the 1D and 3D case
ceff共兲 = c0 2r兩0兩 r 2+ 1/共2l f兲2+0 2 共69兲
and in the 2D case
ceff共兲 = c0 2r兩0兩
r
2+ 1/共2l
f兲2+0
2g共2rlf兲, 共70兲
where g共2rlf兲 is given by Eq. 共52兲. The effective transport
velocity in 2D reduces to Eq. 共69兲 in the weak scattering limit.
We can now investigate the frequency dependence of the diffusion constant for a medium with monodisperse scatterers. We relate the scatterer density n to the average distance between scatterers 共具ds典兲 so that n=具ds典−1 in 1D,
n = 4−1具d
s典−2 in 2D and n = 3共4兲−1具ds典−3 in 3D. Let us
focus on the diffusion constant of the 2D medium. We write
ae共a0兲 =
冑
共a0兲2− 4 冉
a 具ds典冊
2 t0共a0兲, 共71兲so that the dimensionless property ae depends on the
di-mensionless frequency0a共=a / c0兲 and two dimensionless
model parameters, i.e., the velocity contrast␥共=cint/ c0兲 and the average distance between scatterers in number of scat-terer radii共具ds典/a兲. The real and imaginary parts of ae are
needed to obtain the diffusion constant
ar共a0兲 = 兩Re兵ae共a0兲其兩, 共72兲
lf共a0兲
a =
1 2兩Im兵ae共a0兲其兩
. 共73兲
The diffusion constant for a 2D medium is plotted in Fig. 2. The relevant frequency range is from0a共=a / c0兲=0 to
0a⬇/ 2, as for higher frequencies the isotropic scatterer assumption is no longer valid. For the plot, the density of scatterers was determined by setting 具ds典/a=10, increasing
this value shifts the curves up. The shape of the curves is predominantly determined by the mean free path. The effec-tive transport velocity ceffonly deviates considerably from c0 when the scatterer velocity and the frequency are small and the scatterer density high. For the diffusion constants shown in the plot, this is only the case when␥= 0.2. This is also the only case that shows resonances in the relevant frequency range. Further lowering the internal velocity of the scatterers would “pull in” more resonances in the relevant frequency range. These resonances appear because of resonances in the mean free path. When the scatterer-medium velocity ratio is increased, the mean free path 共and thus the diffusion con-stant兲 increases until the ratio is larger than unity and then it drops again. However, increasing ␥ above 10 does not change the diffusion constant much in the frequency range we discuss.
The diffusion constants in 1D and 3D media show the same behavior. Of course, the resonances at low velocity, are caused by the fact that all scatterers are assumed to be of equal size. When solid scatterers in a fluid are considered even more resonances are expected to show up关25兴. When scatterer sizes共or velocities兲 are allowed to vary, the reso-nances are averaged out.
When the scatterer velocity is zero we obtain an impen-etrable model scatterer. This is not a useful model scatterer, as in the low frequency range the mean free path共and thus the diffusion constant兲 differ considerably from the pen-etrable scatterer case. The reason for this is that the limits for
→0 and␥→0 do not commute, as
lim
→0␥→0limlf= const, 共74兲 while
FIG. 2. Diffusion constant of the 2D disor-dered medium in units of c0a, as a function of the
dimensionless frequency 0a for four different
lim
␥→0→0limlf= ⬁ . 共75兲
The effect of this is that at the longer wavelengths lffor the
impenetrable scatterer is several orders of magnitude smaller than lf for nonzero values of␥.
VII. CONCLUSIONS
We have calculated the transport of energy and intensity in disordered 1D, 2D, and 3D共infinite兲 media emitted by a monochromatic source. Using the ladder approximation to the Bethe-Salpeter equation we explicitly show that the total intensity is well approximated by the sum of the coherent and the fully developed diffuse wave field for all source-receiver distances. The results for 3D disordered systems agree with findings previously reported, except for the ex-pression for the intensity in the crossover regime, which has not been reported before. We have obtained more new results studying energy and intensity propagation in the 2D system in detail. When compared to the 3D case, the 2D disordered system shows interestingly different behavior: In 2D, the av-erage energy flux depends on the mean free path and the effective transport velocity depends differently in terms of the scattering parameters. The共gradient of the兲 intensity as a function of the source-receiver distance, on the other hand, behaves similarly in the 2D and the 3D case. The monochro-matic source enables us to investigate the frequency depen-dence of the macroscopic diffusion constant where we par-ticularly focused on the influence of the finite size of the scatterers. For a monodisperse distribution of scatterers shape resonances show up in the relevant frequency range for low internal scatterer velocities 共␥ small兲. In this fre-quency range共where scattering is expected to be isotropic兲 the dependence of the scattering properties on frequency can-not be neglected. This means that descriptions of broadband pulse propagation through these media should in principle incorporate both frequency dependent and multiple scatter-ing effects. The development of a workable Ward identity in this case remains a challenge, however. Finally, we want to point out that our model describes transport of scalar acous-tic waves but results can be extended and many conclusions should also apply to vector wave fields random media.
Note added. After submission of this paper, more
theoret-ical and experimental work on the crossover regime between ballistic and diffusive wave propagation in 2D systems was published 关Y. Lai, S. K. Cheung, and Z. Q. Zhang, Phys. Rev. E 72, 036606共2005兲兴. This paper discusses wave propa-gation from pulsed planar sources in 2D systems and the results on the diffusive-ballistic transition in this case corre-spond to our findings for a monochromatic point source.
ACKNOWLEDGMENTS
We thank Mauro Ferreira and Kees Wapenaar for discus-sions. This work is part of the research program of the “Stichting Technische Wetenschappen” 共STW兲 and the “Stichting Fundamenteel Onderzoek der Materie” 共FOM兲 which is financially supported by the “Nederlandse Organi-satie voor Wetenschappelijk Onderzoek”共NWO兲.
APPENDIX: ENERGY AND INTENSITY IN 2D In this appendix we derive the configuration-averaged in-tensity and energy flux in a disordered 2D medium. Starting point is the 2D Green function propagator
G共r;兲 =
冦
− i 4H0 共1兲关„ r+ i/共2lf兲…r兴 if⬎ 0, i 4H0 共2兲关„ r− i/共2lf兲…r兴 if⬍ 0. 共A1兲We use the properties
H0共2兲关„r− i/共2lf兲…r兴 = H0共1兲关„−r+ i/共2lf兲…r兴 共A2兲 and H0共1兲关„±r+ i/共2lf兲…r兴 = − i 2 K0关„⫿ir+ 1/共2lf兲…r兴, 共A3兲 to express the Hankel functions 共H0共j兲兲 in terms of modified Bessel function of the second kind共K0兲. The Fourier trans-form of the coherent intensity
⌸0共k;兲 = 2
冕
0⬁
dr兩G共r;兲兩2J0共kr兲 共A4兲 is then obtained from Ref. 关27兴 and using properties of the associated Legendre polynomials关28兴 as
⌸0共k;兲 = lf 2 1 1 +共2rlf兲2 P−1/2−1/2共u兲 P1/2−1/2共u兲, 共A5兲 with u =1 −共2rlf兲 2 1 +共2rlf兲2 + 2klf 1 +共2rlf兲2 . 共A6兲
This is can be rewritten as
⌸0共k;兲 = lf2 arcsin
冉
冑
共2rlf兲 2−共kl f兲2冑
1 +共2rlf兲2冊
冑
1 +共klf兲2冑
共2rlf兲2−共klf兲2 . 共A7兲⌸0共k;兲 is real, continuous, and differentiable for all 共real兲
k艌0.
The flux in the 2D system is given by
具n · F共r兲典 =Q 2 0 2 C rn · rˆ, 共A8兲
where C is the constant to be calculated:
C =
冕
0 ⬁ dr Im兵G共r;兲rG*共r;兲其r2 ⌸0−1共k = 0;兲k 2⌸ 0共k;兲兩k=0 . 共A9兲⌸0−1共k = 0;兲k2⌸0共k;兲
冏
k=0 = − lf 2冉
1 − 1 共2rlf兲2 + 1 2rlfarctan共2rlf兲冊
. 共A10兲 The solution to the integral冕
0 ⬁ dr Im兵G共r;兲rG*共r;兲其r2 =−sgn共兲 42lf冕
0 ⬁ drr2Im兵关ir+ 1/共2lf兲兴 ⫻K0关共− ir+ 1/共2lf兲兲r兴K1关共ir+ 1/共2lf兲兲r兴其 =−sgn共兲 42lf 2 Im再
2 共1 + i2rlf兲2 共1 − i2rlf兲4 ⫻F冉
2,2;3;1 −共1 + i2rlf兲 2 共1 − i2rlf兲2冊
冎
, 共A11兲can be found from Ref. 关27兴. However, the proper solution 共on the right Riemann sheet兲 needs to be chosen in order to simplify the hypergeometric series F. One can check numeri-cally that
冕
0 ⬁ dr Im兵G共r;兲rG*共r;兲其r2 =−sgn共兲 42lf 2 arctan共2rlf兲 ⫻冉
1 − 1 共2rlf兲2 + 1 2rlfarctan共2rlf兲冊
. 共A12兲 Hence, C is given by C =sgn共兲 42 arctan共2rlf兲. 共A13兲The intensity is proportional to the propagator ⌸共r;兲, expressed in terms of⌸0共k;兲 by the Fourier transform 共46兲. As only the gradient of the intensity is a well-defined prop-erty, we calculate ⌸共r;兲 = − rˆ
冕
0 ⬁ dk 2 k2J1共kr兲 ⌸0−1共k;兲 − ⌸0−1共k = 0;兲 . 共A14兲 This contains both the coherent and the scattered intensity. As the coherent intensity is known, we focus on the scattered intensity by calculating ⌸sc共r;兲 = − rˆ冕
0 ⬁ dk 2J1共kr兲k 2⌸ sc共k;兲, 共A15兲 with ⌸sc共k;兲 = ⌸0 −1共k = 0;兲⌸ 0共k;兲 ⌸0−1共k;兲 − ⌸0−1共k = 0;兲 . 共A16兲Equation共A15兲 is the integral to calculate numerically when we need to calculate the gradient of the multiply scattered intensity. ⌸sc共k;兲 is a monotonically decaying function with a maximum at k = 0, that vanishes as k→⬁. As the Bessel function is also decaying with r, we know that for
r / lfⰇ1 ⌸td共r;兲 = − rˆ ⌸sc共k = 0;兲 2r 共A17兲 and ⌸td共r兲 = − rˆ arctan共2rlf兲 22 rlf 1 rg −1共2 rlf兲, 共A18兲 where g共2rlf兲 = 1 − 1 共2rlf兲2 + 1 2rlfarctan共2rlf兲 . 共A19兲 td stands for “totally diffusive.”
关1兴 Proc. SPIE 1888 共1993兲, special issue on optical tomography, edited by B. Chance and R. R.Alfrano; S. R. Arridge, Inverse Probl. Eng. 15, R41共1999兲.
关2兴 A. K. Fung, Microwave Scattering and Emission Models and
Their Applications共Artech House, Boston, 1994兲.
关3兴 F. J. Fry, Ultrasound: Its Applications in Medicine and Biology 共Elsevier, Amsterdam, 1978兲; A. Migliori and J. L. Sarrao,
Resonant Ultrasound Spectroscopy, Applications to Physics, Materials Measurement and Nondestructive Evaluation
共Wiley, New York, 1997兲.
关4兴 R. Snieder, in Diffuse Waves in Complex Media, edited by J. P. Fouque共Kluwer, Dordrecht, 1999兲.
关5兴 E. Larose, A. Derode, M. Campillo, and M. Fink, J. Appl.
Phys. 95, 8393共2004兲; A. E. Malcolm, J. A. Scales, and B. A. van Tiggelen, Phys. Rev. E 70, 015601共R兲 共2004兲.
关6兴 P. Sheng, Introduction to Wave Scattering, Localization and
Mesoscopic Phenomena共Academic Press, New York, 1995兲.
关7兴 D. S. Wiersma, P. Bartolini, A. Lagendijk, and R. Righini, Nature共London兲 390, 671 共1997兲.
关8兴 M. Haney and R. Snieder, Phys. Rev. Lett. 91, 093902 共2003兲; S. E. Skipetrov and B. A. van Tiggelen, ibid. 92, 113901 共2004兲; S. K. Cheung, X. Zhang, Z. Q. Zhang, A. A. Cha-banov, and A. Z. Genack, ibid. 92, 173902共2004兲; E. Larose, L. Margerin, B. A. van Tiggelen, and M. Campillo, ibid. 93, 048501共2004兲.
共1990兲; R. H. J. Kop, P. de Vries, R. Sprik, and A. Lagendijk,
ibid. 79, 4369共1997兲; P.-A. Lemieux, M. U. Vera, and D. J.
Durian, Phys. Rev. E 57, 4498 共1998兲; Z. Q. Zhang, I. P. Jones, H. P. Schriemer, J. H. Page, D. A. Weitz, and P. Sheng,
ibid. 60, 4843共1999兲; X. Zhang and Z. Q. Zhang, ibid. 66,
016612共2002兲.
关10兴 M. Fink, Phys. Today 50, 34 共1997兲.
关11兴 P. de Vries, D. V. van Coevorden, and A. Lagendijk, Rev. Mod. Phys. 70, 447共1998兲.
关12兴 M. C. W. van Rossum and Th. M. Nieuwenhuizen, Rev. Mod. Phys. 71, 313共1999兲.
关13兴 G. E. W. Bauer, M. S. Ferreira, and C. P. A. Wapenaar, Phys. Rev. Lett. 87, 113902共2001兲.
关14兴 H. Sato and M. C. Fehler, Seismic Wave Propagation and
Scat-tering in the Heterogeneous Earth共Springer, New York, 1998兲.
关15兴 A. Tourin, M. Fink, and A. Derode, Waves Random Media 10, R31共2000兲.
关16兴 A. Lagendijk and B. A. van Tiggelen, Phys. Rep. 270, 143 共1996兲.
关17兴 R. Hennino, N. Trégourès, N. M. Shapiro, L. Margerin, M. Campillo, B. A. van Tiggelen, and R. L. Weaver, Phys. Rev. Lett. 86, 3447共2001兲; U. Wegler and B.-G. Lühr, Geophys. J.
Int. 145, 579共2001兲.
关18兴 J. de Rosny, A. Tourin, A. Derode, B. van Tiggelen, and M. Fink, Phys. Rev. E 70, 046601共2004兲.
关19兴 A. Tourin, A. Derode, P. Roux, B. A. van Tiggelen, and M. Fink, Phys. Rev. Lett. 79, 3637共1997兲; N. P. Trégourès and B. A. van Tiggelen, Phys. Rev. E 66, 036601共2002兲; G. Samel-sohn and V. Freilikher, ibid. 70, 046612共2004兲.
关20兴 R. Sprik and A. Tourin, Ultrasonics 42, 775 共2004兲.
关21兴 A. Messiah, Quantum Mechanics 共North-Holland, Amsterdam, 1962兲, Vol. II.
关22兴 E. Merzbacher, Quantum Mechanics 共Wiley, New York, 1961兲. 关23兴 P. Exner and P. Šeba, Phys. Lett. A 222, 1 共1996兲.
关24兴 N. G. Einspruch, E. J. Witterholt, and R. T. Truell, J. Appl. Phys. 31, 806共1960兲.
关25兴 A. Derode, A. Tourin, and M. Fink, Phys. Rev. E 64, 036605 共2001兲.
关26兴 M. S. Ferreira and G. E. W. Bauer, Phys. Rev. E 65, 045604共R兲 共2002兲.
关27兴 I. S. Gradshteyn and I. M. Ryzhik, Table of Integrals, Series
and Products共Academic Press, London, 1965兲.
关28兴 M. Abramowitz and I. A. Stegun, Handbook of Mathematical
