Abstract
The radiative transport equation (RTE) is used widely to describe the propagation of multiply scattered light in disordered media. In this tutorial, we present two derivations of the RTE for scalar wave fields. The first derivation is based on diagrammatic perturbation theory, while the second stems from an asymptotic multiscale expansion. Although the two approaches are quite distinct mathematically, some common ground can be found and is discussed.
© 2015 Optical Society of America
1. INTRODUCTION
The propagation of multiply scattered light in complex media, including clouds, colloidal suspensions, and biological tissues, is often described by the radiative transport equation (RTE) [1]. The RTE has its genesis in Boltzmann’s kinetic theory of gases and is now employed in a variety of disciplines, ranging from atmospheric physics [2] to biomedical imaging [3]. There are also parallel applications to neutron transport [4] and mesoscopic physics [5] that have served to enrich and stimulate research in the theory of radiative transport.
In radiative transport theory, the propagation of light is formulated as a conservation law that accounts for gains and losses of electromagnetic energy because of scattering and absorption. The fundamental quantity of interest is the specific intensity , defined as the intensity of light at the position in the direction . The specific intensity obeys the RTE
which we have written in its time-independent form for a non-absorbing medium. Here is the scattering coefficient and is the phase function, which describes the scattering from the direction into the direction . The phase function is taken to obey the reciprocity relation and the normalization condition for all . We note that in an absorbing medium, the coefficient on the left-hand side of the RTE is replaced by the quantity , where is the absorption coefficient.The mathematical description of light propagation varies according to the scale of interest. The Maxwell equations hold on the scale of the optical wavelength , which we will refer to as the microscopic scale. The RTE is valid on the mesoscopic scale, which is characterized by the scattering length . We will be principally concerned with systems in which there is a separation of scales: , where is the distance of propagation. The fundamental question that concerns us is how to understand the origin of this scale separation. That is, how does a microscopic theory of waves give rise to a macroscopic theory that is based on a picture of particle transport? This question is discussed by Mandel and Wolf in their monumental treatise Optical Coherence and Quantum Optics; they note the following [6]:
In spite of the extensive use of the theory of radiative energy transfer, no satisfactory derivation of its basic equation, from electromagnetic theory or even scalar wave theory has been obtained up to now, except in some special cases.
In this tutorial, we provide a partial answer to the above question. We begin by reviewing the theory of wave propagation in random media. We emphasize that randomness is employed in the subjective sense, namely as a proxy for information about the medium. That is, the medium is modeled as a realization of a stochastic process with known statistics. We then proceed to derive the RTE by two complementary approaches. The first derivation is based on diagrammatic perturbation theory, and its elements are relatively well known to physicists [7–10]. The second derivation stems from multiscale asymptotic analysis, a tool that is widely used in applied mathematics [11,12]. Both derivations make use of the Wigner transform of the field in an essential way, but are otherwise quite distinct mathematically. We note that the role of the Wigner transform was anticipated by Wolf nearly 40 years ago [13].
The remainder of this paper is organized as follows. In Section 2, we introduce the Wigner transform and the theory of scalar waves in inhomogeneous media. The derivations of the RTE using diagrammatic and multiscale approaches are then exposed in Sections 3 and 4, respectively. Finally, we make a comparison of the two approaches in Section 5.
2. WIGNER TRANSFORM
We consider the propagation of a monochromatic electromagnetic field in an inhomogeneous dielectric medium, that is assumed to be nonmagnetic. For simplicity, we ignore the effects of polarization and consider a scalar field which obeys the time-independent wave equation
where is the free-space wavenumber, is the source, and is the dielectric susceptibility. For the remainder of this paper, we take to be purely real, which corresponds to the case of a nonabsorbing medium. The conservation of energy is governed by the relation [14] where the energy current is defined asNote that plays the role of the Poynting vector in the scalar theory.
Although the conservation law (3) provides some indication of how the intensity of the field is distributed in space, it does not describe how the intensity propagates. To obtain a local conservation law for the intensity that is resolved over both position and direction, we introduce the Wigner transform of the field:
The Wigner transform, which was first introduced in quantum mechanics to describe the local conservation of probability, will be seen to be related to the specific intensity [15]. The Wigner transform has several important properties. It is real-valued and is related to the intensity by
In addition, the Wigner transform is related to the energy current by
To understand the physical meaning of the Wigner transform and its relation to the specific intensity, it is instructive to consider the case of wave propagation in vacuum. The solution of the wave equation in this case is a plane wave of the form
where is the amplitude of the wave and is its direction. The Wigner transform is then given by and the intensity and current density areIf we identify the Wigner transform with the specific intensity, the above results are then consistent with interpretation of the specific intensity as the amount of energy flowing in direction at point . However, this identification is problematic. The Wigner transform, although real-valued, is not necessarily nonnegative. In contrast, on physical grounds, the specific intensity cannot take on negative values. Evidently, there are two ways to resolve the inconsistency, neither of which is entirely satisfactory: (1) only the high-frequency asymptotic behavior of the Wigner transform is of interest, as discussed in Section 4; (2) the Wigner transform is not directly measurable. However, the intensity and energy current are observable and can be obtained from moments of the Wigner transform. Alternatively, we note that each term in the collision expansion for the RTE is nonnegative and, thus, the solution to the RTE must be generally nonnegative [4].
3. DIAGRAMMATIC APPROACH
In this section, we derive the RTE from diagrammatic perturbation theory. We proceed as follows. First, we use Green’s function methods to express the field in the form of a multiple scattering series. We then use this expansion to obtain an integral equation obeyed by the average Green’s function, which has a convenient diagrammatic form. By following a similar procedure, we average the product of two Green’s functions, which can then be related to the Wigner transform of the field. Finally, we make use of the separation between microscopic and macroscopic scales in the form of the Kubo limit to obtain the RTE.
A. Multiple Scattering Expansion
We begin by considering the problem of wave propagation in a random medium. We assume that the field obeys the wave equation (2), and we take the susceptibility to be a Gaussian random field with correlations
where denotes statistical averaging. We assume that the medium is statistically homogeneous and isotropic. That is, the correlation function depends only upon the quantity . A case of particular interest is that of white noise disorder, in which , where is constant.The solution to the wave equation (2) is given by
Here the Green’s function obeys
Following standard procedures [14], we find that satisfies the integral equation
where the unperturbed free-space Green’s function is given byEquation (15) is a self-consistent equation for . It can be iterated as follows:
We note that Eq. (17) can be interpreted as a multiple scattering series, with each term representing a successively higher order of scattering [14].
B. Dyson Equation
We now turn to the problem of calculating the average Green’s function . As we will see, it will prove useful to represent Eq. (17) in diagrammatic form:
The diagrammatic rules are as follows and are summarized in Table 1:
- • A straight line with a left pointing arrow corresponds to a factor of the Green’s function .
- • A vertex corresponds to a factor of .
- • Integration is carried out over all coordinates corresponding to internal vertices.
To average over disorder, we make use of the following properties of Gaussian random fields. All odd moments of vanish. Even moments can be expressed as products of second moments according to Wick’s theorem:
where the sum is over all distinct pairs of indices. As an example, we note thatAveraging Eq. (17) now consists of pairing the vertices in all possible ways and summing the result. This process is illustrated diagrammatically as follows:
Note that a curved line is associated with a factor of the correlation function . In Eq. (21), we see that there are diagrams which are disconnected in the sense that, by cutting an internal line, a diagram factorizes into two subdiagrams. A diagram that is not disconnected is said to be connected. We define the self-energy to be the sum of all connected diagrams with external lines amputated:
We note that the self-energy can be written as
where , and correspond to the first three self-energy diagrams and are given byBy simple algebraic manipulation of diagrams, we find that the average Green’s function can be expressed in terms of the self-energy as
Equation (27) is known as the Dyson equation. It is represented diagrammatically as
We note that the Dyson equation can be iterated to obtain the series
The Dyson equation can be solved by the Fourier transformation. This follows from the fact that and are invariant under translations, which is a consequence of statistical homogeneity. We then define the Fourier transform of by
which defines . Likewise, we define and byNote that
where is a small positive constant. Upon Fourier transforming the Dyson equation, we obtainSolving for we find that
Evidently, the problem of calculating the self-energy is equivalent to calculating the average Green’s function itself. As a result, various approximations are introduced. The simplest of these is known as the weak scattering approximation (for reasons that will become clear), in which only the first self-energy diagram is retained. Within the accuracy of this approximation, we find upon taking the Fourier transform of Eq. (24) that
Next, we make use of the identities
where denotes the Cauchy principal value. We thus obtain where we have used the fact that the medium is nonabsorbing, which means that is real-valued. Now, can be neglected by introducing a high-frequency cutoff in Eq. (39). Physically, the introduction of such a cutoff corresponds to introducing a minimum scale for the spatial variations of (the size of the smallest scatterer). Let us define the scattering length asWe note that is well defined because statistical isotropy forces to depend only on . Putting everything together and using Eq. (40), we obtain the expression for the Fourier transformed average Green’s function:
whereBy performing an inverse Fourier transform, we find that in the weak scattering limit , the average Green’s function is given by
We conclude that the average field decays exponentially on the scale of the scattering length. We also note that the higher order self-energy diagrams contain two copies of the correlation function and can thus be estimated to be of order .
C. Bethe–Salpeter Equation
We now consider the problem of calculating correlations of the field. Our primary tool is the Bethe–Salpeter equation, which is an integral equation that is obeyed by the average of a product of Green’s functions. We begin by noting that the two-point correlation function of the field is given by
where we have used Eq. (13) and have assumed that the source is deterministic. We are thus led to develop a diagrammatic expansion for the quantity , which is given byNote that the retarded Green’s function is denoted by a left pointing arrow and the advanced Green’s function by a right pointing arrow. Averaging over the disorder as in Section 3.B, we obtain the following diagrammatic expansion for :
As in Eq. (21), the above diagrams fall into connected and disconnected types. In a manner similar to the construction of the self-energy, we define the irreducible vertex as the sum of all connected diagrams with external legs amputated:
The irreducible vertex can be written as the sum
whereBy algebraic manipulation of diagrams, we find that obeys an integral equation involving the irreducible vertex of the form
Equation (53), which is known as the Bethe–Salpeter equation, is represented diagrammatically as follows:
As is the case of self-energy, the irreducible vertex must be calculated approximately. In the weak scattering limit , also known as the ladder approximation, only the first diagram, corresponding to , is retained. Using this approximation, the Bethe–Salpeter equation becomes
The above equation can be iterated to produce the series of so-called ladder diagrams, which are of the form
D. Average Wigner Transform
The average of the Wigner transform is given by
where the statistical average on the left-hand side is not indicated explicitly. By making use of the Bethe–Salpeter equation within the ladder approximation, it is shown in Appendix A that the Fourier transform of the Wigner transform is given by where denotes the wave vector conjugate to the position . Next, we make use of the identities and Eq. (35) to obtain the relationHere and are defined as
It follows from the above and Eq. (45), that Eq. (58) becomes
Finally, upon inverse Fourier transforming with respect to , we obtain the integral equation obeyed by the average Wigner transform within the ladder approximation:
E. Kubo Limit
We now consider Eq. (65) in the Kubo limit . To this end, we replace and by and , respectively. The purpose of the Kubo limit is to separate spatial scales. That is, we assume the average Wigner transform varies slowly on the length scale over which the medium fluctuates. To make further progress, we apply Eqs. (37) and (42) in the weak scattering limit to obtain
We also use Eqs. (40), (41), and (62) to obtain
where the scattering coefficient . Equation (65) thus becomesLet us identify the specific intensity and phase function by
We find that Eq. (68) becomes
where the source term is defined byEvidently, we have just derived the RTE. We note that the coefficient and the phase function are given in terms of correlations of the medium. In the case of white noise disorder, the correlation function , where is constant. We find that
which corresponds to isotropic scattering.In addition to continuous media, radiative transport theory may be applied to random media composed of discrete scatterers. To this end, we consider a medium consisting of identical single-particle scatterers with susceptibility
where is the single-particle susceptibility and are the positions of the particles. In addition, we make the following assumptions: (1) are iid and uniformly distributed (2) is short ranged, and (3) we work in the high density weakly-scattering limit and impose the scaling and such that , where is the number density of the particles and . It can then be seen that the correlation function is given by and where is the differential scattering cross section and is the scattering cross section of a single particle.4. MULTISCALE APPROACH
In this section, we derive the RTE from the asymptotics of high-frequency waves in random media. We begin by rescaling the wave equation to allow for the separation of microscopic and macroscopic scales. We then introduce the scaled Wigner transform and obtain the corresponding Liouville equation. Finally, we average the Wigner transform by making use of a multiscale asymptotic expansion, which ultimately leads to the RTE.
A. High-Frequency Limit
We consider the wave equation in vacuum:
The solutions of Eq. (76) oscillate on the scale of the optical wavelength . On the other hand, we are interested in the behavior of the solutions on the macroscopic scale . To this end, we introduce a small parameter and rescale the position by . Equation (76) thus becomes
where . We will refer to as the high-frequency limit. The solution of Eq. (77) is of the form which is a plane wave with amplitude propagating in the direction . The corresponding Wigner transform is given by which is not well behaved in the limit. Thus, the Wigner transform must be modified to handle oscillatory functions. To solve this problem, we introduce the rescaled Wigner transform which is properly adapted to the high-frequency limit. Using this definition, we see that, for the plane-wave in Eq. (78), which is localized in frequency and independent of , consistent with Eq. (9).We now consider the case of a random medium in the high-frequency limit. We assume that the disorder is sufficiently weak so that the correlation function , defined in Eq. (12), is of the order . Thus, we rescale the susceptibility by and find that the wave equation (2) becomes
B. Liouville Equation
We now turn to the derivation of the Liouville equation, which is a conservation law for the Wigner transform. Let . Since is real-valued, it follows that satisfies the pair of wave equations
Substracting Eq. (84) from Eq. (83) yields
We now perform the change of variables
and make use of the identityEquation (85) thus becomes
Finally, we multiply Eq. (88) by , integrate with respect to , and make use of the definition Eq. (80) of the Wigner transform. We find, as shown in Appendix B, that obeys the Liouville equation
whereWe note that the Liouville equation is analogous to Eq. (65). The analogy is not perfect, however, since the Liouville equation is exact and Eq. (65) is valid within the ladder approximation.
C. Multiscale Expansion
We now consider the asymptotics of the Wigner transform in the high-frequency limit, which allows for the separation of microscopic and macroscopic scales. We will see that averaging over the microscopic scale leads directly to the RTE. Following standard procedures [16], we introduce a multiscale expansion for the Wigner transform of the form
where is a fast variable and is taken to be deterministic. We then regard and as independent variables and make the replacementThe Liouville equation (89) thus becomes
Inserting Eq. (91) into Eq. (93) and collecting terms of the same order in , we find that at
Equation (94) can be solved by Fourier transforms with the result
where and is a small positive regularization parameter that will eventually be set to zero.At we find that
More generally, the recursion relation
is obeyed for . The RTE may be derived by averaging Eq. (97) over realizations of the random medium. To proceed, we make the assumption , which closes the hierarchy Eq. (98) at . We note that this assumption cannot be mathematically justified, but is known to provide the correct limit as for the paraxial wave equation [17]. Equation (97) thus becomesThen, substituting the formula Eq. (95) for into Eq. (99), and using the relation
which follows from Eq. (12), we obtainNote that, in performing the indicated average, we have made use of the fact that . Appendix C presents the necessary details. Next, using the identities
and relating the specific intensity to by we find that Eq. (101) becomesFinally, recognizing the scattering coefficient Eq. (41) and the phase function Eq. (69), we recover the RTE
The above result agrees with the diagrammatic derivation of the RTE in Section 3.
5. DISCUSSION
In this tutorial, two different approaches to the derivation of the RTE have been outlined. The first derivation is based on the diagrammatic perturbation theory, and the second makes use of a multiscale asymptotic expansion. As previously explained, the Wigner transform plays a central role in both approaches, as do the assumptions of Gaussian disorder and statistical homogeneity. An additional common assumption is the separation of microscopic and macroscopic scales. Several points of departure should also be noted. In particular, it is not clear what the counterpart of the closure relation is in the diagrammatic perturbation theory. The same is true of the self-averaging of , the lowest order term in the multiscale expansion of the Wigner transform. Likewise, the analog of the ladder approximation in the multiscale expansion is not understood. Nevertheless, the fact that both approaches lead to the RTE is striking and suggests that it may be possible to produce a more refined derivation, including elements from both approaches.
APPENDIX A: DERIVATION OF EQ. (58)
Here we derive Eq. (58). Using Eq. (45) for the correlation function of the field and Eq. (55), the Bethe–Salpether equation in the ladder approximation, we obtain
Consider the Fourier transform of the average Wigner transform:
First, we evaluate :
Next, we evaluate :
Making the change of variable , we have
Finally, using
and putting everything together, we obtain which is Eq. (58).APPENDIX B: DERIVATION OF EQ. (89)
Here we derive the Liouville equation (89). The starting point is Eq. (88), which we repeat for convenience:
First, we multiply Eq. (B1) by and integrate over . Let us first deal with the left-hand side of the above:
Next, we consider the first term on the right-hand side of Eq. (B1):
The second term on the right-hand side of Eq. (B1) is obtained similarly.
Finally, putting everything together, we obtain the Liouville equation (89).
APPENDIX C: DERIVATION OF EQ. (101)
Here we derive Eq. (101). To proceed, we substitute Eq. (95) into the right-hand side of Eq. (99) and obtain
FUNDING INFORMATION
National Science Foundation (NSF) (DMR 1120923, DMS 1108969, DMS 1115574).
ACKNOWLEDGMENTS
J.C. Schotland is grateful to Guillaume Bal, Lenya Ryzhik, Sam Schotland, and Emil Wolf for valuable discussions.
REFERENCES
1. S. Chandrasekhar, Radiative Transfer (Dover, 1960).
2. A. Ishimaru, Wave Propagation and Scattering in Random Media (Academic, 1978).
3. S. Arridge and J. C. Schotland, “Optical tomography: forward and inverse problems,” Inverse Probl. 25, 123010 (2009). [CrossRef]
4. K. M. Case and P. F. Zweifel, Linear Transport Theory (Addison-Wesley, 1967).
5. E. Akkermans and G. Montambaux, Mesoscopic Physics of Electrons and Photons (Cambridge University, 2007).
6. L. Mandel and E. Wolf, Optical Coherence and Quantum Optics (Cambridge University, 1995).
7. Y. N. Barabanenkov and V. M. Finkelberg, “Radiation transport equation for correlated scatterers,” Sov. Phys. J. Exp. Theor. Phys. 26, 587 (1968).
8. D. Vollhardt and P. Wolfle, “Diagrammatic self-consistent treatment of the Anderson localization problem in dimensions,” Phys. Rev. B 22, 4666–4679 (1980). [CrossRef]
9. U. Frisch, Probabilistic Methods in Applied Mathematics, A. T. Barucha-Reid, ed. (Academic, 1968).
10. F. C. MacKintosh and S. John, “Diffusing-wave spectroscopy and multiple scattering of light in correlated random media,” Phys. Rev. B 40, 2383–2406 (1989). [CrossRef]
11. L. Ryzhik, G. Papanicolaou, and J. B. Keller, “Transport equations for elastic and other waves in random media,” Wave Motion 24, 327–370 (1996). [CrossRef]
12. G. Bal, T. Komorowski, and L. Ryzhik, “Kinetic limits for waves in a random medium,” Kinet. Relat. Models 3, 529–644 (2010). [CrossRef]
13. E. Wolf, “New theory of radiative transfer in free electromagnetic fields,” Phys. Rev. D 13, 869–886 (1976). [CrossRef]
14. M. Born and E. Wolf, Principles of Optics (Cambridge University, 1997).
15. E. Wigner, “On the quantum correction for thermodynamic equilibrium,” Phys. Rev. 40, 749–759 (1932). [CrossRef]
16. A. Nayfeh, Perturbation Methods (Wiley, 1973).
17. L. Erdos and H. T. Yau, “Linear Boltzmann equation as the weak coupling limit of a random Schrodinger equation,” Comm. Pure Appl. Math. 53, 667–735 (2000). [CrossRef]