## Abstract

We consider the scattering and absorption of light in discrete random media of densely packed spherical particles. In what we term “radiative transfer with reciprocal transactions” (${\text{R}}^{2}{\text{T}}^{2}$), we introduce a volume element of the random medium, derive its scattering and absorption characteristics using the superposition $\text{T}$-Matrix method (STMM), and compute its frequency-domain incoherent volume-element scattering characteristics. Using an order-of-scattering approach, we then compute a numerical Monte Carlo solution for the scattering problem with an exact treatment of the interaction between two volume elements. We compute both the direct and reciprocal contributions along a sequence of volume elements, allowing us to evaluate the coherent backscattering effects. We show that the ${\text{R}}^{2}{\text{T}}^{2}$ and exact STMM solutions are in mutual agreement for large finite systems of densely packed spherical particles. We conclude that the ${\text{R}}^{2}{\text{T}}^{2}$ method provides a viable numerical solution for scattering by asymptotically infinite systems of particles.

© 2018 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

Electromagnetic scattering and absorption in particle systems much larger than the wavelength constitute an open computational problem in classical electromagnetics. In what follows, a novel approximate multiple scattering method is introduced for discrete random media of particles, that is, for media composed of separate particles distributed randomly within a given volume. The multiple scattering method utilizes rigorous electromagnetic interactions among the particles and is applicable to large particle systems. For the validation of the method, exact electromagnetic methods are utilized.

Earlier studies on scattering by dense discrete random media have been mostly based on the Percus–Yevick approximation (e.g., [1,2]). We start from the radiative transfer coherent backscattering method ($\text{RT}\text{}$-CB) [3,4] developed for sparse media of particles. Here, for dense media, we incorporate full incoherent extinction, scattering, and absorption properties of a volume element larger than the wavelength, extending the earlier studies (see [5,6]). We term the present approach “radiative transfer with reciprocal transactions” (${\text{R}}^{2}{\text{T}}^{2}$). The Superposition $T$-matrix method ([STMM], [7,8]) gives us the means to obtain asymptotically exact solutions for scattering by small finite media of spherical particles, allowing for the validation of the ${\text{R}}^{2}{\text{T}}^{2}$.

This Letter is motivated by two ubiquitous astronomical phenomena observed at small solar phase angles (the Sun-object-observer angle) for the Moon, asteroids, Saturn’s rings, transneptunian objects, and atmosphereless solar system objects at large. First, a nonlinear increase of brightness, commonly called *the opposition effect* (e.g., [9]), is observed toward the zero phase angle in the magnitude scale. Secondly, the scattered light is observed to be partially linearly polarized parallel to the sun-object-bserver plane that is commonly called *negative polarization* ([10]).

In this Letter, we describe the theoretical framework of the ${\text{R}}^{2}{\text{T}}^{2}$. Then we describe the numerical method. After that, we compare the results from the ${\text{R}}^{2}{\text{T}}^{2}$ and STMM methods, as well as those from a formal application of the sparse-medium $\text{RT}$-CB method, for systems of spherical particles with varying system sizes and particle volume fractions. The last section provides conclusions and future prospects.

Consider electromagnetic scattering by a finite, spherical discrete random medium (radius $R$) of $N$ identical spherical particles (radius $a$, complex refractive index $m$; Fig. 1). The Foldy–Lax equations (FLEs; e.g., [11]) describe the electromagnetic scattering by the discrete medium in the frequency domain and are mathematically equivalent to the Maxwell equations. The Neumann series solution to the FLEs can be written as

We express the ensemble-averaged scattered field, in other words, the mean or coherent scattered field, as

The incoherent scattered field of realization $i$ is then obtained by subtracting the coherent scattered field from the scattered field of the realization:

Within the framework of the Neumann series solution of the FLEs in Eq. (1), the incoherent intensity $\u27e8{|{\mathbf{E}}^{\mathrm{s},\mathrm{ic}}|}^{2}\u27e9$ can be computed as an order-of-scattering expansion among any chosen partition of the random medium into particle groups, with each particle belonging to one and only one of the groups. The exact scattering and absorption characteristics are then utilized in the FLEs for each particle group.

Consider next the scattering problem for a discrete random medium approaching infinity in size. The exact treatment with FLEs is inapplicable due to unsurmountable computational challenges. In the approximate multiple scattering method, it is our goal to compute the ensemble-averaged incoherent extinction, scattering, and absorption characteristics for the entire random medium. We incorporate rigorous electromagnetic interaction among the particle groups and introduce what we term incoherent interaction: the incoherent field of particle group A acts as the exciting field for particle group B and, after ensemble-averaging and generating the realization for particle group B, we obtain the incoherent scattered field of particle group B.

In more detail, first, we define a particle group with a size of the order of a few wavelengths, including large numbers of particles. Secondly, we generate sample groups of particles and compute their free-space scattering and absorption characteristics. Thirdly, we compute the ensemble-averaged coherent scattered field of the particle group. The incoherent scattered field of a single particle-group realization follows by subtracting the coherent field from the free-space field according to Eq. (3). The ensemble-averaged incoherent extinction characteristics follow from the ensemble-averaged incoherent scattered and absorbed flux of the particle group. We must monitor that the particle group remains substantially smaller than the resulting mean free path length for the incoherent extinction. Should this condition be violated, the aforedescribed procedure must be repeated with a smaller particle group size, nevertheless fulfilling the condition that there must be large numbers of particles in the group. Fourthly, in the numerical Monte Carlo solution for multiple scattering by discrete random media, we compute the interaction between two particle groups with the help of the incoherent scattered fields. Finally, the multiple scattering solution needs to be independent of changes in particle group size.

Consider $\u27e8{|{\mathbf{E}}^{\mathrm{s}}|}^{2}\u27e9$ available from the absolute-squared Neumann series of Eq. (1). In the resulting series, we can identify the so-called ladder and cyclical interaction diagrams, and they will appear in the sum with a certain weight corresponding to the discrete random medium under consideration. For the current problem of large media with large numbers of particle groups, these ladder and cyclical interaction diagrams will constitute the leading terms (e.g., [2]). Accounting only for the leading terms dramatically simplifies the solution of the FLEs. As Mishchenko *et al.* [11] have shown to be valid for sparse discrete random media, we here hypothesize that, also in the case of dense discrete random media, the effects of the remaining interaction diagrams are accounted for in the incoherent extinction. This implies introducing new weights for the ladder and cyclical interaction diagrams in the original FLEs: these weights are here derived from incoherent extinction that is assumed exponential in its form. In summary, we re-formulate the original FLEs using the particle groups, introducing incoherent extinction within the medium.

The ${\text{R}}^{2}{\text{T}}^{2}$ is a numerical method following the theoretical framework outlined above. In this pilot Letter, we utilize equal-sized spherical volume elements of particles, accepting the caveat of not being able to partition the full space into the given volume elements. We compute the contributions from the ladder and cyclical diagrams by assuming exponential extinction within the random medium. In ${\text{R}}^{2}{\text{T}}^{2}$, we compute the interaction between any two different volume elements rigorously, including all near-field and far-field effects.

The numerical ray tracing proceeds as follows. First, the extinction coefficient ${\kappa}_{\mathrm{e}}^{\mathrm{ic}}$ (in units of inverse distance) and extinction mean free path length ${\ell}_{\mathrm{e}}=1/{\kappa}_{\mathrm{e}}^{\mathrm{ic}}$ are computed from incoherent scattering and absorption by ensemble-averaging over a large number of volume-element realizations. Secondly, incident rays are traced into the random medium based on exponential extinction, allowing for the directly transmitted flux to escape from the medium and for the remaining flux to be traced within the medium. Thirdly, the interaction between the incident field and the volume element is computed exactly by using the STMM method, accompanied by absorption and peel-off of radiation emerging from the medium in discrete directions on a unit sphere. The coherent scattered field is subtracted, and the incoherent scattering characteristics are computed. Fourthly, a scattered ray and the next interaction point are generated from the far-field incoherent scattering characteristics and exponential extinction. Fifthly, the interaction between two volume elements is computed exactly, and, as above, the resulting coherent scattered field is subtracted in order to obtain the incoherent scattering characteristics. Thereafter, absorption and peel-off take place analogously to the procedure for the first scattering. The rays are traced until their flux decreases below a pre-defined cutoff.

For accurate ensemble-averaging over the volume elements of particles, we generate large, periodically continued random media of spherical particles. The fundamental cubic cell is composed of millions of particles, resulting in a negligible variance for the number of particles in the cubic cell. We draw a sample spherical volume element of particles from the center of the cubic cell, repeating the generation of the cell for each sample volume element. Finally, the entire ${\text{R}}^{2}{\text{T}}^{2}$ software has been written in Fortran with parallelization for a supercomputing environment.

We illustrate the ${\text{R}}^{2}{\text{T}}^{2}$ method with four examples of scattering by a spherical discrete random medium of spherical particles. In all cases, the size parameter of the particles is $ka=2$, and the real part of the refractive index is $\mathrm{Re}(m)=1.31$, corresponding to water ice. The volume fraction of the particles is either $v=12.5\%$ or $v=25\%$. In the first two microscopic examples, the size parameters of the entire medium and the volume element are $kR=100$ and $k{R}_{0}=10$, respectively, and the particles are nonabsorbing with $\mathrm{Im}(m)=0$. In the two macroscopic examples, we have $kR=10000$, $k{R}_{0}=10$, and $\mathrm{Im}(m)=1.0\times {10}^{-4}$.

Figure 2 shows the incoherent far-field scattering-matrix elements for the volume elements in the microscopic cases for illustrating the characteristics of incoherent scattering, recalling that the ${\text{R}}^{2}{\text{T}}^{2}$ utilizes rigorous interactions. It is evident that removing the coherent scattered field affects the forward-scattering characteristics in the ${S}_{11}$ element describing the angular dependence of scattered intensity for unpolarized incident light. As to the degree of linear polarization for unpolarized incident light $-{S}_{21}/{S}_{11}$, as well as to the other ratios ${S}_{ij}/{S}_{11}$, the characteristics are not changed.

Figures 3 and 4 show the incoherent scattering matrix elements ${S}_{11}$ and $-{S}_{21}/{S}_{11}$ for the entire random media as computed using the ${\text{R}}^{2}{\text{T}}^{2}$ and the STMM methods. Also included are the formal sparse-medium RT-CB results with the pure Mie scattering characteristics as input. First, it is clear that the sparse-medium RT-CB fails to reproduce the exact results in the present case of densely packed particles. Secondly, the ${\text{R}}^{2}{\text{T}}^{2}$ accurately reproduces both ${S}_{11}$ and $-{S}_{21}/{S}_{11}$ for the backscattering regime of $\alpha \le 90\xb0$. For $\alpha >90\xb0$, there is a deviation in $-{S}_{21}/{S}_{11}$ presumably deriving from collective incoherent surface effects not accounted for in the ${\mathrm{R}}^{2}{\mathrm{T}}^{2}$ method. The STMM results for ${S}_{11}$ show forward diffraction: its approximate modeling is beyond the scope of this Letter.

Figure 5 shows the results for multiple scattering by a large discrete random medium of spherical particles with $kR=10000$, $ka=2$, and $m=1.31+\mathrm{i}{10}^{-4}$, with particle volume fractions of either $v=0.125$ or $v=0.25$. There is striking similarity in the angular dependences. Near the backscattering direction, the angular widths for the coherent backscattering effects are marginally narrower for $v=0.125$ due to the longer extinction mean free path length. These results are realistic: for example, in the case of sparse semi-infinite plane-parallel discrete random media, the pure radiative transfer results are independent of the particle volume fraction.

The time usage of the ${\text{R}}^{2}{\text{T}}^{2}$ depends on multiple parameters, such as the number of rays, flux cutoff value, the size of the volume element, and the mean free path length. For large random media, the convergence is worse especially in the forward scattering regime. The computation of the incoherent $T$-matrices took 1–6 days with a single CPU for each of the four cases. It took an additional 1–3 days and 240 days, respectively, to obtain the ${\text{R}}^{2}{\text{T}}^{2}$ results for the $kR=100$ (10000 rays) and $kR=10000$ cases (27000 rays). The STMM results for the $kR=100$ cases averaged over 128 realizations and 64 scattering planes took around 250 days for $v=0.125$ and 600 days for $v=0.25$, while the RT-CB required only 10–15 min. Due to the high time usage, a computing cluster is required although, by adjusting the parameters, it is possible to run the ${\text{R}}^{2}{\text{T}}^{2}$ on a personal computer in reasonable times when the studied media are small. In addition, the time consumption can be lowered if precomputed incoherent $T$-matrices are available from previous runs.

We have described a novel numerical Monte Carlo method for the computation of multiple scattering by densely packed discrete random media. Utilizing finite systems of spherical particles, we have validated the ${\text{R}}^{2}{\text{T}}^{2}$ method by a comparison to the corresponding exact results available from the STMM method. The agreement among the ${\text{R}}^{2}{\text{T}}^{2}$ and STMM methods is excellent in all cases studied.

There are several pathways for future work. First, whereas we are presently assuming exponential extinction in the random media, it is possible to remove this assumption by incorporating additional steps in the computation. Secondly, for extinction mean free path lengths that are large, compared to the wavelength, we can devise an efficient version of the ${\text{R}}^{2}{\text{T}}^{2}$ based on far-field interactions. Thirdly, it is feasible to extend the ${\text{R}}^{2}{\text{T}}^{2}$ to random media of nonspherical particles. Fourthly, we can study the effect of the surface on the incoherent scattering and absorption characteristics of the volume elements as well as the random media.

## Funding

H2020 European Research Council (ERC) (320773); computational resources provided by CSC—IT Centre for Science Ltd., Finland.

## REFERENCES

**1. **L. Tsang, J. A. Kong, and R. T. Shin, *Theory of Microwave Remote Sensing* (Wiley, 1985).

**2. **L. Tsang and A. Ishimaru, J. Electromagn. Waves Appl. **1**, 59 (1987). [CrossRef]

**3. **K. Muinonen, Waves Random Media **14**, 365 (2004). [CrossRef]

**4. **K. Muinonen, M. I. Mishchenko, J. M. Dlugach, E. Zubko, A. Penttilä, and G. Videen, Astrophys. J. **760**, 118 (2012). [CrossRef]

**5. **L. M. Zurk, L. Tsang, K. H. Ding, and D. P. Winebrenner, J. Opt. Soc. Am. A **12**, 1772 (1995). [CrossRef]

**6. **C. C. Lu, W. C. Chew, and L. Tsang, Radio Sci. **30**, 25 (1995). [CrossRef]

**7. **J. Markkanen and A. J. Yuffa, J. Quant. Spectrosc. Radiat. Transfer **189**, 181 (2017). [CrossRef]

**8. **D. W. Mackowski and M. I. Mishchenko, J. Quant. Spectrosc. Radiat. Transfer **112**, 2182 (2011). [CrossRef]

**9. **N. P. Barabashev, Astron. Nachr. **217**, 445 (1922). [CrossRef]

**10. **B. Lyot, *Recherches sur la polarisation de la lumière des planètes et de quelques substances terrestres* (Observatoire de Paris, 1929), Vol. 8.

**11. **M. I. Mishchenko, L. D. Travis, and A. A. Lacis, *Multiple Scattering of Light by Particles* (Cambridge University, 2006).