Space-time domain velocity distributions in isotropic radiative transfer in two dimensions
SSpace-time domain velocity distributions in isotropic radiative transfer in twodimensions
Vincent Rossetto
Université Grenoble Alpes / CNRSLaboratoire de physique et modélisation des milieux condensésMaison des Magistères - CNRS, BP 16625, avenue des Martyrs,38042 Grenoble CEDEX, France
We compute the exact solutions of the radiative transfer equation in two dimensions for isotropicscattering. The intensity and the radiance are given in the space-time domain when the sourceis punctual and isotropic or unidirectional. These analytical results are compared to Monte-Carlosimulations in four particular situations.
PACS numbers: 42.25.Dd, 05.20.-y
Transport in disordered and random media is a widelyaddressed physical question which plays an importantrole in several domains of Physics. Motivated by the ki-netic theory of gases, the Boltzmann equation has beenstudied since more than a century. The radiative transferequation is a Boltzmann equation where speed is fixed. Itwas derived by Chandrasekhar [3] to study the radiationtransport in a scattering atmosphere. Although radiativetransfer is mostly used in three-dimensional systems, thetwo-dimensional radiative transfer is of interest in severaldomains, such as seismology, where surface waves carrymost of the energy.Some solutions of the two-dimensional radiative trans-fer equation are known analytically. For isotropic scat-tering and isotropic source, the energy distribution hasbeen found by Shang and Gao, and Sato [13, 14] andPaasschens improved these results by providing the ra-diance distribution [10]. Recent progress in numericaland analytical solutions have been made by Liemert andKienle [4, 5]. These numerical methods efficiently ex-tend to the three-dimensional case [6, 7]. The situationin two dimensions is more favorable to analytical resultsfor many reasons, geometrical and analytical; let us onlymention that the rotation group has a single parameterand that the Green’s function has an algebraic Fourier-Laplace transform.Let us write q ( r , t, θ ) the space-time density of energyflux at position r at the time t with direction angle θ . q ( r , t, · ) is called the radiance in the standard terminol-ogy in optics. If we integrate the radiance with respectto θ , we get the spatial distribution of energy flux at r and t . In systems with no absorption, the energy fluxdistribution integrated over space is constant and nor-malized to c(cid:96) in this work.The differential equation for q is ∂ t q ( r , t, θ ) + c ˆ u ( θ ) · ∇ q ( r , t, θ ) + c(cid:96) q ( r , t, θ ) = c(cid:96) (cid:90) S ϕ ( θ − θ (cid:48) ) q ( r , t, θ (cid:48) ) d θ (cid:48) . (1)The phase function ϕ is an even real valued function de-scribing the distribution of scattering angle. The case where all scattering angles are equally likely, ϕ ( θ ) =1 / π , is the isotropic case.To solve the equation (1) we introduce “unscattered”distributions that are the spatial distribution of proba-bility of particles that have not been scattered. Thesedistributions are distinguished by a subscript . “Scat-tered” distributions receive the same notations withoutthis subscript. The probability to meet a scatterer onits trajectory at a distance r from the source is e − r/(cid:96) ,where (cid:96) is the mean free path. The distribution of par-ticles starting from the origin at time t = 0 and movingwith speed c with angle θ that have not been scatteredat time t is G ( r , t, θ ) = c(cid:96) δ (2) ( r − ct ˆ u ( θ )) e − ct/(cid:96) . (2) G is the unscattered energy distribution from an unidi-rectional point source . δ (2) is a two-dimensional Diracdelta function. In the absence of scattering, the prop-agation angle θ is preserved, its distribution is a Diracdelta-function δ ( θ − θ ) . This defines the unscattered ra-diance distribution from an unidirectional point source G ( r , t, θ ; θ ) = G ( r , t, θ ) δ ( θ − θ ) . (3)We remark that G ( r , t, θ ; θ ) is invariant if one ex-changes θ and θ . As a consequence G ( r , t, θ ) is alsothe unscattered radiance from an isotropic source . Fi-nally, integrating G over the angle θ , we obtain the un-scattered energy distribution from an isotropic source g ( r , t ) = 12 π (cid:90) S G ( r , t, θ ) d θ = c(cid:96) δ ( r − ct )2 πr e − ct/(cid:96) . (4)The distributions defined by Equations (2), (3) and (4)constitute the building blocks for the multiple scatteringtheory presented in this paper. Notations and analytic transforms
We use the units c = 1 , (cid:96) = 1 . We denote by (cid:101) f ( k ) the spatial Fouriertransform of the radial function f ( r ) and (cid:98) f ( s ) the timeLaplace transform of f ( t ) . We will also use the Hankeltransform as defined in the appendix A.The Fourier-Laplace transform of f ( r, t ) is denotedby f ( k, s ) . The leading exponential factor e − t in g , a r X i v : . [ c ond - m a t . s t a t - m ec h ] J a n G and G results in the shift of the variable s by .All Fourier-Laplace transformed functions carry this shiftand are written only with their angular dependences, asin G ( θ ) ≡ G ( k , s − , θ ) . The Fourier-Laplace trans-forms of the unscattered distributions admit the followingexpressions g = 1 √ k + s , (5) G ( θ ) = ( s + i k · ˆ u ( θ )) − , (6) G ( θ ; θ ) = δ ( θ − θ ) G ( θ ) . (7) I. ANALYTICAL DERIVATION
The energy distribution g of the two-dimensionalisotropic radiative transfer has been provided by Gao &Shang [14], Sato [13] and Paasschens [10], the latter hav-ing also given the radiance solution G . This section isdedicated to the analytical computation of the scatteredradiance from an unidirectional point source G . On theway to this result we compute the scattered Green’s func-tions g and G . A. Solutions in the Fourier-Laplace domain
The radiative transfer equation (1) governing G in theFourier-Laplace domain rewrites for isotropic scattering ( s + i k · ˆ u ( θ )) G ( θ ) = (cid:101) G ( k , , θ )+ 12 π (cid:90) S G ( θ (cid:48) ) d θ (cid:48) . (8)The initial condition G ( r , , θ ) = δ (2) ( r ) enters into theequation (8) as (cid:101) G ( k , , θ ) = 1 . We multiply the equation(8) by G ( θ ) as given by (6) and we obtain G ( θ ) = G ( θ ) + G ( θ ) 12 π (cid:90) S G ( θ (cid:48) ) d θ (cid:48) = G ( θ ) + G ( θ ) g . (9)We notice that G will be known as soon as we know g .From the integration of equation (9) over θ we get, usingthe definition (4), the Green-Dyson relation g = g + g g (10)and deduce from the expression of g in equation (5) the scattered energy distribution from an isotropic source g = g − g = g + g + g + · · · = 1 √ s + k − . (11)Each order of the expansion corresponds to a given num-ber of scattering events the particle has experienced. From the expression (9) we find the scattered radiancedistribution from an isotropic source to be G ( θ ) = G ( θ )1 − g = G ( θ ) + G ( θ ) g + G ( θ ) g + · · · (12)in which we can se that only the last scattering eventdepends on the angle. Conversely, the scattered energydistribution from an unidirectional point source is givenby the same formula written as G ( θ ) = G ( θ ) + g G ( θ ) + g G ( θ ) + · · · in which direction is lost after the first scattering event.If the source is unidirectional, the scattered radiance dis-tribution from an unidirectional point source is thereforegiven by the relation G ( θ ; θ ) = G ( θ ; θ ) + G ( θ ) G ( θ ) . (13)We can now use the solutions (11), (12) and (13) to givethe expressions of this functions in the space-time do-main. B. Solutions in the space-time domain
To find the expression of g in the space-time domain weuse the simultaneous Hankel-Laplace inverse transform oforder zero (see the appendix) with the function (cid:98) f ( s ) = s/ ( s − (and thus f ( t ) = δ ( t ) + e t Θ( t ) ) and we obtain g ( r, t ) = c e − ct/(cid:96) π(cid:96)r δ ( ct − r ) + e − ct/(cid:96) + cT/(cid:96) π(cid:96) T Θ( ct − r ) , (14)the energy distribution from an isotropic source as al-ready found by Shang and Gao and by Sato. We haveused the notation T = (cid:112) t − r /c . Note that the si-multaneous inverse transform was performed thanks tothe fact that g is a function of g − . For t (cid:29) r , this so-lution approaches the Gaussian distribution of diffusion,with D = = c(cid:96)/ . To compute G we can use the spaceand time convolution defined by the Equation (9) whichyields the expression G ( r , t, θ ) = G ( r , t, θ ) + e − ct/(cid:96) π(cid:96) e cT/(cid:96) t − r · ˆ u ( θ ) /c Θ( ct − r ) . (15)The expression (15) was first derived by Paass-chens [10]. Our work extends his results to the radiancedistribution from an unidirectional point source, G . Tocompute G , we use the result (15) together with the re-lation (13) which defines a convolution in the space-timedomain. After integration with respect to the space co-ordinate we get G ( r , t, θ ; θ ) = G ( r , t, θ ; θ ) + G ( r , t, θ ; θ ) + Θ( t − r ) e − t π (cid:90) t − r t − r · ˆ u ( θ e √ t − r − τ ( t − r · ˆ u ( θ )) t − τ − r · ˆ u ( θ ) + τ ˆ u ( θ ) · ˆ u ( θ ) d τ, FIG. 1. Description of the geometric variables in the situa-tion where r and θ are fixed. The trajectory with a singlescattering event is drawn as a dashed line. Around the point r , the shaded region shows the angle interval in which G contributes to the intensity. where G is the single scattering contribution arisingfrom the convolution of G with itself (see the figure 1and the appendix B).If r = t ˆ u ( θ ) , the integral vanishes (all the energy iscontained in the ballistic term G ), otherwise we can per-form the change of variable t − r · ˆ u ( θ )) τ = t − r − y and we get for θ (cid:54) = θ G = G + G + e − ct/(cid:96) π t − r )1 − cos( θ − θ ) (cid:90) T y e y X + y d y, where X is defined by X = 1 (cid:96) (cid:115) ct − r · ˆ u ( θ )) ( ct − r · ˆ u ( θ ))1 − cos( θ − θ ) − c T . (16)We finally obtain our main result for isotropic scattering(with θ (cid:54) = θ ) G ( r , t, θ ; θ ) = G ( r , t, θ ; θ ) + G ( r , t, θ ; θ ) + c e − ct/(cid:96) π(cid:96) ct − r )1 − cos( θ − θ ) Re (cid:104) E (i X ) e i X − E (cid:16) i X − c(cid:96) T (cid:17) e i X (cid:105) . (17)The function E n is the n th order exponential integral function as defined in [1, chap. 5]. The expression (17) has beenobtained using the antiderivative 5.1.44 in this reference. In the case where θ = θ , the integral gives the result G ( r , t, θ ; θ ) = G ( r , t, θ ) + Θ( ct − r ) c e − ct/(cid:96) π(cid:96) cT /(cid:96) −
1) e cT/(cid:96) ( ct − r · ˆ u ( θ )) . (18)The term G is the unscattered contribution while thesecond term is a scattered contribution of second order(at least two scattering events have occured). There areno single scattering contributions from G in (18). C. Steady-state solutions
The time-dependent scattered solutions measured at agiven point r exhibit a variety of behaviours that can beexploited when using pulse sources. However, some ex-perimental setups may require the use of a steady source.Hence, we discuss here the steady-state solutions of theradiative transfer equation in two dimensions. We haveto first remark that the large time regime is diffusive andas Brownian motion in two dimensions is recurrent, asteady source would yield a diverging energy density astime goes to infinity. However, in the presence of an ab-sorption rate µ > , all unscattered and scattered Green’sfunctions get a leading regularizing factor e − µt . Such aconstant rate could come from energy dissipation underanother form (like, typically, heat) or account for lossesinto the third dimension. Since the dimension two is thecritical dimension for Brownian recurrence, we expect the steady-state distribution to diverge logarithmically as µ or r goes to zero.In the presence of absorption, the steady-state coun-terpart f ss ( r ) of a Green’s function f ( r, t ) is well definedand we have f ss ( r ) = (cid:90) ∞ e − µt f ( r, t ) d t = (cid:98) f ( r, µ ) . (19)It could also be obtained as the inverse Fourier transformof (cid:101) f ss ( k ) = f ( k, µ ) . Denoting by α = (cid:96) − + µ/c thetotal extinction rate, the unscattered energy distributionis g ss0 ( r ) = e − αr / (2 π(cid:96)r ) and we show in the appendix Cthat g ss = g ss0 + g ss ∞ where g ss ∞ ( r ) = 12 π ∞ (cid:88) n =0 ( − κr ) n n ! E n +1 (cid:16) µr c (cid:17) . (20)with κ = (cid:96) − + µ/ c . The expression (20) is exact andis convenient for small r expansion, where it convergesquickly. Using the steepest descent method, we obtainan approximation of g ss ( r ) for large r ( r (cid:29) c/µ ) as g ss ∞ ( r ) ≈ r →∞ π ) / (cid:96) e − r √ κµ/c (cid:112) r (2 κµ/c ) / . FIG. 2. Radiance for fixed θ = 0 , t = 1 . and r = (0 , (thus θ r = π ) as a function of the propagation angle θ . The Monte-Carlo simulations have been performed until trajectorieswith at least two scattering events are found satisfying thefixed conditions for t and r within ∆ t = 0 . and ∆ r = 0 . .A total of . × trajectories have been computed. For large µ ( µ (cid:29) c/(cid:96) ) we find g ss ∞ ( r ) ≈ K ( µr ) / (2 π(cid:96) ) . Inboth cases, we observe a slower energy decay away fromthe source than for the unscattered energy distribution.The unscattered radiance distributions are propor-tional to g ss0 . We easily find G ss0 ( r , θ ) = g ss0 ( r ) δ ( θ − θ r ) and G ss0 ( r , θ ; θ ) = G ss0 ( r , θ ) δ ( θ − θ ) . Since the equation(9) states that G = G + G g we have (cid:101) G ss = (cid:101) G ss0 + (cid:101) G ss0 (cid:101) g ss .The steady-state distributions have the same convolutionrelations as the time dependent ones. The distribution G ss cannot be computed exactly, but we should remarkthat near the source, the lowest order of scattering dom-inates the distribution. The unscattered term, propor-tional to g ss0 , dominates everywhere it is not equal tozero. A single scattering contribution appears in G ss , wedenote it by G ss1 . We can therefore decompose G ss into G ss = G ss0 + G ss1 + G ss ∞ . On the figure 1, the shaded regioncorresponds to the geometric configuration where the ra-diance distribution from an unidirectional point sourcehas a contribution from single scattering. If single scat-tering does not contribute, the main contribution is fromdouble scattering, G ss2 (we do not provide an expressionfor this contribution). The distribution G ss decomposesinto G ss0 + G ss1 + G ss2 + G ss ∞ . The distributions G ss1 and G ss1 are given in the appendix B. II. NUMERICAL SIMULATIONS
We compare the solution (17) to statistics obtainedfrom a Monte-Carlo simulation of the two dimensionalisotropic Boltzmann equation. We start random walksfrom the origin at t = 0 with propagation angle θ = 0 .The step length x is distributed exponentially accordingto the probability distribution p ( x ) = e − x . After eachstep, we chose a random angle ≤ θ < π for the prop-agation.The figure 1 shows a situation where θ r = π usedin the Monte-Carlo simulations. If the random walker FIG. 3. Radiance in the direction orthogonal to the initialpropagation direction as a function of time at the position r = (1 , . Single scattering contributions have been removed.The Monte-Carlo simulations have been performed until trajectories with at least two scattering events are found sat-isfying the fixed conditions for r and θ within ∆ r = 0 . and ∆ θ = π . A total of . × trajectories have beencomputed. approaches the “target” position r by a distance lessthan ∆ r between the times t − ∆ t/ and t + ∆ t/ ,we store the value taken by θ during the correspond-ing step. The statistics of θ follow the distribution G ( r , t, θ, t × π (∆ r ) . The distribution exhibits apeak at θ = 2 tan − ( tr ) , corresponding to the single scat-tering trajectory. The results of these simulations is dis-played in the figure 2, they compare the distributions ofthe propagation angle θ for fixed position r and time t ob-tained by Monte-Carlo simulations to the predicted for-mula (17).The figure 3 displays the radiance at angle θ = π/ as a function of time at the fixed point r = (1 , . Ifthe random walker approaches the “target” position r bya distance less than ∆ r and the propagation angle θ (cid:48) issuch that cos( θ (cid:48) − θ ) > cos(∆ θ ) , we store the value of t corresponding to the closest point along the matchingstep.The figures 4 and 5 show the radiance at angle θ = θ = 0 , the pathological case where Equation (17) has tobe replaced by (B1), as a function of time at the fixedpoints r = (1 , and r = (1 , respectively. The methodis the same as explained for figure 3.In the figures, the normalization of the numerical dis-tributions has been adjusted, no other parameters havebeen tuned. III. OUTLOOK
We have computed the exact solutions of the radiativetransfer equation in two dimensions with angular reso-lution both at the source and the receiver. The time-dependent solutions are useful for the signal analysiswhen the source is modulated. The steady-state solu-tions could only be estimated: The approximations wehave obtained are asymptotically close to the exact solu-
FIG. 4. Radiance in the initial direction of propagation asa function of time at the position r = (1 , . Non-scatteredcontributions have been removed. The Monte-Carlo simula-tions have been performed until trajectories are foundsatisfying the fixed conditions for r and θ within ∆ r = 0 . and ∆ θ = π . A total of . × trajectories have beencomputed.FIG. 5. Radiance in the initial direction of propagation asa function of time at the position r = (1 , . Non-scatteredcontributions have been removed. The Monte-Carlo simula-tions have been performed until trajectories are foundsatisfying the fixed conditions for r and θ within ∆ r = 0 . and ∆ θ = π . A total of . × trajectories have beencomputed. tions when the absorption is strong or when the distancefrom the source is large. When absorption is low andthe distance from the source is small, the single scatter-ing contribution grows logarithmically and dominates theradiance.We now briefly discuss the applications of these re-sults. The angular resolution of the theoretical radiancewill be useful for analysing data collected with full orpartial angular dependences. Angular dependences areeasily acessible in optics: Some light sources, such aslaser beams, are inherently unidirectional, and collimatedreceivers can be used to measure the radiance. In acous-tics, the so-called beam-forming methods use the signalsrecorded by an array of aligned receivers to select thesound from an incoming direction. This technique alsoworks with an arrays of sources to produce a unidirec- tional source of sound. These beam-forming techniquesare also frequently used in geophysics as well, with seis-mic waves. In these fields, using the angular dependencesof scattered waves would represent a substantial increaseof the available amount of data. We expect that such anincrease will help improving imaging methods. Appendix A: Simultaneous Hankel-Laplacetransform
We give here a proof of the simultaneous Hankel-Laplace transformation formula of arbitrary order n .This transformation is more general than the case n = 0 used in the main text. The two-dimensional Fouriertransform of a function f ( r ) = (cid:80) n ∈ Z f n ( r )e i nθ r is bydefinition (cid:80) n i n e i nθ k ˇ f n ( k ) , with ˇ f n the Hankel tranformof | n | th order of f n .The n th order Hankel transform of δ ( r − x ) / πr is given by J n ( kx ) . Replacing x by √ t − u andmultiplying by ( t − u ) n/ ( t + u ) − n/ we recognizethe Laplace transform 29.3.97 in [1] which is equalto exp (cid:0) − u √ s + k (cid:1) ξ n with ξ = k/ ( s + √ s + k ) .The expression exp (cid:0) − u √ s + k (cid:1) ξ n is therefore the n -Hankel-Laplace transform of ( t − u ) n/ ( t + u ) − n/ δ ( r −√ t − u ) / πr . Multiplying both these expressions byan arbitrary function f ( u ) and integrating from u = 0 toinfinity we find that √ s + k (cid:18) ks + √ s + k (cid:19) n (cid:98) f (cid:16)(cid:112) s + k (cid:17) (A1)is the n th order Hankel-Laplace transform of π Θ( t − r ) √ t − r (cid:18) rt + √ t − r (cid:19) n f (cid:16)(cid:112) t − r (cid:17) . (A2)We call the transformation between (A1) and (A2) asimultaneous double transform. Simultaneous doubletransforms have been introduced recently in Ref. [12].It should be noticed that the simultaneous inversion foranisotropic scattering is possible because the expansionof G is a rational function of g − = √ k + s . Appendix B: The single scattering functions
The term G in equation (17) is purely geometric. Itscomputation is straighforward and yields G ( r , t, θ ; θ ) = e − ct/(cid:96) δ (2) (cid:2)(cid:0) r − ct ˆ u ( θ ) (cid:1) × (cid:0) ˆ u ( θ ) − ˆ u ( θ ) (cid:1)(cid:3) Θ (cid:2)(cid:0) r − ct ˆ u ( θ ) (cid:1) · (cid:0) ˆ u ( θ ) − ˆ u ( θ ) (cid:1)(cid:3) Θ (cid:2) ct (cid:0) − cos( θ − θ ) (cid:1) − r · (cid:0) ˆ u ( θ ) − ˆ u ( θ ) (cid:1)(cid:3) . (B1)In the steady state regime, it reduces to G ss1 ( r , θ ; θ ) =Θ [cos( θ r − θ ) − cos( θ r − θ ) cos( θ − θ )]4 π (cid:96) | sin( θ − θ ) | exp (cid:18) − αr (cid:12)(cid:12)(cid:12)(cid:12) sin( θ r − θ ) − sin( θ r − θ )sin( θ − θ ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:19) . (B2)In the steady-state regime, the contribution of singlescattering to the radiance distribution from an isotropicsource arises from the convolution of g ss0 and G ss0 , it isequal to G ss1 ( r , θ ) = 14 π (cid:96) exp (cid:16) − αr cos( θ r − θ ) (cid:17)(cid:17) E (cid:16) αr (cid:16) − cos( θ r − θ ) (cid:17)(cid:17) . (B3) Appendix C: Steady-state approximations
The steady-states solutions and approximations arebased on the following expressions. We perform thechange of variable t = r c (cid:0) x + x (cid:1) in the time integral(19) of g given by equation (14) and we get g ss ∞ ( r ) = 12 π (cid:90) ∞ exp (cid:16) − µr x − κrx (cid:17) d xx . We obtain the equation (20) by expanding the exponen-tial e − κr/x into the series (cid:80) n ( − κr/x ) n n ! . The large r ap-proximation is obtained by second order polynomial ex-pansion of µr x + κrx around the minimum at x = (cid:112) κ/µ . [1] Milton Abramowitz and Irene A. Stegun, Handbook ofmathematical functions , 10th ed., Dover, 1972.[2] Jérémie Boulanger, Nicolas Le Bihan, and Vincent Ros-setto,
Stochastic description of geometric phase for polar-ized waves in random media , J. Phys. A. Math. Theor. (2013), 035203.[3] Subrahmanyan Chandrasekhar, Radiative transfer ,Dover, 1960.[4] André Liemert and Alwin Kienle,
Radiative transfer intwo dimensional infinitely extended scattering media , J.Phys. A. : Math. Theor. (2011), 505205.[5] , Green functions for the two dimensional radia-tive transfer equation in bounded media , J. Phys. A.:Math. Theor. (2012), 175201.[6] , Light transport in three-dimensional semi-infinite scattering media , J. Opt. Soc. Am. (2012),1475.[7] , Explicit solutions of the radiative transport equa-tion in the P3 approximations , Med. Phys. (2014),111916.[8] Manabu Machida, George Panasyuk, John Schotland,and Vadim Markel, The Green’s function for the radia-tive transport equation in the slab geometry , J. Phys. A.: Math. Theor. (2010), 065402.[9] Nori Nakata, Pierre Boué, Florent Brenguier, PhilippeRoux, Valérie Ferrazzini, and Michel Campillo, Body andsurface wave reconstruction from seismic noise correla-tions between arrays at Piton de la Fournaise volcano ,Geophys Res. Lett. (2016), 1047–1054.[10] J. C. J. Paasschens, Solution ot the time-dependent Boltz-mann equation , Phys. Rev. E (1997), 1135–1141.[11] Vincent Rossetto, A general framework for multiplescattering of polarized waves including anisotropies andBerry phase , Phys. Rev. E (2009), 056605.[12] Vincent Rossetto, Simultaneous double transforma-tions for function depending on space and time ,arXiv:1311.3140, 2013.[13] Haruo Sato,
Energy transport in one- and two-dimensional scattering media: Analytic solutions of themultiple isotropic scattering model , Geo. Phys. J. Int. (1993), 487–494.[14] T. Shang and L. Gao,
Transportation theory of multiplescattering and its application to seismic coda waves ofimpulse source , Scientia Sinica B31