Transfer matrix method for four-flux radiative transfer

TL;DR

This paper introduces a transfer matrix formalism for four-flux radiative transfer models, enabling efficient analysis of multilayer scattering media with accurate predictions validated against Monte Carlo simulations.

Contribution

A novel transfer matrix approach for four-flux radiative transfer models tailored for multilayer scattering media, validated against numerical simulations.

Findings

Excellent agreement with Monte Carlo simulations across various parameters.

Accurate prediction of reflection and transmission in multilayer films.

Applicable to diffuse and collimated incident light.

Abstract

We develop a transfer matrix formalism for four-flux radiative transfer models, which is ideally suited for studying transport through multiple scattering layers. The model, derived for spherical particles within the diffusion approximation, predicts the specular and diffuse reflection and transmission of multilayer composite films for diffuse or collimated incidence. The model shows remarkable agreement with numerical Monte Carlo simulations for a range of absorption and film thicknesses, as well as for an example multilayer slab.

Figures

Equations

$$\frac{d{F}_c^{+}}{dz}=-(k+S_f+S_b)F_c^{+},$$

$$\frac{d{F}_c^{-}}{dz}=(k+S_f+S_b)F_c^{-},$$

$$\frac{d{F}_d^{+}}{dz}=-(K+S)F_d^{+}+S{F}_d^{-}+S_f{F}_c^{+}+S_b{F}_c^{-},$$

$$\frac{d{F}_d^{-}}{dz}=(K+S)F_d^{-}-S{F}_d^{+}-S_f{F}_c^{-}-S_b{F}_c^{+},$$

$$F_c^{+}(z)=C_1{e}^{-\lambda z},F_c^{-}(z)=C_4{e}^{\lambda z}$$

$$F_d^{+}(z)=C_1{A}_1{e}^{-\lambda z}+C_2{e}^{-\alpha z}+C_3{e}^{\alpha z}+C_4{A}_4{e}^{\lambda z}$$

$$F_d^{-}(z)=C_1{B}_1{e}^{\lambda z}+C_2{A}_2{e}^{-\alpha z}+C_3{A}_3{e}^{\alpha z}+C_4{B}_4{e}^{\lambda z}$$

$$A_1=B_4=\frac{S{S}_b+S_f(K+S+\lambda )}{{\alpha }^2-{\lambda }^2},$$

$$A_2=A_3^{-1}=\frac{K+2S-\alpha }{K+2S+\alpha },$$

$$A_4=B_1=\frac{S{S}_f+S_b(K+S-\lambda )}{{\alpha }^2-{\lambda }^2},$$

$$\lambda =k+S_f+S_b,\alpha =[k(K+2S){]}^{1/2}.$$

$$S_f={\sigma }_{f}N,S_b=({\sigma }_s-{\sigma }_f)N,$$

$$k=({\sigma }_e-{\sigma }_s)N,$$

$$S=\frac{3}{4}(1-g)(S_f+S_b)-\frac{1}{4}k,K=2k.$$

$$\left[\begin{array}{c}{F}_c^{+}(d)\\ F_c^{-}(d)\\ F_d^{+}(d)\\ F_d^{-}(d)\end{array}\right]=\left[\begin{array}{cccc}{M}_{f,11}& M_{f,12}& M_{f,13}& M_{f,14}\\ M_{f,21}& M_{f,22}& M_{f,23}& M_{f,24}\\ M_{f,31}& M_{f,32}& M_{f,33}& M_{f,34}\\ M_{f,41}& M_{f,42}& M_{f,43}& M_{f,44}\end{array}\right]\left[\begin{array}{c}{F}_c^{+}(0)\\ F_c^{-}(0)\\ F_d^{+}(0)\\ F_d^{-}(0)\end{array}\right],$$

$$M_{f,11}=e^{-\lambda d},M_{f,12}=0,M_{f,13}=0,M_{f,14}=0,$$

$$M_{f,21}=0,M_{f,22}=e^{\lambda d},M_{f,23}=0,M_{f,24}=0,$$

$$M_{f,31}=A_1{e}^{-\lambda d}+\frac{A_1(A_2{e}^{\alpha d}-A_3{e}^{-\alpha d})-2B_1\sinh (\alpha d)}{A_3-A_2},$$

$$M_{f,32}=A_4{e}^{\lambda d}+\frac{A_4(A_2{e}^{\alpha d}-A_3{e}^{-\alpha d})-2B_4\sinh (\alpha d)}{A_3-A_2},$$

$$M_{f,33}=\frac{A_3{e}^{-\alpha d}-A_2{e}^{\alpha d}}{A_3-A_2},M_{f,34}=\frac{2\sinh (\alpha d)}{A_3-A_2},$$

$$M_{f,41}=B_1{e}^{-\lambda d}+\frac{2A_1\sinh (\alpha d)+B_1(A_2{e}^{-\alpha d}-A_3{e}^{\alpha d})}{A_3-A_2},$$

$$M_{f,42}=B_4{e}^{\lambda d}+\frac{A_3{e}^{\alpha d}(A_2{A}_4-B_4)-A_2{e}^{-\alpha d}(A_3{A}_4-B_4)}{A_3-A_2},$$

$$M_{f,43}=-\frac{2\sinh (\alpha d)}{A_3-A_2},M_{f,44}=\frac{A_3{e}^{\alpha d}-A_2{e}^{-\alpha d}}{A_3-A_2}.$$

$$\left[\begin{array}{c}{F}_c^{-}(0)\\ F_c^{+}(d)\\ F_d^{-}(0)\\ F_d^{+}(d)\end{array}\right]=\left[\begin{array}{cccc}{R}_s& T_s& 0& 0\\ T_s& R_s& 0& 0\\ 0& 0& R_d^{+}& T_d^{-}\\ 0& 0& T_d^{+}& R_d^{-}\end{array}\right]\left[\begin{array}{c}{F}_c^{+}(0)\\ F_c^{-}(d)\\ F_d^{+}(0)\\ F_d^{-}(d)\end{array}\right]$$

$$R_d^{+-}=\frac{R_{\varphi }^{+-}+R_J^{+-}}{2-R_{\varphi }^{+-}+R_J^{+-}},T_d^{+-}=1-R_d^{+-},$$

$$R_{\varphi }^{+-}=2{\int }_0^{\pi /2}d\theta \sin (\theta )\cos (\theta )R^{+-}(\theta ),$$

$$R_J^{+-}=3{\int }_0^{\pi /2}d\theta \sin (\theta ){\cos }^2(\theta )R^{+-}(\theta ),$$

$$M_{i,11}=\frac{T_s^2-R_s^2}{T_s},M_{i,12}=\frac{R_s}{T_s},M_{i,13}=M_{i,14}=0,$$

$$M_{i,21}=-\frac{R_s}{T_s},M_{i,22}=\frac{1}{T_s},M_{i,23}=M_{i,24}=0,$$

$$M_{i,31}=M_{i,32}=0,M_{i,33}=\frac{T_d^{+}{T}_d^{-}-R_d^{+}{R}_d^{-}}{T_d^{-}},M_{i,34}=\frac{R_d^{-}}{T_d^{-}},$$

Full text

I Introduction

Modeling light propagation and scattering in multilayer media is relevant to a variety of applications, including atmospheric science Chandrasekhar (2013); Ricchiazzi et al. (1998); Sokoletsky et al. (2014); Meador and Weaver (1980); Emde et al. (2016) and astronomy Whitney et al. (2003); Heng et al. (2014); Leenaarts et al. (2012), remote sensing Platnick (2000); Hedley et al. (2009); Kotchenova et al. (2006, 2008), and biomedical optics Tarvainen et al. (2005); Wang and Wang (2007); Lubbers and Hoffmann (2014); Cheong et al. (1990). For example, clouds and aerosols in the Earth’s atmosphere Meador and Weaver (1980); Emde et al. (2016), and the photosphere and chromosphere in the solar atmosphere Leenaarts et al. (2012), can be modeled as multiple plane-parallel layers to study the impact of radiative transfer on climate and temperature variations. In remote sensing, plane-parallel radiative transfer models are used to correct for atmospheric scattering to improve image clarity Kotchenova et al. (2006, 2008), while in biomedical optics, layered models of blood Lubbers and Hoffmann (2014) and tissue Cheong et al. (1990) can be used to improve imaging techniques such as diffuse optical tomography used in neuroscience and cancer detection Edstrom (2007); Eggebrecht et al. (2014).

Several methods have been developed to solve the radiative transfer equation for the plane-parallel geometry. A common approach is to use Monte Carlo or raytracing simulations Wang and Wang (2007); Iwabuchi (2015), in which the photon trajectory is modeled as a random walk with a step length and direction determined probabilistically by the mean free path and phase function, respectively. However, Monte Carlo methods are computationally expensive and time consuming, particularly for strongly scattering media. On the other hand, analytical models decompose the scattered light into a finite number of fluxes. The simplest model is the two-flux Kubelka-Munk model Kubelka (1948); Vargas and Niklasson (1997a), which decomposes the flux into two isotropic diffuse fluxes propagating in the forward and backward directions perpendicular to the slab. The Kubelka-Munk model has found widespread use in paint design Orel et al. (1997), radiation barriers Slovick et al. (2015), and computer graphics Haase and Meyer (1992). However, because it assumes diffuse incident light, it cannot predict specular reflectance and is only accurate for thick films Vargas and Niklasson (1997a), much thicker than the mean free path.

Alternatively, four-flux models decompose the total flux into collimated and diffuse components traveling in the forward and backward directions Ishimaru (1997); Maheu et al. (1984); Vargas (1998); Simonot et al. (2016). In principle, four-flux models are more accurate for thin films in which there is a substantial unscattered component. Four-flux models have been applied recently to design thermal barrier coatings Wang et al. (2014), x-ray shielding Haber et al. (2012), dental composites Bayou et al. (2013), and nanotechnology for energy applications Laaksonen et al. (2014). However, a major drawback of the current formalism is that it cannot be easily extended to multiple layers Roze et al. (2001). The reason for this can be traced to the way the interface reflections are included. In general, the solution to the four-flux model contains four integration constants, which are determined by imposing boundary conditions on the flux values at the interfaces of the slab Ishimaru (1997). While this process is straightforward for a single layer, each additional layer adds four integration constants, and so the complexity increases rapidly with the number of layers.

In this article, we resolve this issue by developing a transfer matrix formalism for the four-flux model. We derive closed-form analytical expressions for the $4\times 4$ transfer matrices of both the film and interfaces. Multiplication of these matrices results in a total transfer matrix, from which the diffuse and specular reflection and transmission can be calculated. Our model is derived within the diffusion approximation, which assumes the diffuse flux is nearly isotropic Wang and Wang (2007); Gemert et al. (1987), and for spherical particles, though the formalism is applicable to other particle shapes. The model shows remarkable agreement with Monte Carlo simulations for a range of absorption values and film thicknesses, as well as for an example multilayer slab. The closed-form and accuracy of the model suggest it can be an efficient and reliable tool for studying radiative transfer for a variety of applications.

II Four-flux model

In the four-flux model, the total flux in the film is decomposed into four fluxes: two collimated or specular fluxes traveling perpendicular to the slab in the forward and backward directions, $F_{c}^{+}(z)$ and $F_{c}^{-}(z)$, and two diffuse fluxes traveling in the forward and backward directions, $F_{d}^{+}(z)$ and $F_{d}^{-}(z)$ Ishimaru (1997). The collimated fluxes decay due to absorption and scattering as Ishimaru (1997); Maheu et al. (1984); Vargas (1998)

[TABLE]

[TABLE]

where $k$ is the absorption coefficient and $S_{f}$ and $S_{b}$ are the scattering coefficients into the forward and backward hemispheres, respectively. The diffuse fluxes obey the following equations Ishimaru (1997); Maheu et al. (1984); Vargas (1998):

[TABLE]

[TABLE]

where $K$ and $S$ are the absorption and scattering coefficients for the diffuse fluxes. The forward traveling diffuse flux decays due to absorption and scattering, and grows due to scattering of the backward traveling diffuse flux, forward scattering of the forward traveling collimated flux, and backscattering of the backward traveling collimated flux. The same process applies to the backward traveling diffuse flux. The solutions to Eqs. (1)-(4) are Ishimaru (1997); Maheu et al. (1984)

[TABLE]

[TABLE]

[TABLE]

where $C_{1}$, $C_{2}$, $C_{3}$, and $C_{4}$ are integration constants determined by the boundary conditions and

[TABLE]

[TABLE]

[TABLE]

[TABLE]

The forward and backward scattering coefficients for the collimated flux are related to the total scattering cross section $\sigma_{s}$ and the scattering cross section into the forward hemisphere $\sigma_{f}$ by Chylek (1973); Vargas and Niklasson (1997b)

[TABLE]

where $N$ is the particle density. Note that the total scattering coefficient is $S_{f}+S_{b}=\sigma_{s}N$. The absorption coefficient for the collimated flux is given by

[TABLE]

where $\sigma_{e}$ is the extinction coefficient. The scattering and absorption coefficients for the diffuse flux $S$ and $K$ in general depend on the angular distribution of the flux. We assume the Eddington or diffusion approximation Wang and Wang (2007); Gemert et al. (1987), in which the diffuse flux is written as an isotropic term with an anisotropic cosine-dependent perturbation. In the diffusion approximation Gemert et al. (1987); Star et al. (1988)

[TABLE]

III Transfer matrix

First, we derive the transfer matrix for the film. To do so, it is necessary to write the flux values at the back interface ($z=d$) in terms of the flux values at the front interface ($z=0$). This procedure involves evaluating the fluxes in Eqs. (5)-(7) at $z=0$, solving for $C_{1}$, $C_{2}$, $C_{3}$, and $C_{4}$ and substituting these values into the fluxes evaluated at $z=d$. This results in the following transfer matrix equation:

[TABLE]

where the transfer matrix elements are given by

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

The transfer matrix for the entry and exit interfaces is obtained from the known values of the scattering matrix, which relates the outgoing fluxes to the incoming fluxes. Assuming the diffuse fluxes do not reflect or transmit into specular fluxes, and vice versa, the scattering matrix equation for the interfaces is

[TABLE]

where $T_{s}$ and $R_{s}$ are the specular Fresnel transmission and reflection coefficients for normal incidence and $T_{d}^{+-}$ and $R_{d}^{+-}$ are the transmission and reflection coefficients for the diffuse flux, where $+$ and $-$ denote incidence from $z=0$ and $z=d$, respectively. In the diffusion approximation, these coefficients are given by Haskell et al. (1994); Wang and Wang (2007)

[TABLE]

where

[TABLE]

[TABLE]

where $R^{+-}(\theta)$ is the polarization-averaged angular dependent Fresnel reflectivity Haskell et al. (1994); Wang and Wang (2007). Using the $S$ to $M$ matrix conversion procedure described in the appendix, the interface transfer matrix elements are

[TABLE]

[TABLE]

[TABLE]

[TABLE]

The primary advantage of the transfer matrix formalism is the ability to obtain an overall transfer matrix for a multilayer slab by multiplying the transfer matrices of the individual layers and interfaces $\textbf{M}_{j}$ as

[TABLE]

For example, the overall transfer matrix of a single film with two similar interfaces is given by $\textbf{M}_{i}\textbf{M}_{f}\textbf{M}_{i}$. The diffuse and specular reflection and transmission are then calculated using the $M$ to $S$ matrix conversion procedure described in the appendix.

IV Monte Carlo Validation

To validate the four-flux transfer matrix model, we performed Monte Carlo simulations using the commercial raytracing software LightTools by Optical Research Associates. All simulations trace 60,000 rays with a relative ray power threshold of $10^{-5}$ and employ a probabilistic ray split at the interfaces. The first case we considered is intended to evaluate the accuracy of the four-flux model for different amounts of absorption in the particles. Figure 1 shows the calculated specular and diffuse transmission and reflection for a 100-$\mu$m thick slab with a matrix index of 1.5 and 30% loading of 500-nm particles with a real refractive index of 2.5 and different values of the imaginary refractive index.

In the case of zero absorption, shown in Fig. 1(a) and (b), there is remarkable agreement between the models. As the absorption increases, shown in Fig. 1 (c)-(f), the agreement remains excellent but decreases slightly, the most notable difference being the slightly larger diffuse reflection predicted by the four-flux model. The likely reason for this can be traced to the assumed angular distribution of the diffuse flux in the four-flux model. The four-flux model assumes the diffuse flux is nearly isotropic, and thus generally more diffuse than the actual flux. As a result, the average penetration depth is shallower in the four-flux model, leading to less absorption and more diffuse reflection.

Next, we consider the accuracy of the model for different film thicknesses, assuming no absorption. In general, we expect the four-flux model to be more accurate for thicker films, as the assumption of nearly isotropic diffuse flux is valid. Figure 2 shows the calculated specular and diffuse transmission and reflection for the composite in Fig. 1 for slab thicknesses of 10 $\mu$m [(a),(b)] and 1 $\mu$m [(c),(d)]. Overall, we find excellent agreement between the four-flux model and the Monte Carlo simulations, though, as expected, the accuracy of the four-flux model decreases as the film thickness decreases. In particular, in the strong Mie scattering region around 1 $\mu$m the four-flux model predicts higher diffuse reflection and lower diffuse transmission. This can be explained by the directionality of the scattering in the Mie region. When the scattering is directional and the film is thin compared to the mean free path, the scattered flux will be nearly collimated, which contradicts the assumption of nearly isotropic diffuse flux in the four-flux model.

Finally, to illustrate the functionality of the transfer matrix formalism, we considered a two-layer composite consisting of two 50-$\mu$m slabs with 30% loading of particles with an index of 2.5. The matrix indices and particle sizes in the first and second layers are 1.5 and 2 and 500 nm and 1000 nm, respectively. Figure 3 shows a comparison of the four-flux model to the Monte Carlo simulation for the multilayer film. Again we find excellent agreement between the two models, though the four-flux model predicts slightly higher diffusion reflection and lower diffuse transmission for wavelengths between 4 and 8 $\mu$m, which can again be explained by the assumption of a nearly isotropic diffuse flux in the four-flux model.

V Summary

We have developed a closed-form transfer matrix formalism for the four-flux theory of radiative transfer, which is ideally suited for studying multilayered turbid media. The model provides expressions for both the transfer matrix of the inhomogeneous film as well as the interfaces. The model has been derived within the diffusion limit, which approximates the scattered flux as nearly isotropic, and for spherical particles, though the formalism is generally applicable to other shapes provided that the scattering parameters can be calculated. The model shows remarkable agreement with numerical Monte Carlo simulations for a wide range of absorption values and film thicknesses, as well as for an example multilayer slab. The closed form and accuracy of the model suggest it can be an efficient and reliable tool for studying radiative transfer for a variety of applications.

VI acknowledgment

This work was funded by the Materials and Manufacturing Directorate of the Air Force Research Laboratory (prime contract # FA8650-11-D-5400) through UES subcontract #S-925-004-0013.

VII Appendix

VII.1 Conversion between S and M matrix

In this section, we provide the procedure to convert from the $M$ matrix to the $S$ matrix, whose elements correspond to the reflection and transmission coefficients. By definition, the $S$ matrix relates the outgoing fluxes to the incoming fluxes as

[TABLE]

To obtain the $S$ matrix elements in terms of the $M$ matrix elements, Eq. (8) is solved for the outgoing fluxes as a function of the incoming fluxes. The $S$ matrix elements are then obtained by equating the coefficients of the incoming fluxes those in Eq. (9). We demonstrate this procedure by calculating the collimated and diffuse transmission and reflection coefficients for collimated incidence. These expressions are considerably simplified when calculated for light incident from the $z=d$ side of the slab, in which case

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where, for example, $T_{cd}$ corresponds to diffuse transmission for collimated incidence. A similar procedure can be followed to convert the $S$ matrix elements to $M$ matrix elements. In this case, Eq. (9) is solved to obtain expressions for the fluxes at $z=d$ as a function of the $S$ parameters, which can be compared to Eq. (8) to determine the $M$ matrix elements.

VII.2 Mie Scattering parameters

Mie theory provides an exact solution for the scattering parameters of spherical particles in terms of the multipole Mie scattering coefficients Bohren and Huffman (2008)

[TABLE]

[TABLE]

where $n_{p}$ and $n_{m}$, respectively, are the refractive index of the particle and matrix, $k_{m,p}=n_{m,p}2\pi/\lambda$, $\psi_{n}(x)=xj_{n}(x)$ and $\chi_{n}(x)=-xy_{n}(x)$ are the Riccati-Bessel functions where $j_{n}(x)$ and $y_{n}(x)$ are the spherical Bessel functions, and $\xi_{n}(x)=\psi_{n}(x)-i\chi_{n}(x)$. The primes denote differentiation with respect to the argument. The scattering cross section is given by

[TABLE]

and the extinction cross section by

[TABLE]

The anisotropy parameter, which measures the average cosine of the scattering angle, is given by Wang and Wang (2007)

[TABLE]

Finally, the scattering cross section for light into the forward hemisphere is given by Chylek (1973); Vargas and Niklasson (1997b)

[TABLE]

where the prime and double prime, respectively, denote odd and even integers, and

[TABLE]

[TABLE]