1. Introduction

The problem of electromagnetic scattering by a macroscopic medium composed of randomly distributed particles is a subject of great importance to many science and engineering disciplines. If the number of particles is not too large and the overall size of the scattering medium is sufficiently small then this problem can be addressed by means of a direct numerical solution of the Maxwell equations [1–4]. This solution allows one to compute any optical observable anywhere in space. However, the computation of electromagnetic scattering by a medium consisting of a very large number of particles, such as a cloud or a particulate surface, still has to be based on a simplified approximate approach. One of such frequently used approaches is the radiative transfer theory (RTT) intended to describe the transport of electromagnetic energy in a medium consisting of sparsely and randomly distributed particles [5–23].

Until quite recently, the traditional way to introduce the radiative transfer equation (RTE) had been purely phenomenological and essentially required the postulation of the RTE as an artificial supplement to basic physical laws controlling the interaction of macroscopic electromagnetic fields with particles. The numerous inconsistencies and the overall inadequacy of the phenomenological approach have been exposed and discussed thoroughly in [24–27]. Recently, this uncomfortable situation has been rectified by deriving the most general form of the RTE (i.e., applicable to arbitrarily shaped and arbitrarily oriented particles and fully accounting for the vector nature of light) directly from the macroscopic Maxwell equations (MMEs) [28,29]. In contrast to the traditional phenomenological introduction of the RTE, the derivation in [28,29] can be called microphysical since it uses no ad hoc concepts or quantities not already contained in the MMEs [26].

Still, that microphysical derivation was based on the computation of the ensemble-averaged coherency dyad of the total electric field at an observation point, which is not an actual optical observable. Therefore, the physical meaning of the specific intensity column vector as well as that of the very RTE had to be inferred indirectly, through an a posteriori analysis of the resulting RTE. While it was argued on the basis of the integral form of the RTE that the specific intensity column vector had a well-defined physical meaning consistent with that in [5], that argument suffered from the coherency dyad being a purely mathematical construction without explicit physical content.

Because of the use of the coherency dyad of the electric field only, the derivation of the RTE in [28,29] can still be considered incomplete. Indeed, it is well known from classical electromagnetics that the instantaneous local directional flow of electromagnetic energy is described by the Poynting vector involving both the electric and the magnetic field at the observation point [30–32]. Therefore, it would be fundamentally important to establish the physical relationship between the specific intensity column vector and the time-averaged (or, assuming full ergodicity, ensemble-averaged) Poynting vector.

This overarching problem was formulated by Rudolph Preisendorfer 55 years ago ([33], Chapter XIV “Connections with the mainland”), but has never been solved. Preisendorfer himself defined an ad hoc specific intensity that could be linked to the Poynting vector and expected this quantity to satisfy the RTE. However, it remains questionable whether the specific intensity, as defined by Preisendorfer, can even be measured. Also, Preisendorfer’s derivation suffered from the incorrect assumption that the instantaneous electric and magnetic field vectors at any point inside a discrete random medium are always mutually orthogonal (p. 393). It is, therefore, not surprising that Preisendorfer failed to solve the fundamental problem that he himself so eloquently posed.

Unfortunately, by virtue of being a vector product of the electric and magnetic fields, the Poynting vector does not carry explicit information on the polarization state of the scattered electromagnetic field and the field itself. As a consequence, it cannot be used to describe the scattering of electromagnetic waves by particles (particles scatter the electromagnetic field rather than the Poynting vector) and, thus, to derive the RTE.

In this paper we identify and use a more general quantity such that it has the dimension of electromagnetic energy flux, on the one hand, and carries sufficient information about the electric and magnetic fields in order to describe multiple scattering and calculate the resulting Poynting vector at any observation point, on the other hand. The quantity satisfying these requirements in the framework of the frequency-domain formalism is the dyadic product of the magnetic field and the complex conjugate of the electric field. We will call this quantity the Poynting–Stokes tensor (PST).

The main objective of this paper is to obtain a complete and self-consistent microphysical derivation of the RTE directly from the MMEs. We show that the RTE emerges as a bi-product of the theoretical computation of the time (or, equivalently, the ensemble) average of the PST at the observation point and, ultimately, of the expression of the PST in terms of the angular integral of the specific intensity column vector. Our derivation demonstrates that the specific intensity column vector has no fundamental physical meaning and is just an intermediate mathematical quantity in the derivation of a closed-form analytical expression for the ensemble-averaged PST and Poynting vector at the observation point. In particular, the specific intensity is not interpretable in terms of quantifying the amount of electromagnetic energy propagating at a given point in a given direction. Obviously, this outcome is fundamentally at odds with the standard phenomenological notion of the specific intensity as formulated in [5,34] as well as in virtually all subsequent monographs on radiative transfer (RT). Although this result does not negate the RTE, it provides a profoundly different perspective on its physical meaning in the case of electromagnetic scattering by a random particulate medium.

Despite the fact that the specific intensity column vector has no definite physical meaning, this quantity can still be a useful optical characteristic of a turbid medium directly observable with a well-collimated detector of electromagnetic energy. We will demonstrate this by combining the microphysical approach to RT with the physical representation of a well-collimated radiometer as a filter that passes only quasi-plane wavefronts coming from particles located within the acceptance solid angle of the instrument. This result will also explain why the integration of the reading of a well-collimated radiometer over all viewing directions yields the average Poynting vector (and, thus, the net electromagnetic energy flux) at an observation point located inside or outside a sparse random particulate medium.

Consistent with the above rationale, one needs to analyze three fundamental and interrelated aspects of the RTT. The first one deals with the theoretical evaluation of the radiation budget of the entire turbid medium or any part of it. The second one concerns the identification of specific measurement approaches that could be used to quantify the radiation budget experimentally and thereby supplement and/or verify the RTT prediction. The third one is related to specific physical information on the scattering medium that can be imbedded in the solution of the RTE and can potentially be retrieved by using the RTE in the inversion of certain laboratory, in situ, or remote-sensing observations. The following discussion will deal with these three fundamental aspects of the RTT in succession.

In order to avoid redundancy and save space, we use consistently the terminology and notation introduced in [20,28,35]. We denote vectors using the Times bold and Times bold italic fonts and matrices using the Arial bold font. Unit vectors are denoted by a caret, whereas dyads and tensors are denoted by the symbol ↔. The Times italic font is reserved for scalar quantities.

2. General framework

In contrast to various phenomenological approaches to RT, the microphysical theory of electromagnetic scattering by a discrete random medium rests on well-defined assumptions intended to formulate the overall problem in strict physical terms. These assumptions are as follows:

In addition, we will assume that any significant changes in the scattering object (e.g., changes in particle positions and/or orientations with respect to the laboratory reference frame) occur

These two basic assumptions imply that over time intervals long compared to $2\pi /\omega$ but short compared to typical periods of fluctuations of the amplitudes E(r) and H(r), all fields and sources of fields can be considered to be perfectly time-harmonic. As a consequence, the electromagnetic field at any moment in time everywhere in space can be found by solving the frequency-domain differential MMEs [30–32] subject to certain boundary conditions. The specific dependence of the optical constants on spatial coordinates and the corresponding boundary conditions at any moment t are fully defined by the instantaneous geometrical configuration of the N particles (Fig. 1).

Specifically, the frequency-domain monochromatic Maxwell curl equations describing the scattering problem in terms of the time-independent electric and magnetic field amplitudes E(r) and H(r) can be written as follows:

In these equations, $V_{\text{INT}}$ is the cumulative “interior” volume occupied by the particulate scattering object (Fig. 1); $V_{\text{EXT}}$ is the infinite exterior region such that $V_{\text{INT}}\cup V_{\text{EXT}}$ $={\Re }^3,$ where ${\Re }^3$ denotes the entire three-dimensional space; the host medium and the scattering object are assumed to be nonmagnetic; ${\mu }_0$ is the permeability of a vacuum; ${\epsilon }_1$ is the real-valued electric permittivity of the host medium; and ${\epsilon }_2(r,\omega )$ is the complex permittivity of the object.

It can be proven that given the standard boundary conditions for the electric and magnetic fields defined by the specific spatial distribution of the refractive index as well as the so-called radiation condition at infinity, Eqs. (1) have a solution, this solution being unique [37,38]. This fundamental factor makes the MMEs a self-sufficient basis of the electromagnetic scattering theory in general and of the microphysical RTT in particular.

It is important to recognize that there are two sources of randomness of the radiation field in a turbid scattering medium. The first one is the potential quasi-monochromaticity of the incident radiation, as exemplified by sunlight. The second one is the randomness of the particle configuration caused by random changes in particle positions, morphologies, orientations, sizes, and/or refractive indices. However, conventional radiation-budget and remote-sensing applications deal with what can be called the “static” component of the radiation field and are based on the averaging of relevant scattering and absorption characteristics of a turbid medium over time intervals T much longer than T ₂. Therefore, the hierarchy T >> T ₂ >> T ₁ >> $2\pi /\omega$ allows one to split the theoretical computation of electromagnetic scattering by a turbid medium into the following three consecutive steps:

The linearity of this expression makes it applicable to quasi-monochromatic as well as perfectly monochromatic light over time intervals much shorter than T ₂ but much longer than T ₁. Specifically, the linear transformation operator remains constant, whereas the observable characteristics of the incident and scattered radiation are replaced by their averages over a time interval $\text{Δ}t$ such that T ₁ << $\text{Δ}t$ << T ₂:

where the subscript R denotes averaging over all particle coordinates and the subscript ξ denotes averaging over all particle states (i.e., morphologies, orientations, sizes, and refractive indices). Quite often this procedure results in a closed-form equation for ${〈\overline{O_{\text{\hspace{0.05em} sca}}}〉}_{R,\xi }$ or ${〈\widehat{T}〉}_{R,\xi }$ that is much easier to solve than the original MMEs, the RTE being a prime example.

3. Radiation budget of a macroscopic volume element of turbid medium

Let us first consider the radiation-budget problem. In order to characterize the local directional flow of electromagnetic energy resulting from scattering by a complex particulate object such as a cloud, one must calculate the Poynting vector of the total electromagnetic field at the observation point r. The instantaneous value of the Poynting vector is given by the vector product of the real-valued electric and magnetic fields: $S(r,t)=E(r,t)\times H(r,t).$ In the framework of the frequency-domain scattering formalism, the quasi-instantaneous (i.e., averaged over a time interval $2\pi /\omega$<<$\text{Δ}t$<< T ₁) value of the real-valued Poynting vector is equal to the real part of the complex Poynting vector given by

 

Fig. 2 (a) Quasi-instantaneous radiation budget of a volume element $\text{Δ}V$ bounded by a closed surface $\text{Δ}S.$ The arrows represent the distribution of Re[S(r)] over the closed boundary $\text{Δ}S$ corresponding to the specific multi-particle configuration. (b) Configuration-averaged radiation budget of the same volume element evaluated for a statistically uniform spatial distribution of particle positions inside V.

Download Full Size | PDF

However, both the direction and the magnitude of S(r) change in time owing to temporal changes of the multi-particle configuration, thus resulting in a complex speckle pattern rapidly fluctuating in time (Fig. 2a). To suppress the speckle and thereby isolate a static pattern relevant to radiation-budget applications, one needs to average S(r) over a sufficiently long time interval or, equivalently, over all particle positions and states as described in the preceding section. Ideally this would be done by deriving and then solving a closed-form RTE-type equation for ${〈S(r)〉}_{R,\xi }$ and thereby avoiding the impossible task of solving the MMEs directly for a large and statistically representative set of different multi-particle configurations.

4. The Poynting–Stokes tensor

Unfortunately, the Poynting vector carries no information about the polarization state of the field and cannot be used to fully characterize electromagnetic scattering and derive a self-contained equation such as the RTE. A standard descriptor of polarization is the Stokes column vector [39], but it contains no explicit information about the direction of energy propagation and is defined only for transverse electromagnetic waves, whereas the total electromagnetic field at any observation point inside a turbid medium is never a transverse wave.

Another quantity used to describe electromagnetic scattering is the coherency dyad $E(r)\otimes {[E(r)]}^{\ast },$ where ⊗ denotes the dyadic product of two vectors [20,40,41]. This quantity does preserve polarization information, can be applied to an arbitrary field, provides a convenient characterization of scattering by time-variable objects, and can be used to analyze situations in which an object is illuminated by two or more sources of radiation. However, it is defined in terms of the local electric field only and as such does not provide a definitive characterization of the propagation direction and may cause unphysical results in some cases [20].

It is, therefore, necessary to define an alternative quantity which also has the dimension of electromagnetic energy flux while providing a complete and self-contained description of electromagnetic scattering by a turbid medium in the context of practical optical analysis. It is rather obvious that a quantity combining the attributes of all the previously mentioned descriptors of electromagnetic radiation is the PST defined as

Indeed, by definition, the PST is applicable to an arbitrary time-harmonic electromagnetic field. Furthermore, it can be used to find both the Poynting vector and, whenever applicable, the Stokes parameters. Indeed, one has in laboratory Cartesian coordinates (Fig. 1):

where $\widehat{x},$ $\widehat{y},$ and $\widehat{z}$ are the corresponding unit vectors. Also, for a transverse electromagnetic wave propagating in the direction of a unit vector $\widehat{q},$ $H(r)={({\epsilon }_1/{\mu }_0)}^{1/2}\widehat{q}\times E(r).$ Therefore, one has in local spherical coordinates (Fig. 3 ):

(9)$$S=(P_{\varphi \theta }^{\ast }-P_{\theta \varphi }^{\ast })\widehat{q},$$

(10)$$J=\frac{1}{2}\sqrt{\frac{{\epsilon }_{\text{\hspace{0.05em} 1}}}{{\mu }_0}}\left[\begin{array}{c}{E}_{\theta }{E}_{\theta }^{\ast }\\ E_{\theta }{E}_{\varphi }^{\ast }\\ E_{\varphi }{E}_{\theta }^{\ast }\\ E_{\varphi }{E}_{\varphi }^{\ast }\end{array}\right]=\left[\begin{array}{c}{P}_{\varphi \theta }^{\ast }\\ -P_{\theta \theta }^{\ast }\\ P_{\varphi \varphi }^{\ast }\\ -P_{\theta \varphi }^{\ast }\end{array}\right],$$

(11)$$\mathbf{I}=\frac{1}{2}\sqrt{\frac{{\epsilon }_{\text{\hspace{0.05em} 1}}}{{\mu }_0}}\left[\begin{array}{c}{E}_{\theta }{E}_{\theta }^{\ast }+E_{\varphi }{E}_{\varphi }^{\ast }\\ E_{\theta }{E}_{\theta }^{\ast }-E_{\varphi }{E}_{\varphi }^{\ast }\\ -E_{\theta }{E}_{\varphi }^{\ast }-E_{\varphi }{E}_{\theta }^{\ast }\\ \text{i}(E_{\varphi }{E}_{\theta }^{\ast }-E_{\theta }{E}_{\varphi }^{\ast }\text{})\end{array}\right]=\left[\begin{array}{c}{P}_{\varphi \theta }^{\ast }-P_{\theta \varphi }^{\ast }\\ P_{\varphi \theta }^{\ast }+P_{\theta \varphi }^{\ast }\\ P_{\theta \theta }^{\ast }-P_{\varphi \varphi }^{\ast }\\ \text{i}(P_{\theta \theta }^{\ast }+P_{\varphi \varphi }^{\ast })\end{array}\right],$$

where J is the coherency column vector and I is the Stokes column vector [20,35].

 

Fig. 3 Local spherical coordinates. The local Cartesian frame $\{x^{\prime },y^{\prime },z^{\prime }\}$ has the same spatial orientation as the laboratory Cartesian frame {x, y, z}.

Download Full Size | PDF

Thus, the radiation-budget problem for a turbid medium can be fully solved if one can compute theoretically the configuration-averaged PST according to

(12)$${〈\overleftrightarrow{P}(r)〉}_{R,\xi }=\frac{1}{2}{〈H(r)\otimes [E(r)]{}^{\ast }〉}_{R,\xi }$$

since Eq. (8) still applies. Alternatively, as follows from the first equality of Eq. (1), the configuration-averaged PST can be expressed as

(13)$${〈\overleftrightarrow{P}(r)〉}_{R,\xi }=\frac{1}{2\text{\hspace{0.05em} i}\omega {\mu }_0}{[{\nabla }_{r^{\prime }}\times \overleftrightarrow{C}(r^{\prime },r)]|}_{r^{\prime }=r},$$

where

(14)$$\overleftrightarrow{C}(r^{\prime },r)={〈E(r^{\prime })\otimes [E(r)]{}^{\ast }〉}_{R,\xi }$$

is the configuration-averaged so-called dyadic correlation function.

5. Radiative transfer equation

Let us now assume that

• • the incident field is a plane electromagnetic wave propagating in the direction of the unit vector $\widehat{s}$ (Fig. 1);
• • the number of particles N is very large (i.e., tends to infinity);
• • the position and state of each particle are statistically independent of each other and of those of all the other particles;
• • the spatial distribution of the particles throughout the volume V is random and statistically uniform;
• • the scattering medium is convex (which assures that a wave exiting the medium cannot re-enter it);
• • all scattering paths going through a particle more than once can be neglected (the so-called Twersky approximation [42]);
• • all diagrams with crossing connectors in the diagrammatic expansion of the configuration-averaged dyadic correlation function can be neglected (the so-called ladder approximation).

Then the configuration-averaged dyadic correlation function is given by the diagrammatic formula on page 191 of [20] (see also Fig. 16 of [28]). The explicit expanded form of this formula is as follows:

$$\overleftrightarrow{C}(r^{\prime },r)=E_{\text{c}}(r^{\prime })\otimes {[E_{\text{c}}^{}(r)]}^{\ast }$$

$$+n_0{\displaystyle {\int }_{}\text{\hspace{0.17em} d}{R}_1\text{\hspace{0.05em} d}{\xi }_1}\frac{\overleftrightarrow{\eta }\text{\hspace{0.05em} (}{\widehat{r}}_1^{\prime },r_1^{\prime })}{r_1^{\prime }}\cdot {\overleftrightarrow{A}}_1({\widehat{r}}_1^{\prime },\widehat{s})$$$$\cdot {\overleftrightarrow{C}}_{\text{c}}(R_1)\cdot {[{\overleftrightarrow{A}}_1({\widehat{r}}_1,\widehat{s})]}^{\text{T}\ast }\cdot \frac{{[\overleftrightarrow{\eta }\text{\hspace{0.05em} (}{\widehat{r}}_1,r_1)]}^{\text{T}\ast }}{r_1}$$

$$+n_0^2{\displaystyle {\int }_{}\text{\hspace{0.17em} d}{R}_1\text{\hspace{0.05em} d}{\xi }_1}{\displaystyle {\int }_{}\text{\hspace{0.17em} d}{R}_2\text{\hspace{0.05em} d}{\xi }_2}\frac{\overleftrightarrow{\eta }\text{\hspace{0.05em} (}{\widehat{r}}_1^{\prime },r_1^{\prime })}{r_1^{\prime }}\cdot {\overleftrightarrow{A}}_1({\widehat{r}}_1^{\prime },{\widehat{R}}_{12})\cdot \frac{\overleftrightarrow{\eta }\text{\hspace{0.05em} (}{\widehat{R}}_{12},R_{12})}{R_{12}}\cdot {\overleftrightarrow{A}}_2({\widehat{R}}_{12},\widehat{s})$$

$$\cdot {\overleftrightarrow{C}}_{\text{c}}(R_2)\cdot {[{\overleftrightarrow{A}}_2({\widehat{R}}_{12},\widehat{s})]}^{\text{T}\ast }\cdot \frac{{[\overleftrightarrow{\eta }\text{\hspace{0.05em} (}{\widehat{R}}_{12},R_{12})]}^{\text{T}\ast }}{R_{12}}\cdot {[{\overleftrightarrow{A}}_1({\widehat{r}}_1,{\widehat{R}}_{12})]}^{\text{T}\ast }\cdot \frac{{[\overleftrightarrow{\eta }\text{\hspace{0.05em} (}{\widehat{r}}_1,r_1)]}^{\text{T}\ast }}{r_1}$$

$$+n_0^3{\displaystyle {\int }_{}\text{\hspace{0.17em} d}{R}_1\text{\hspace{0.05em} d}{\xi }_1}{\displaystyle {\int }_{}\text{\hspace{0.17em} d}{R}_2\text{\hspace{0.05em} d}{\xi }_2}{\displaystyle {\int }_{}\text{\hspace{0.17em} d}{R}_3\text{\hspace{0.05em} d}{\xi }_3}\frac{\overleftrightarrow{\eta }\text{\hspace{0.05em} (}{\widehat{r}}_1^{\prime },r_1^{\prime })}{r_1^{\prime }}\cdot {\overleftrightarrow{A}}_1({\widehat{r}}_1^{\prime },{\widehat{R}}_{12})\cdot \frac{\overleftrightarrow{\eta }\text{\hspace{0.05em} (}{\widehat{R}}_{12},R_{12})}{R_{12}}$$

$$\cdot {\overleftrightarrow{A}}_2({\widehat{R}}_{12},{\widehat{R}}_{23})\cdot \frac{\overleftrightarrow{\eta }\text{\hspace{0.05em} (}{\widehat{R}}_{23},R_{23})}{R_{23}}\cdot {\overleftrightarrow{A}}_3({\widehat{R}}_{23},\widehat{s})\cdot {\overleftrightarrow{C}}_{\text{c}}(R_3)\cdot {[{\overleftrightarrow{A}}_3({\widehat{R}}_{23},\widehat{s})]}^{\text{T}\ast }$$

(15)$$\cdot \frac{{[\overleftrightarrow{\eta }\text{\hspace{0.05em} (}{\widehat{R}}_{23},R_{23})]}^{\text{T}\ast }}{R_{23}}\cdot {[{\overleftrightarrow{A}}_2({\widehat{R}}_{12},{\widehat{R}}_{23})]}^{\text{T}\ast }\cdot \frac{{[\overleftrightarrow{\eta }\text{\hspace{0.05em} (}{\widehat{R}}_{12},R_{12})]}^{\text{T}\ast }}{R_{12}}\cdot {[{\overleftrightarrow{A}}_1({\widehat{r}}_1,{\widehat{R}}_{12})]}^{\text{T}\ast }\cdot \frac{{[\overleftrightarrow{\eta }\text{\hspace{0.05em} (}{\widehat{r}}_1,r_1)]}^{\text{T}\ast }}{r_1}$$$$+\cdots ,$$

where the following notation is used (see Fig. 4 ):

 

Fig. 4 Geometry showing the quantities appearing in Eq. (15).

Download Full Size | PDF

• • T denotes a transposed dyadic;
• • $n_0=N/V$ is the number of particles per unit volume;
• • the indices 1, 2, 3,… number individual particles;
• • R _i is the position vector of the ith particle;
• • ${\xi }_i$ denotes collectively the state of the ith particle (i.e., its morphology, orientation, size, and refractive index);
• • the vector $r_1$ connects the origin of particle 1 with the observation point r;
• • the vector ${\mathbf{r}}_1^{\prime }$ connects the origin of particle 1 with the observation point $r^{\prime };$
• • $r_1=\left|r_1\right|,$ $r_1^{\prime }=\left|r_1^{\prime }\right|,$ ${\widehat{r}}_1=r_1/r_1,$ ${\widehat{r}}_1^{\prime }=r_1^{\prime }/r_1^{\prime };$
• • the vector $R_{ij}$ connects the origin of particle j with that of particle i;
• • $R_{ij}=\left|R_{ij}\right|,$ ${\widehat{R}}_{ij}=R_{ij}/R_{ij};$
• • ${\overleftrightarrow{A}}_i({\widehat{q}}^{\prime },\widehat{q})$ is the far-field scattering dyadic of the ith particle computed in the particle reference frame for the incidence and scattering directions $\widehat{q}$ and ${\widehat{q}}^{\prime },$ respectively [20,35];
• • $\overleftrightarrow{\eta }(\widehat{q},l)=\exp [\text{\hspace{0.17em} i}{k}_{1}l\overleftrightarrow{I}+\text{\hspace{0.17em} i}2\pi n_0{k}_1^{-1}l{〈\overleftrightarrow{A}(\widehat{q},\widehat{q})〉}_{\xi }]$ is the coherent transmission dyadic, where $k_1=$ $\omega {({\epsilon }_1{\mu }_0)}^{1/2}$ is the wave number, $\overleftrightarrow{I}$ is the unit dyadic, and ${〈\overleftrightarrow{A}(\widehat{q},\widehat{q})〉}_{\xi }$ is the forward-scattering dyadic averaged over all particle states;
• • $E_{\text{c}}(r)=\overleftrightarrow{\eta }[\widehat{s},s(r)]\cdot E^{\text{inc}}(A)$ is the so-called coherent field, where s(r) is the distance between an observation point r and the corresponding entrance point A measured along the direction of incidence (Fig. 4) and $E^{\text{inc}}(A)$ is the incident field at the entrance point; and
• • ${\overleftrightarrow{C}}_{\text{c}}(r)=E_{\text{c}}(r)\otimes {[E_{\text{c}}^{}(r)]}^{\ast }.$

Let us now take into account the far-field condition $k_1{r}_1^{\prime }$ >> 1 as well as assume that ${\overleftrightarrow{A}}_1({\widehat{r}}_1^{\prime },\widehat{q})$ is a rather slowly varying function of $r_1^{\prime }$ and all elements of the dyadic $2\pi n_0{k}_1^{-1}{〈\overleftrightarrow{A}(\widehat{q},\widehat{q})〉}_{\xi }$ are much smaller (in the absolute sense) than the wave number $k_1.$ Then

(16)$${\nabla }_{r_1^{\prime }}\times \frac{\overleftrightarrow{\eta }\text{\hspace{0.05em} (}{\widehat{r}}_1^{\prime },r_1^{\prime })}{r_1^{\prime }}\cdot {\overleftrightarrow{A}}_1({\widehat{r}}_1^{\prime },\widehat{q})$$ ≈$$\text{i}{k}_1{\widehat{r}}_1^{\prime }\times \frac{\overleftrightarrow{\eta }\text{\hspace{0.05em} (}{\widehat{r}}_1^{\prime },r_1^{\prime })}{r_1^{\prime }}\cdot {\overleftrightarrow{A}}_1({\widehat{r}}_1^{\prime },\widehat{q}),$$

where we have used the formulas $\nabla \times (f\overleftrightarrow{C})=(\nabla f)\times \overleftrightarrow{C}+f\nabla \times \overleftrightarrow{C}$ and

$$\nabla [r^{-1}\exp (\text{i}kr)]=$$$$(\text{i}k-r^{-1})r^{-1}\exp (\text{i}kr)\widehat{r}\underset{kr>\text{}>1}{\to }\text{\hspace{0.17em} i}k{r}^{-1}\exp (\text{i}kr)\widehat{r}.$$

Let us also integrate over all positions of particle 1 using a local coordinate system with origin at the observation point r, integrate over all positions of particle 2 using a local coordinate system with origin at the origin of particle 1, integrate over all positions of particle 3 using a local coordinate system with origin at the origin of particle 2, etc. Using the notation introduced in Fig. 4 and taking into account that $\text{d}p=r_1^2\text{d}p\text{\hspace{0.05em} d}\widehat{p},$ $\text{d}{R}_{21}=R_{12}^2\text{\hspace{0.05em} d}{R}_{21}\text{\hspace{0.05em} d}{\widehat{R}}_{21},$ and so on, we get from Eqs. (13), (15), and (16):

(17)$${〈\overleftrightarrow{P}(r)〉}_{R,\xi }=-\frac{1}{2}\sqrt{{\epsilon }_{\text{\hspace{0.05em} 1}}/{\mu }_0}{\displaystyle {\int }_{4\pi }\text{d}\widehat{p}\widehat{p}\times {\overleftrightarrow{\Sigma }}_{\text{\hspace{0.05em} L}}(r,-\widehat{p}),}$$

where ${\overleftrightarrow{\Sigma }}_{\text{\hspace{0.05em} L}}(r,-\widehat{p})$ is the ladder specific coherency dyadic defined by

$${\overleftrightarrow{\Sigma }}_{\text{\hspace{0.05em} L}}(r,-\widehat{p})=\text{\hspace{0.17em} δ}(\widehat{p}+\widehat{s}){\overleftrightarrow{C}}_{\text{c}}(r)$$

$$+n_0{\displaystyle {\int }_{}\text{\hspace{0.17em} d}p\text{\hspace{0.17em} d}{\xi }_1}\overleftrightarrow{\eta }\text{\hspace{0.05em} (}-\widehat{p},p)\cdot {\overleftrightarrow{A}}_1(-\widehat{p},\widehat{s})\cdot {\overleftrightarrow{C}}_{\text{c}}(r+p)\cdot {[{\overleftrightarrow{A}}_1(-\widehat{p},\widehat{s})]}^{\text{T}\ast }\cdot {[\overleftrightarrow{\eta }\text{\hspace{0.05em} (}-\widehat{p},p)]}^{\text{T}\ast }$$

$$+n_0^2{\displaystyle {\int }_{}\text{\hspace{0.17em} d}p\text{\hspace{0.17em} d}{\xi }_1}{\displaystyle {\int }_{}\text{\hspace{0.17em} d}{R}_{21}\text{\hspace{0.17em} d}{\widehat{R}}_{21}\text{\hspace{0.17em} d}{\xi }_2\overleftrightarrow{\eta }\text{\hspace{0.05em} (}-\widehat{p},p)\cdot {\overleftrightarrow{A}}_1(-\widehat{p},-{\widehat{R}}_{21})}\cdot \overleftrightarrow{\eta }\text{\hspace{0.05em} (}-{\widehat{R}}_{21},R_{21})$$

$$\cdot {\overleftrightarrow{A}}_2(-{\widehat{R}}_{21},\widehat{s})\cdot {\overleftrightarrow{C}}_{\text{c}}(r+p+R_{21})\cdot {[{\overleftrightarrow{A}}_2(-{\widehat{R}}_{21},\widehat{s})]}^{\text{T}\ast }\cdot {[\overleftrightarrow{\eta }\text{\hspace{0.05em} (}-{\widehat{R}}_{21},R_{21})]}^{\text{T}\ast }$$

$$\cdot {[{\overleftrightarrow{A}}_1(-\widehat{p},-{\widehat{R}}_{21})]}^{\text{T}\ast }\cdot {[\overleftrightarrow{\eta }\text{\hspace{0.05em} (}-\widehat{p},p)]}^{\text{T}\ast }$$

$$+n_0^3{\displaystyle {\int }_{}\text{\hspace{0.17em} d}p\text{\hspace{0.17em} d}{\xi }_1}{\displaystyle {\int }_{}\text{\hspace{0.17em} d}{R}_{21}}\text{d}{\widehat{R}}_{21}\text{\hspace{0.17em} d}{\xi }_2{\displaystyle {\int }_{}\text{\hspace{0.17em} d}{R}_{32}\text{\hspace{0.17em} d}{\widehat{R}}_{32}\text{\hspace{0.17em} d}{\xi }_3\overleftrightarrow{\eta }\text{\hspace{0.05em} (}-\widehat{p},p)}\cdot {\overleftrightarrow{A}}_1(-\widehat{p},-{\widehat{R}}_{21})$$

$$\cdot \overleftrightarrow{\eta }\text{\hspace{0.05em} (}-{\widehat{R}}_{21},R_{21})\cdot {\overleftrightarrow{A}}_2(-{\widehat{R}}_{21},-{\widehat{R}}_{32})\cdot \overleftrightarrow{\eta }\text{\hspace{0.05em} (}-{\widehat{R}}_{32},R_{32})\cdot {\overleftrightarrow{A}}_3(-{\widehat{R}}_{32},\widehat{s})$$

$$\cdot {\overleftrightarrow{C}}_{\text{c}}(r+p+R_{21}+R_{32})\cdot {[{\overleftrightarrow{A}}_3(-{\widehat{R}}_{32},\widehat{s})]}^{\text{T}\ast }\cdot {[\overleftrightarrow{\eta }\text{\hspace{0.05em} (}-{\widehat{R}}_{32},R_{32})]}^{\text{T}\ast }$$

(18)$$\cdot {[{\overleftrightarrow{A}}_2(-{\widehat{R}}_{21},-{\widehat{R}}_{32})]}^{\text{T}\ast }\cdot {[\overleftrightarrow{\eta }\text{\hspace{0.05em} (}-{\widehat{R}}_{21},R_{21})]}^{\text{T}\ast }\cdot {[{\overleftrightarrow{A}}_1(-\widehat{p},-{\widehat{R}}_{21})]}^{\text{T}\ast }\cdot {[\overleftrightarrow{\eta }\text{\hspace{0.05em} (}-\widehat{p},p)]}^{\text{T}\ast }$$$$+\cdots .$$

Note that p ranges from zero at the observation point r to the corresponding value at the point where the straight line in the $\widehat{p}$ direction crosses the boundary of the medium (point $C_1$ in Fig. 4), $R_{21}$ ranges from zero at the origin of particle 1 to the corresponding value at point $C_2,$ etc. Importantly, the ladder specific coherency dyadic has the dimension of specific intensity or radiance (Wm^–2sr^–1) rather than that of intensity (Wm^–2).

It can be readily verified that the ladder specific coherency dyadic satisfies the following closed-form integral RTE:

$${\overleftrightarrow{\Sigma }}_{\text{L}}(r,-\widehat{p})=\text{\hspace{0.17em} δ}(\widehat{p}+\widehat{s}){\overleftrightarrow{C}}_{\text{c}}(r)+n_0{\displaystyle {\int }_{}\text{\hspace{0.17em} d}p}\text{d}{\widehat{p}}^{\prime }\text{\hspace{0.17em} d}\xi \overleftrightarrow{\eta }\text{\hspace{0.05em} (}-\widehat{p},p)\cdot \overleftrightarrow{A}(\xi ;-\widehat{p},-{\widehat{p}}^{\prime })$$

(19)$$\cdot {\overleftrightarrow{\Sigma }}_{\text{\hspace{0.05em} L}}(r+p,-{\widehat{p}}^{\prime })\cdot {[\overleftrightarrow{A}(\xi ;-\widehat{p},-{\widehat{p}}^{\prime })]}^{\text{T}\ast }\cdot {[\overleftrightarrow{\eta }\text{\hspace{0.05em} (}-\widehat{p},p)]}^{\text{T}\ast }.$$

Indeed, using $\text{δ}(\widehat{p}+\widehat{s}){\overleftrightarrow{C}}_{\text{c}}(r)$ as an initial approximation for ${\overleftrightarrow{\Sigma }}_{\text{\hspace{0.05em} L}}(r,-\widehat{p}),$ we can substitute it in the integral on the right-hand side of Eq. (19) and obtain an improved approximation. By continuing this iterative process, we arrive at Eq. (18), which is simply the Neumann order-of-scattering expansion of the ladder specific coherency dyadic with the coherent field serving as the effective initial source of multiple scattering.

The interpretation of Eq. (19) is very transparent: the ladder specific coherency dyadic for a direction $-\widehat{p}$ at a point r consists of a coherent part and an incoherent part. The latter is a cumulative contribution of all particles located along the straight line in the $\widehat{p}\text{\hspace{0.05em} - direction}$ and scattering radiation coming from all directions $-{\widehat{p}}^{\prime }$ into the direction $-\widehat{p}.$

It is easily verified that the well-known transversality of the scattering dyadic causes the transversality of the ladder specific coherency dyadic [20]: $\widehat{q}\cdot {\overleftrightarrow{\Sigma }}_{\text{\hspace{0.05em} L}}(r,\widehat{q})$ $={\overleftrightarrow{\Sigma }}_{\text{\hspace{0.05em} L}}(r,\widehat{q})\cdot \widehat{q}$ $=0,$ where 0 is a zero vector. This allows one to introduce the specific coherency column vector consisting of the four non-zero elements of ${\overleftrightarrow{\Sigma }}_{\text{\hspace{0.05em} L}}(r,\widehat{q}),$

(20)$$\tilde{J}(r,\widehat{q})=\frac{1}{2}\sqrt{\frac{{\epsilon }_{\text{\hspace{0.05em} 1}}}{{\mu }_0}}\left[\begin{array}{c}{\tilde{J}}_1(r,\widehat{q})\\ {\tilde{J}}_2(r,\widehat{q})\\ {\tilde{J}}_3(r,\widehat{q})\\ {\tilde{J}}_4(r,\widehat{q})\end{array}\right]=\left[\begin{array}{c}\widehat{\theta }(\widehat{q})\cdot {\overleftrightarrow{\Sigma }}_{\text{\hspace{0.05em} L}}(r,\widehat{q})\cdot \widehat{\theta }(\widehat{q})\\ \widehat{\theta }(\widehat{q})\cdot {\overleftrightarrow{\Sigma }}_{\text{\hspace{0.05em} L}}(r,\widehat{q})\cdot \widehat{\phi }(\widehat{q})\\ \widehat{\phi }(\widehat{q})\cdot {\overleftrightarrow{\Sigma }}_{\text{\hspace{0.05em} L}}(r,\widehat{q})\cdot \widehat{\theta }(\widehat{q})\\ \widehat{\phi }(\widehat{q})\cdot {\overleftrightarrow{\Sigma }}_{\text{\hspace{0.05em} L}}(r,\widehat{q})\cdot \widehat{\phi }(\widehat{q})\end{array}\right],$$

where $\widehat{\theta }(\widehat{q})$ and $\widehat{\phi }(\widehat{q})$ are unit vectors in the local coordinate system corresponding to the propagation direction $\widehat{q}$ (Fig. 3), as well as the specific intensity column vector

(21)$$\tilde{\mathbf{I}}(r,\widehat{q})=\left[\begin{array}{c}\tilde{I}(r,\widehat{q})\\ \tilde{Q}(r,\widehat{q})\\ \tilde{U}(r,\widehat{q})\\ \tilde{V}(r,\widehat{q})\end{array}\right]=\left[\begin{array}{c}{\tilde{J}}_1(r,\widehat{q})+{\tilde{J}}_4(r,\widehat{q})\\ {\tilde{J}}_1(r,\widehat{q})-{\tilde{J}}_4(r,\widehat{q})\\ -{\tilde{J}}_2(r,\widehat{q})-{\tilde{J}}_3(r,\widehat{q})\\ \text{i}[{\tilde{J}}_3(r,\widehat{q})-{\tilde{J}}_2(r,\widehat{q})]\end{array}\right],$$

where $\tilde{I}(r,\widehat{q})$ is traditionally called the specific intensity (cf. Eqs. (10) and (11)). Note that we use tildes to denote quantities having the dimension of radiance. Equation (19) implies that either specific column vector satisfies an RTE, the classical integro-differential form of the RTE for ${\tilde{\mathbf{I}}}_{}(r,\widehat{q})$ being as follows:

(22)$$\widehat{q}\cdot \nabla \tilde{\mathbf{I}}(r,\widehat{q})=-n_0{〈\mathbf{K}(\widehat{q})〉}_{\xi }\tilde{\mathbf{I}}(r,\widehat{q})+n_0{\displaystyle {\int }_{4\pi }}\text{d}{\widehat{q}}^{\prime }{〈\mathbf{Z}(\widehat{q},{\widehat{q}}^{\prime })〉}_{\xi }\tilde{\mathbf{I}}(r,{\widehat{q}}^{\prime }),$$

where ${〈\mathbf{K}(\widehat{q})〉}_{\xi }$ and ${〈\mathbf{Z}(\widehat{q},{\widehat{q}}^{\prime })〉}_{\xi }$ are the $4\times 4$ extinction and phase matrices, respectively, averaged over all particle states [20,35]. Importantly, the uniqueness of solution of the RTE (22), when supplemented by appropriate boundary conditions, and the fact that both ${〈\mathbf{K}(\widehat{q})〉}_{\xi }$ and ${〈\mathbf{Z}(\widehat{q},{\widehat{q}}^{\prime })〉}_{\xi }$ are real-valued implies that all elements of ${\tilde{\mathbf{I}}}_{}(r,\widehat{q})$ are real-valued. This is straightforward to demonstrate using an order-of-scattering expansion of the RTE. Furthermore, ${〈\mathbf{K}(\widehat{q})〉}_{\xi }$ and ${〈\mathbf{Z}(\widehat{q},{\widehat{q}}^{\prime })〉}_{\xi }$ are sums of pure Mueller matrices [43,44], which implies that $\tilde{I}(r,\widehat{q})$ is always positive.

All results of this section remain valid if one replaces the incident plane electromagnetic wave with a quasi-monochromatic parallel beam of light of infinite lateral extent (cf. Section 8.15 of [20]).

6. Solution of the radiation-budget problem

Rewriting Eq. (17) as

(23)$${〈\overleftrightarrow{P}(r)〉}_{R,\xi }=\frac{1}{2}\sqrt{{\epsilon }_{\text{\hspace{0.05em} 1}}/{\mu }_0}{\displaystyle {\int }_{4\pi }\text{d}\widehat{q}\widehat{q}\times {\overleftrightarrow{\Sigma }}_{\text{\hspace{0.05em} L}}(r,\widehat{q})},$$

using local spherical coordinates to compute the integrand, and recalling Eqs. (20) and (21) yields:

(24)$${〈S(r)〉}_{R,\xi }={\displaystyle {\int }_{4\pi }}\text{d}\widehat{q}\widehat{q}\tilde{I}(r,\widehat{q}).$$

To the extent that the specific intensity $\tilde{I}(r,\widehat{q})$ is real-valued and positive, ${〈S(r)〉}_{R,\xi }$ is also real-valued, while $\widehat{q}\tilde{I}(r,\widehat{q})$ represents a Poynting-vector component in the direction of the unit vector $\widehat{q}.$

Equation (24) completes the microphysical solution of the radiation-budget problem posed above. Indeed, the RTE (22) follows directly from the MMEs upon making specific assumptions about the scattering particulate medium. Therefore, Eq. (24) is also a direct corollary of the MMEs and implies that to compute the local Poynting vector averaged over a sufficiently long period of time one can solve the RTE for the direction-dependent specific intensity column vector and then integrate the direction-weighted first element of ${\tilde{\mathbf{I}}}_{}(r,\widehat{q})$ over all directions. This result is valid for quasi-monochromatic incident radiation as well.

7. Physical meaning of specific intensity

The traditional definition of the specific intensity in the phenomenological RTT states that $\tilde{I}(r,\widehat{q})$ gives the amount of electromagnetic energy transported in the direction $\widehat{q}$ per unit area normal to $\widehat{q}$ per unit time per unit solid angle (e.g., Chandrasekhar 1950). This notion of the specific intensity implies that at the observation point r, electromagnetic energy propagates simultaneously in all directions and does it according to the angular distribution function $\tilde{I}(r,\widehat{q}).$

Our microphysical derivation of Eqs. (23) and (24) directly from the MMEs reveals that in the case of radiative transfer in a turbid medium this interpretation of $\tilde{I}(r,\widehat{q})$ is profoundly incorrect. Indeed, the instantaneous local flow of electromagnetic energy is given by a monodirectional real Poynting vector $S(r,t)=E(r,t)\times H(r,t).$ Averaging over a time interval $2\pi /\omega$<<$\text{Δ}t$<<T ₁ in the framework of the frequency-domain formalism also yields a monodirectional vector $\mathrm{Re}[S(r)]=\text{ (}\text{1}/\text{2}\text{)}\mathrm{Re}\{E(r)\times [H(r)]{}^{\ast }\}$ (Fig. 2a). Averaging over a time interval $\text{Δ}t$>>T ₂ (or, equivalently, over all particle positions and states) still yields a monodirectional Poynting vector given by Eq. (24). Thus under no circumstances is the local flow of electromagnetic energy polydirectional.

This fundamental disagreement between the microphysical and phenomenological approaches is not surprising. Indeed, the computation of the local electromagnetic energy flow in the framework of the microphysical RTT is based on the solution of the MMEs. As a consequence, $\tilde{I}(r,\widehat{q})$ emerges as a purely mathematical quantity entering Eq. (24) and is found by solving an auxiliary intermediate equation, viz., the RTE. As such, $\tilde{I}(r,\widehat{q})$ has no independent physical meaning. To state otherwise is equivalent to claming independent physical meaning for expansion coefficients appearing in an ad hoc mathematical expansion of a function in, e.g., Legendre or Chebyshev polynomials. The same is true of the RTE itself: it is just an intermediate mathematical equation that one has to solve in order to complete the computation of ${〈S(r)〉}_{R,\xi }$ in the framework of macroscopic electromagnetics.

The phenomenological approach turns everything upside down: the actual existence of the specific intensity as a fundamental physical quantity is postulated, the RTE is “derived” as an outcome of verbal “energy balance considerations”, and the MMEs are invoked on an ad hoc basis only at the very last stage in order to compute the single-particle scattering and absorption characteristics entering the RTE. In other words, the MMEs are treated as a supplement to the phenomenological theory rather than as primordial physical equations fully controlling all aspects of the interaction of macroscopic electromagnetic fields with particulate media.

8. Collimated detector of electromagnetic energy flux

Although the specific intensity is not a fundamental physical quantity, it proves to be very useful in practice. To demonstrate that, we need to discuss the physical nature of a measurement with a typical well-collimated detector of electromagnetic energy flux, as shown schematically in Fig. 5a . Let us first consider two plane electromagnetic waves propagating in directions $\widehat{q}$ and ${\widehat{q}}^{\prime },$ respectively. The objective lens of the well-collimated radiometer transforms both plane wavefronts into converging spherical wavefronts with their respective focal points located in the plane of the diaphragm. The pink wavefront passes through the pinhole and is eventually relayed onto the sensitive surface of the photodetector, whereas the blue wavefront gets extinguished by the diaphragm. Thus the combination {objective lens, diaphragm} serves to select only wavefronts propagating in directions very close to the optical axis of the instrument and falling within its small “acceptance” solid angle $\text{Δ}\Omega .$ It is fundamentally important that this directional filter operates in the “wave domain” rather than in the “Poynting-vector domain”.

 

Fig. 5 (a) Optical scheme of a well-collimated detector of electromagnetic energy flux. (b) The well-collimated detector does not always respond to the Poynting vector directed along the optical axis of the instrument.

Download Full Size | PDF

This analysis implies that the well-collimated radiometer does not necessarily react to the local Poynting vector at a point on the objective lens even if this vector is directed along the detector axis. Indeed, let us consider two plane waves propagating in directions ${\widehat{q}}_1$ and $\widehat{q}{}_2,$ respectively (Fig. 5b). Let the instantaneous values of their real magnetic vectors be H ₁ and H ₂, as shown by the magenta arrows, while their instantaneous real electric vectors E ₁ and E ₂ are normal to the paper and are directed towards the reader. The cumulative instantaneous field is represented by the vectors E = E ₁ + E ₂ and H = H ₁ + H ₂, the former again being normal to the paper and directed towards the reader. The resulting instantaneous real Poynting vector S = E × H is shown by the green arrow and is directed along the optical axis of the instrument. However, neither plane wavefront will be passed by the {objective lens, diaphragm} combination, and the instantaneous reading of the detector will be zero.

Of course, the instantaneous electric and magnetic vectors of the two plane waves rotate around the respective propagation directions, and the instantaneous Poynting vector can oscillate in both the magnitude and the direction. Nevertheless, the well-collimated radiometer does not respond even when S is directed along its optical axis and does not accumulate the corresponding component of the time-averaged Poynting vector.

The failure of the well-collimated radiometer to react to the instantaneous Poynting vector in Fig. 5b can be traced to the following fundamental fact: although the Poynting vector is sought at points on the surface of the objective lens, the actual photodetector is invariably located very far from those points. The only circumstance under which the optical system of the well-collimated radiometer can relay the Poynting vector from the surface of the objective lens onto the sensitive surface of the photodetector is when the incident plane wavefront propagates along (or almost along) the optical axis of the radiometer (Fig. 5a). This is true of any well-collimated radiometer irrespective of its specific optical scheme.

Let us now assume that there are several plane wavefronts incident on the objective lens of the well-collimated radiometer in directions ${\widehat{q}}_1,{\widehat{q}}_2,$… (Fig. 6 ). The above discussion implies that the radiometer will select only those wavefronts whose propagation directions fall within the small acceptance solid angle $\text{Δ}\Omega$ (i.e., ${\widehat{q}}_3,{\widehat{q}}_4,$ and ${\widehat{q}}_5),$ and the reading of the detector will be given by the integral of the Poynting vector resulting from the superposition of these “qualifying” wavefronts over the surface of the objective lens.

 

Fig. 6 Response of a well-collimated radiometer to multidirectional illumination.

Download Full Size | PDF

The implications of this analysis for the case of a well-collimated radiometer measuring electromagnetic scattering by a random particulate medium (Fig. 7 ) are profound. Indeed, the instantaneous real electric and magnetic vectors at the surface of the objective lens can be represented as superpositions of the respective vectors of the incident field and the partial fields coming from all the individual particles:

(25)$$E(r,t)=E^{\text{inc}}(r,t)+{\displaystyle \sum _{i=1}^N{E}_i(r,t)},H(r,t)=H^{\text{inc}}(r,t)+{\displaystyle \sum _{i=1}^N{H}_i(r,t)},$$

where the index i numbers the particles. Since the observation point is assumed to be in the far-field zone of any particle, each pair $\{E_i(r,t),H_i(r,t)\}$ represents an outgoing spherical wavelet originating at the center of the respective particle. At a large distance from the particle this wavelet can be considered locally flat. Therefore, the radiometer will select only those wavelets that come from the particles located within the “acceptance volume” of the radiometer $\text{Δ}{V}_{\widehat{q}}$ defined by its acceptance solid angle $\text{Δ}{\Omega }_{\widehat{q}}$ (Fig. 7) and will integrate the resulting Poynting vector over the surface of the objective lens $S_{\text{d}}.$ Assuming for simplicity that $\text{Δ}{\Omega }_{\widehat{q}}$ does not subtend the propagation direction of the incident light $\widehat{s},$ we conclude that the instantaneous reading of the well-collimated radiometer is given by $S_{\text{d}}{\displaystyle {\sum }_{l^{\prime }}{\displaystyle {\sum }_{m^{\prime }}{E}_{l^{\prime }}(r,t)\times H_{m^{\prime }}(r,t)}},$ where the primed indices $l^{\prime }$ and $m^{\prime }$ number particles located inside the acceptance volume $\text{Δ}{V}_{\widehat{q}}.$ Note that since $\text{Δ}\Omega$ is very small, each term in this sum is a vector directed essentially along the unit vector $\widehat{q}$ in Fig. 7. Accumulating this reading over a sufficiently long time interval $\text{Δ}t,$ dividing the result by $\text{Δ}t,$ and assuming ergodicity yields the following average signal per unit area of the objective lens:

(26)$$\overline{{\displaystyle {\sum }_{l^{\prime }}{\displaystyle {\sum }_{m^{\prime }}{E}_{l^{\prime }}(r,t)\times H_{m^{\prime }}(r,t)}}}=\frac{1}{2}\mathrm{Re}{〈{\displaystyle {\sum }_{l^{\prime }}{\displaystyle {\sum }_{m^{\prime }}{E}_{l^{\prime }}(r)\times {[H_{m^{\prime }}(r)]}^{\ast }}}〉}_{R,\xi },$$

where the subscripts R and ξ denote averaging over coordinates and states of all the N particles constituting the scattering medium and not just those located inside $\text{Δ}{V}_{\widehat{q}}.$

 

Fig. 7 A well-collimated radiometer placed inside a random particulate medium. The sizes of the detector and the particles relative to that of the medium are exaggerated for demonstration purposes.

Download Full Size | PDF

9. Practical meaning of specific intensity and specific intensity column vector

The right-hand side of Eq. (26) can be expressed in terms of the Poynting–Stokes dyadic

(27)$$\overleftrightarrow{P}(r;\text{\hspace{0.17em} Δ}{V}_{\widehat{q}})=\frac{1}{2}{〈{\displaystyle {\sum }_{l^{\prime }}{\displaystyle {\sum }_{m^{\prime }}{H}_{l^{\prime }}(r)\otimes {[E_{m^{\prime }}(r)]}^{\ast }}}〉}_{R,\xi }.$$

The far-field order-of-scattering expansion of the total field at on observation point r (see Section 8.1 of [20]) implies that each partial field $\{E_i(r),H_i(r)\}$ with $1\le i\le N$ is a superposition of contributions corresponding to all possible sequences of particles ending at particle i. We can compute $\overleftrightarrow{P}(r;\text{\hspace{0.17em} Δ}{V}_{\widehat{q}})$ by making the standard assumptions invoked previously to calculate ${〈\overleftrightarrow{P}(r)〉}_{R,\xi }$ and in addition requiring that the end particle of any scattering sequence be located inside the acceptance volume $\text{Δ}{V}_{\widehat{q}}.$ The outcome of this lengthy yet straightforward computation is as follows:

(28)$$\overleftrightarrow{P}(r;\text{\hspace{0.17em} Δ}{V}_{\widehat{q}})=\text{\hspace{0.17em} Δ}\Omega \sqrt{{\epsilon }_{\text{\hspace{0.05em} 1}}/{\mu }_0}\widehat{q}\times {\overleftrightarrow{\Sigma }}_{\text{\hspace{0.05em} L}}(r,\widehat{q}),$$

where the unit vector $\widehat{q}$ is directed along the optical axis of the well-collimated radiometer (Fig. 7). It can be verified that Eq. (28) remains valid even if $\widehat{q}=\widehat{s},$ which implies that $\text{Δ}{\Omega }_{\widehat{q}}$ subtends the incidence direction. Furthermore, Eq. (28) applies to the case of quasi-monochromatic as well as monochromatic incident radiation.

The importance of Eq. (28) is difficult to overstate. Indeed, comparing it with Eq. (23) shows, first and foremost, that if all the above assumptions about the scattering particulate medium are valid then the well-collimated radiometer such as that shown in Fig. 7 measures the specific intensity $\tilde{I}(r,\widehat{q})$ provided that the reading of the radiometer is averaged over a sufficiently long period of time. Therefore, by sampling all incoming directions $\widehat{q},$ one can determine the local electromagnetic energy flow according to Eq. (24) and thereby solve the radiation budget problem experimentally.

Secondly, the angular dependence of the measured specific intensity can be analyzed by solving the RTE for a representative range of physical models of the scattering medium, which may yield certain information about the medium. Furthermore, since the well-collimated optical instrument selects only locally plane wavefronts propagating in essentially the same direction $\widehat{q},$ one can add special optical elements and convert the radiometer into a photopolarimeter capable of measuring the entire specific intensity column vector ${\tilde{\mathbf{I}}}_{}(r,\widehat{q}).$ This measurement can contain additional implicit information about particle microphysics which can often be retrieved since ${\tilde{\mathbf{I}}}_{}(r,\widehat{q})$ can also be calculated theoretically by solving the RTE (23).

10. Discussion and conclusions

Equations (17) and (20)–(22) establish, for the first time, the fundamental physical relation between the MMEs and the RTT by clarifying unequivocally how the solution of the RTE enters the local energy budget computation. Furthermore, Eq. (28) demonstrates that well-collimated radiometers (or, more generally, photopolarimeters) measure the specific intensity (specific intensity column vector) entering the RTE and thereby provides the ultimate physical justification for the use of such instruments in radiation-budget and particle-characterization applications. These are the main results of our paper. They are based on specific assumptions listed in Sections 2, 3, and 5 (see also the discussion in [26]) and can be readily extended to the case of an external observation point along the lines of [29]. In what follows, we will discuss these results and their implications, especially in relation to old phenomenological RT concepts and more recent developments.

10.1. The Poynting–Stokes tensor versus the coherency dyadic

Following [40], the microphysical approach to RT pursued in [20,28,29] was based on the calculation of the configuration-averaged coherency dyadic $\overleftrightarrow{C}(r)={〈E(r)\otimes [E(r)]{}^{\ast }〉}_{R,\xi }$ rather than the configuration-averaged PST defined by Eq. (12). The outcome of that calculation,

(29)$${〈\overleftrightarrow{C}(r)〉}_{R,\xi }={\displaystyle {\int }_{4\pi }\text{d}\widehat{p}{\overleftrightarrow{\Sigma }}_{\text{\hspace{0.05em} L}}(r,-\widehat{p}),}$$

involved the same specific coherency dyadic ${\overleftrightarrow{\Sigma }}_{\text{\hspace{0.05em} L}}$ satisfying the same RTE (19), but provided no explicit recipe how to compute the configuration-averaged Poynting vector. As a consequence, the discussion of the physical meaning of the specific intensity in [20,28,29] remained qualitative and, in fact, questionable.

Equation (17) eliminates this remaining ambiguity and allows one to derive the definitive link (24) between the solution of the conventional RTE (22) and the configuration-averaged Poynting vector in the case of electromagnetic scattering by an ergodic turbid medium. Unlike Eq. (17), Eq. (29) lacks the vector pre-multiplication of the specific coherency dyadic by the unit vector $-\widehat{p},$ which does not allow for an explicit derivation of Eq. (24). Obviously, this deficiency of Eq. (29) can be traced back to the use of only the local electric field at the observation point and the resulting lack of explicit information on the local direction of electromagnetic energy flow.

10.2. What is the specific intensity?

We have already emphasized in Section 7 that although the local configuration-averaged Poynting vector ${〈S(r)〉}_{R,\xi }$ can be calculated by integrating the product $\widehat{q}\tilde{I}(r,\widehat{q})$ over all directions $\widehat{q}$ according to Eq. (24), the specific intensity $\tilde{I}(r,\widehat{q})$ obtained by solving the RTE cannot be said to represent the amount of energy per unit solid angle propagating in the direction $\widehat{q}.$ Indeed, it does not matter whether one considers a specific moment in time or averages over a long period of time: the local flow of electromagnetic energy at an observation point always remains monodirectional rather than being distributed over all “propagation directions”.

It is important to recognize that actual physical significance can be attributed only to the integral of $\widehat{q}\tilde{I}(r,\widehat{q})$ over all directions rather than to the individual values of $\tilde{I}(r,\widehat{q})$ corresponding to different directions. For example, one can add to $\tilde{I}(r,\widehat{q})$ any function $f(r,\widehat{q})$ such that ${\displaystyle {\int }_{4\pi }\text{d}\widehat{q}\widehat{q}f(r,\widehat{q})}=0$ and obtain another “specific intensity” yielding the same ${〈S(r)〉}_{R,\xi }.$ A simple example would be any symmetric function such that $f(r,-\widehat{q})=f(r,\widehat{q}).$

To illustrate this point, let us discuss the definition of radiance proposed in Section 124 of [33]. This discussion invokes a hypothetical radiometer, called by Preisendorfer the “radiance flux meter” (Fig. 8 ), which now replaces the well-collimated radiometer in Fig. 7. The sensitive surface of this instrument is assumed to react to the local instantaneous Poynting vector $S(r^{\prime },t)$ only if this vector is directed normally or almost normally to the surface. Specifically, if the direction of $S(r^{\prime },t)$ at any point $r^{\prime }$ of the sensitive surface S _d is within the small acceptance solid angle $\text{Δ}{\Omega }_{\widehat{q}}$ of the detector (e.g., vector S ₁ in Fig. 8) then this vector will contribute to the cumulative reading of the instrument. Instantaneous local Poynting vectors with directions outside $\text{Δ}{\Omega }_{\widehat{q}}$ (e.g., vectors S ₂ and S ₃) are ignored by the instrument and do not contribute to its cumulative reading. Note that unlike the well-collimated radiometer discussed previously, the hypothetical radiance flux meter is a directional filter operating in the Poynting-vector domain. Preisendorfer’s radiance is then defined as follows:

$$\tilde{N}(r,\widehat{q})=\frac{1}{S_{\text{d}}}\underset{\text{Δ}{\Omega }_{\widehat{q}}\to 0}{\lim }\frac{1}{\text{Δ}{\Omega }_{\widehat{q}}}\underset{T\to \infty }{\lim }\frac{1}{T}{\displaystyle {\int }_{S_{\text{d}}}^{}\text{d}{S}^{\prime }{\displaystyle {\int }_t^{t+T}\text{d}{t}^{\prime }|S(r^{\prime },t^{\prime })|\chi [\text{\hspace{0.05em} Δ}{\Omega }_{\widehat{q}},\widehat{s}(r^{\prime },t^{\prime })]}}$$

(30)$$=\underset{\text{Δ}{\Omega }_{\widehat{q}}\to 0}{\lim }\frac{1}{\text{Δ}{\Omega }_{\widehat{q}}}\underset{T\to \infty }{\lim }\frac{1}{T}{\displaystyle {\int }_t^{t+T}\text{d}{t}^{\prime }|S(r,t^{\prime })|\chi [\text{\hspace{0.05em} Δ}{\Omega }_{\widehat{q}},\widehat{s}(r,t^{\prime })]},$$

where $\widehat{q}$ is the unit vector specifying the orientation of the optical axis of the instrument and normal to its sensitive surface S _d centered at r (Fig. 8), $\widehat{s}(r,t^{\prime })=S(r,t^{\prime })/|S(r,t^{\prime })|$ is the unit vector in the direction of $S(r^{\prime },t^{\prime }),$ and

 

Fig. 8 A hypothetical radiance flux meter accumulating local Poynting vectors with directions falling within the acceptance solid angle $\text{Δ}{\Omega }_{\widehat{q}}.$

Download Full Size | PDF

(31)$$\chi [\text{\hspace{0.05em} Δ}{\Omega }_{\widehat{q}},\widehat{s}(r^{\prime },t^{\prime })]=\{\begin{array}{c}1\text{\hspace{0.17em} if}\widehat{s}(r^{\prime },t^{\prime })\in \text{Δ}{\Omega }_{\widehat{q}},\\ \text{0 otherwise . }\end{array}$$

Note that the dependences on $r^{\prime }$ over the detector surface S _d and on t vanish upon averaging over a sufficiently long period of time T. It is, of course, assumed that the diameter of the sensitive surface S _d is much smaller than any dimension of the scattering medium.

It follows from the definition of $\tilde{N}(r,\widehat{q})$ and from the assumption of ergodicity that

(32)$$\overline{S(r,t)}={\displaystyle {\int }_{4\pi }\text{\hspace{0.17em} d}\widehat{q}\widehat{q}\tilde{N}(r,\widehat{q})}={〈S(r)〉}_{R,\xi }={\displaystyle {\int }_{4\pi }\text{\hspace{0.17em} d}\widehat{q}\widehat{q}\tilde{I}(r,\widehat{q})}.$$

Yet, $\tilde{N}(r,\widehat{q})\ne \tilde{I}(r,\widehat{q}).$ Indeed, on one hand, the microphysical derivation of the RTE implies that pairs of multi-particle sequences ending, respectively, at particles 1 and 2 in Fig. 7 do not contribute to $\tilde{I}(r,\widehat{q}).$ On the other hand, the above discussion of Fig. 5b implies that they can contribute to $\tilde{N}(r,\widehat{q})$ at certain moments in time. This result suggests, of course, that $\tilde{N}(r,\widehat{q})$ is not a solution of the RTE (22).

At the first glance, this conclusion may appear a bit paradoxical. We should recall, however, that the microphysical derivation of the RTE is based on the complex-vector frequency-domain formalism, which allows one to conveniently factor out the time-harmonic dependence of the electric and magnetic fields. However, this forces one to define, at the very outset, the Poynting vector and related optical observables as averages over a sufficiently large number of time-harmonic oscillations (Section 3). As a consequence, the frequency-domain formalism cannot be used to model the measurement with an instrument such as the hypothetical Peisendorfer’s radiance flux meter shown in Fig. 8. Indeed, time-harmonic oscillations of the electric and magnetic vectors cause a rapidly oscillating Poynting vector at the detector surface. At certain moments in time $S(r,t)$ can be directed along $\widehat{q}$ and thereby contribute to $\tilde{N}(r,\widehat{q}),$ while at other moments it can be directed along $-\widehat{q}$ and contribute to $\tilde{N}(r,-\widehat{q}).$ While these “opposing” contributions are accounted for and accumulated in the computation of the Presisendorfer’s radiance $\tilde{N}(r,\widehat{q}),$ they substantially (but, of course, not completely) cancel each other in the frequency-domain computation of the specific intensity $\tilde{I}(r,\widehat{q})$ owing to the primordial averaging over a large number of time-harmonic oscillations.

The frequency-domain modeling of the reading of a well-collimated radiometer, such as the one shown in Fig. 7, is free of this problem since all contributing wavefronts propagate in essentially the same direction $\widehat{q}$ at all times.

10.3. Physical nature of radiance measurements

The hypothetical Preisendorfer’s radiance flux meter is an instrument universally applicable, by its very definition, to the measurement of the local time-averaged Poynting vector via $\overline{S(r,t)}={\displaystyle {\int }_{4\pi }\text{d}\widehat{q}\widehat{q}\tilde{N}(r,\widehat{q})}.$ If it were feasible, this measurement would not require many of the assumptions one has to make in order to derive the RTE such as ergodicity and low packing density of the scattering medium. Furthermore, this instrument would serve as an ideal bridge between fundamental electromagnetics and the semi-empirical field of radiometry.

Unfortunately, to the best of the author’s knowledge, such an instrument has never been built, and it remains unclear whether it can be built in principle. Obviously, designing an instrument with sensitivity to the local instantaneous Poynting vector (Fig. 8) would require a detailed quantum-mechanical analysis of the localized and directional light-matter interaction.

Designing well-collimated radiometers and photopolarimeters is a much simpler problem since the selection of macroscopic quasi-plane wavefronts with specific propagation directions is much more straightforward than the directional selection of instantaneous local Poynting vectors and involves easy-to-manufacture macroscopic optical elements such as lenses. In this case, however, the measurement of the local time-averaged Poynting vector via $\overline{S(r,t)}$ ≈${\displaystyle {\int }_{4\pi }\text{d}\widehat{q}\widehat{q}\tilde{I}(r,\widehat{q})}$ is based on numerous assumptions listed in Sections 2, 3, and 5. These assumptions can be expected to be valid in the case of a sparse turbid medium, thereby enabling radiation-budget and particle characterization applications. However, most of these assumptions become quite questionable in the case of densely packed scattering media. Although the well-collimated radiometer would still perform a directional selection of macroscopic quasi-plane wavefronts, the relation of this measurement to $\overline{S(r,t)}$ remains uncertain, and the theoretical modeling of this measurement is likely to require a more complicated approach than solving an RTE.

In summary, our detailed analysis has revealed that the specific intensity is not what it is postulated to be in the phenomenological RTT, and, contrary to a widespread belief, a well-collimated radiometer is not a radiance flux meter universally capable of yielding the time-averaged local Poynting vector. It turns out, however, that in the case of electromagnetic scattering by an ergodic sparse discrete random medium,

• • the well-collimated radiometer (or photopolarimeter) measures $\tilde{I}(r,\widehat{q})$ (or $\tilde{\mathbf{I}}(r,\widehat{q}));$
• • $\tilde{\mathbf{I}}(r,\widehat{q})$ can be computed by solving the conventional RTE; and
• • the time-averaged local Poynting vector can be found by evaluating the integral in Eq. (24).

These firmly established facts make the combination of the RTE and a well-collimated radiometer/photopolarimeter useful in a wide range of applications.

Acknowledgments

The author thanks Brian Cairns, Joop Hovenier, Nikolai Khlebtsov, Daniel Mackowski, Pinar Mengüç, and Victor Tishkovets for numerous illuminating discussions. This research was supported by the NASA Radiation Sciences Program managed by Hal Maring and by the NASA Glory Mission project.

References and links

1. M. I. Mishchenko, L. Liu, D. W. Mackowski, B. Cairns, and G. Videen, “Multiple scattering by random particulate media: exact 3D results,” Opt. Express 15(6), 2822–2836 (2007). [CrossRef]   [PubMed]  

2. M. I. Mishchenko, J. M. Dlugach, and L. Liu, “Azimuthal asymmetry of the coherent backscattering cone: theoretical results,” Phys. Rev. A 80(5), 053824 (2009). [CrossRef]  

3. M. I. Mishchenko, J. M. Dlugach, L. Liu, V. K. Rosenbush, N. N. Kiselev, and Yu. G. Shkuratov, “Direct solutions of the Maxwell equations explain opposition phenomena observed for high-albedo solar system objects,” Astrophys. J. 705(2), L118–L122 (2009). [CrossRef]  

4. M. I. Mishchenko, V. K. Rosenbush, N. N. Kiselev, et al., Polarimetric Remote Sensing of Solar System Bodies (Akedemperoidyka, Kyiv, 2010). http://www.giss.nasa.gov/staff/mmishchenko/books.html

5. S. Chandrasekhar, Radiative Transfer (Dover, New York, 1960).

6. J. E. Hansen and L. D. Travis, “Light scattering in planetary atmospheres,” Space Sci. Rev. 16(4), 527–610 (1974). [CrossRef]  

7. V. V. Sobolev, Light Scattering in Planetary Atmospheres (Pergamon, Oxford, 1975).

8. A. Ishimaru, Wave Propagation and Scattering in Random Media (Academic Press, New York, 1978).

9. H. C. van de Hulst, Multiple Light Scattering (Academic Press, New York, 1980).

10. J. Lenoble, ed., Radiative Transfer in Scattering and Absorbing Atmospheres (A. Deepak, Hampton, Va., 1985).

11. R. M. Goody, and Y. L. Yung, Atmospheric Radiation: Theoretical Basis (Oxford U. Press, Oxford, 1989).

12. C. D. Mobley, Light and Water: Radiative Transfer in Natural Waters (Academic Press, San Diego, Ca., 1994).

13. G. L. Stephens, Remote Sensing of the Lower Atmosphere (Oxford U. Press, New York, 1994).

14. E. G. Yanovitskij, Light Scattering in Inhomogeneous Atmospheres (Springer, Berlin, 1997).

15. G. E. Thomas, and K. Stamnes, Radiative Transfer in the Atmosphere and Ocean (Cambridge U. Press, New York, 1999).

16. K. N. Liou, An Introduction to Atmospheric Radiation (Academic Press, San Diego, 2002).

17. M. Modest, Radiative Heat Transfer (Academic Press, San Diego, Ca., 2003).

18. J. W. Hovenier, C. van der Mee, and H. Domke, Transfer of Polarized Light in Planetary Atmospheres – Basic Concepts and Practical Methods (Springer, Berlin, 2004).

19. A. Marshak, and A. B. Davis, eds., 3D Radiative Transfer in Cloudy Atmospheres (Springer, Berlin, 2005).

20. M. I. Mishchenko, L. D. Travis, and A. A. Lacis, Multiple Scattering of Light by Particles: Radiative Transfer and Coherent Backscattering (Cambridge U. Press, Cambridge, UK, 2006).

21. G. W. Petty, A First Course in Atmospheric Radiation (Sundog Publishing, Madison, Wi., 2006).

22. W. Zdunkowski, T. Trautmann, and A. Bott, Radiation in the Atmosphere (Cambridge U. Press, Cambridge, UK, 2007).

23. A. B. Davis and A. Marshak, “Solar radiation transport in the cloudy atmosphere: a 3D perspective on observations and climate impacts,” Rep. Prog. Phys. 73(2), 026801 (2010). [CrossRef]  

24. Yu. N. Barabanenkov, “Multiple scattering of waves by ensembles of particles and the theory of radiation transport,” Sov. Phys. Usp. 18(9), 673–689 (1975). [CrossRef]  

25. L. A. Apresyan, and Yu. A. Kravtsov, Radiation Transfer (Gordon and Breach, Basel, 1996).

26. M. I. Mishchenko, “Multiple scattering, radiative transfer, and weak localization in discrete random media: the unified microphysical approach,” Rev. Geophys. 46(2), RG2003 (2008). [CrossRef]  

27. M. I. Mishchenko, “Gustav Mie and the fundamental concept of electromagnetic scattering by particles: a perspective,” J. Quant. Spectrosc. Radiat. Transf. 110(14-16), 1210–1222 (2009). [CrossRef]  

28. M. I. Mishchenko, “Vector radiative transfer equation for arbitrarily shaped and arbitrarily oriented particles: a microphysical derivation from statistical electromagnetics,” Appl. Opt. 41(33), 7114–7134 (2002). [CrossRef]   [PubMed]  

29. M. I. Mishchenko, “Microphysical approach to polarized radiative transfer: extension to the case of an external observation point,” Appl. Opt. 42(24), 4963–4967 (2003). [CrossRef]   [PubMed]  

30. J. A. Stratton, Electromagnetic Theory (McGraw-Hill, New York, 1941).

31. J. D. Jackson, Classical Electrodynamics (Wiley, New York, 1999).

32. E. J. Rothwell, and M. J. Cloud, Electromagnetics (CRC Press, Boca-Raton, Fl., 2009).

33. R. W. Preisendorfer, Radiative Transfer on Discrete Spaces (Pergamon Press, Oxford, 1965).

34. E. A. Milne, “Thermodynamics of the stars,” Handbuch der Astrophysik 3, 65–255 (1930).

35. M. I. Mishchenko, L. D. Travis, and A. A. Lacis, Scattering, Absorption, and Emission of Light by Small Particles (Cambridge U. Press, Cambridge, UK, 2002). http://www.giss.nasa.gov/staff/mmishchenko/books.html

36. M. I. Mishchenko, “Multiple scattering by particles embedded in an absorbing medium. 2. Radiative transfer equation,” J. Quant. Spectrosc. Radiat. Transf. 109(14), 2386–2390 (2008). [CrossRef]  

37. S. Silver, ed., Microwave Antenna Theory and Design (McGraw-Hill, New York, 1949).

38. C. Müller, Foundations of the Mathematical Theory of Electromagnetic Waves (Springer, Berlin, 1969).

39. G. G. Stokes, “On the composition and resolution of streams of polarized light from different sources,” Trans. Cambridge Philos. Soc. 9, 399–416 (1852).

40. A. G. Borovoy, ““Method of iterations in multiple scattering: the transfer equation,” Izv. Vuzov Fizika, 6, 50–54 (1966).

41. A. G. Borovoi, “Multiple scattering of short waves by uncorrelated and correlated scatterers,” Light Scattering Rev. 1, 181–252 (2006). [CrossRef]  

42. V. Twersky, “On propagation in random media of discrete scatterers,” Proc. Symp. Appl. Math. 16, 84–116 (1964).

43. C. V. M. van der Mee, “An eigenvalue criterion for matrices transforming Stokes parameters,” J. Math. Phys. 34(11), 5072–5088 (1993). [CrossRef]  

44. J. W. Hovenier and C. V. M. van der Mee, “Testing scattering matrices: a compendium of recipes,” J. Quantum. Spectrosc. Radiat. Transf. 55(5), 649–661 (1996). [CrossRef]