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 $I(\mathbf{r},\widehat{\mathbf{s}})$, defined as the intensity of light at the position $\mathbf{r}$ in the direction $\widehat{\mathbf{s}}$. The specific intensity obeys the RTE

The mathematical description of light propagation varies according to the scale of interest. The Maxwell equations hold on the scale of the optical wavelength $\lambda$, which we will refer to as the microscopic scale. The RTE is valid on the mesoscopic scale, which is characterized by the scattering length ${\ell }_s=1/{\mu }_s$. We will be principally concerned with systems in which there is a separation of scales: $\lambda \ll {\ell }_s\ll L$, where $L$ 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 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 $U$ which obeys the time-independent wave equation

Note that $\mathbf{J}$ 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 $I={|U|}^2$ by

In addition, the Wigner transform is related to the energy current $\mathbf{J}$ 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

If 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 $\widehat{\mathbf{s}}$ at point $\mathbf{r}$. 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 $U$ obeys the wave equation (2), and we take the susceptibility $\eta$ to be a Gaussian random field with correlations

The solution to the wave equation (2) is given by

Here the Green’s function $G$ obeys

Following standard procedures [14], we find that $G$ satisfies the integral equation

Equation (15) is a self-consistent equation for $G$. 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 $〈G〉$. As we will see, it will prove useful to represent Eq. (17) in diagrammatic form:

(18)

The diagrammatic rules are as follows and are summarized in Table 1:

Table 1. Diagrammatic Rules

View Table

To average over disorder, we make use of the following properties of Gaussian random fields. All odd moments of $\eta$ vanish. Even moments can be expressed as products of second moments according to Wick’s theorem:

Averaging Eq. (17) now consists of pairing the vertices in all possible ways and summing the result. This process is illustrated diagrammatically as follows:

(21)

Note that a curved line is associated with a factor of the correlation function $C$. 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 $\mathrm{\Sigma }$ to be the sum of all connected diagrams with external lines amputated:

(22)

We note that the self-energy can be written as

By simple algebraic manipulation of diagrams, we find that the average Green’s function $〈G〉$ can be expressed in terms of the self-energy as

Equation (27) is known as the Dyson equation. It is represented diagrammatically as

(28)

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 $〈G〉$ and $\mathrm{\Sigma }$ are invariant under translations, which is a consequence of statistical homogeneity. We then define the Fourier transform of $〈G〉$ by

Note that

Solving for $〈G(\mathbf{k})〉$ 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 ${\mathrm{\Sigma }}_1$ 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

We note that ${\ell }_s$ is well defined because statistical isotropy forces $\tilde{C}(k_0(\widehat{\mathbf{s}}-{\widehat{\mathbf{s}}}^{\prime }))$ to depend only on $\widehat{\mathbf{s}}·{\widehat{\mathbf{s}}}^{\prime }$. Putting everything together and using Eq. (40), we obtain the expression for the Fourier transformed average Green’s function:

By performing an inverse Fourier transform, we find that in the weak scattering limit $k_0{\ell }_s\gg 1$, 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 $C$ and can thus be estimated to be of order $O(1/{\ell }_s^2)$.

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

(46)

Note that the retarded Green’s function $G_0$ is denoted by a left pointing arrow and the advanced Green’s function $G_0^{*}$ by a right pointing arrow. Averaging over the disorder as in Section 3.B, we obtain the following diagrammatic expansion for $〈G{G}^{*}〉$:

(47)

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 $\mathrm{\Gamma }$ as the sum of all connected diagrams with external legs amputated:

(48)

The irreducible vertex can be written as the sum

By algebraic manipulation of diagrams, we find that $〈G{G}^{*}〉$ 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:

(54)

As is the case of self-energy, the irreducible vertex must be calculated approximately. In the weak scattering limit $k_0{\ell }_s\gg 1$, also known as the ladder approximation, only the first diagram, corresponding to ${\mathrm{\Gamma }}_1$, 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

(56)

D. Average Wigner Transform

The average of the Wigner transform is given by

Here $\mathrm{\Delta }\mathrm{\Sigma }$ and $\mathrm{\Delta }G$ are defined as

It follows from the above and Eq. (45), that Eq. (58) becomes

Finally, upon inverse Fourier transforming $\tilde{W}(\mathbf{k},\mathbf{Q})$ with respect to $\mathbf{Q}$, 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 $\mathbf{Q}\to \mathbf{0}$. To this end, we replace $\mathrm{\Delta }G(\mathbf{Q},\mathbf{k})$ and $\mathrm{\Delta }\mathrm{\Sigma }(\mathbf{Q},\mathbf{k})$ by $\mathrm{\Delta }G(\mathbf{0},\mathbf{k})$ and $\mathrm{\Delta }\mathrm{\Sigma }(\mathbf{0},\mathbf{k})$, respectively. The purpose of the Kubo limit is to separate spatial scales. That is, we assume the average Wigner transform $W$ 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 $k_0{\ell }_s\gg 1$ to obtain

We also use Eqs. (40), (41), and (62) to obtain

Let us identify the specific intensity $I$ and phase function $A$ by

We find that Eq. (68) becomes

Evidently, we have just derived the RTE. We note that the coefficient ${\mu }_s$ and the phase function $A$ are given in terms of correlations of the medium. In the case of white noise disorder, the correlation function $C(\mathbf{r})=C_0\delta (\mathbf{r})$, where $C_0$ is constant. We find that

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

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 $\lambda =2\pi /k_0$. On the other hand, we are interested in the behavior of the solutions on the macroscopic scale $L\gg \lambda$. To this end, we introduce a small parameter $ϵ=1/(k_{0}L)$ and rescale the position $\mathbf{r}$ by $\mathbf{r}\to \mathbf{r}/ϵ$. Equation (76) thus becomes

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 $C$, defined in Eq. (12), is of the order $O(ϵ)$. Thus, we rescale the susceptibility $\eta$ by $\eta \to \sqrt{ϵ}\eta$ 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 ${\mathrm{\Phi }}_{ϵ}({\mathbf{r}}_1,{\mathbf{r}}_2)=U_{ϵ}({\mathbf{r}}_1)U_{ϵ}^{*}({\mathbf{r}}_2)$. Since $\eta$ is real-valued, it follows that ${\mathrm{\Phi }}_{ϵ}({\mathbf{r}}_1,{\mathbf{r}}_2)$ satisfies the pair of wave equations

Substracting Eq. (84) from Eq. (83) yields

We now perform the change of variables

Equation (85) thus becomes

Finally, we multiply Eq. (88) by $e^{i\mathbf{k}·{\mathbf{r}}^{\prime }}/{(2\pi )}^3$, integrate with respect to ${\mathbf{r}}^{\prime }$, and make use of the definition Eq. (80) of the Wigner transform. We find, as shown in Appendix B, that $W_{ϵ}$ obeys the Liouville equation

We 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

The Liouville equation (89) thus becomes

Inserting Eq. (91) into Eq. (93) and collecting terms of the same order in $ϵ$, we find that at $O(\sqrt{ϵ})$