Numerical simulations of radiation transfer in partially-ordered stratified media using Monte Carlo methods

This work includes the research method of terahertz radiation scattering process in partially-ordered stratified media. This method was set on algorithm of computational electrodynamics and Monte Carlo simulations. The method allows to determine differential and integral cross sections and scattering phase functions. The algorithm of modeling radiation transfer in stratified media using Monte Carlo approach was created. Rigorous and Henyey- Grinstein phase functions were considered. Analysis of simulated results was developed. The technique for accounting single scattering was designed.


Introduction
Multiple scattering problem is essential for such fields as atmosphere remote sensing, spectroscopy, biochemistry and others. Different aspects of this problem have been considered in literature [1][2][3][4][5].
According to this work we consider one of the essential cases of multiply scattering: scattering in partially-ordered medium with group formation, when particles with dimensions of the order of wavelength are combined densely in groups, and these groups are randomly distributed in the scattering volume.
The radiation transfer theory is one of the widely used approximated methods to describe multiple scattering processes [6]. It operates with the radiation transfer equation (RTE). In order to solve RTE numerically we need to use the Monte-Carlo (MC) approach [7].This method allows the parameters of inhomogeneous scattering media to be varied easily and different experimental geometries to be simulated.
When medium parameters differ within the volume, averaging over the volume results in inaccuracy of MC simulation. In this case we must divide the volume into layers with different parameters. When every layer consists of structured group of particles, it is rather complicated to understand what is the influence of the scattering phase function on radiation transfer results. The Henny-Greenstein phase function [8] is suitable in the event that medium consists of randomly distributed single particles. However, it is essential to take into account the impact of the interference of scattered waves on the group scattering phase function. The primary aim of the work is to compare Monte Carlo simulation results of using rigorous and Henyey-Greenstein phase functions.

Scattering volume
We consider the arbitrary two-layered medium: each layer consists of randomly distributed groups of nine identical particles with circular cross section in Y Z plane (figure 1b) and permittivity ε P = 1.96. Dimensions of particles (table 1) in every layer are of the order of wavelength λ of TEM plane wave polarized radiation, incident in OZ direction (figure 1a). Table 1. Geometric parameters of particle groups.

Single scattering phase function
The single scattering phase function is the angular distribution of light scattered on object when it is illuminated. In order to find single scattering phase function χ(θ) more accurately we use finite-difference time-domain method (FDTD) [9] of Maxwells equations numerical solution. As a result, we obtain the rigorous group single scattering phase function (figure 2a): where σ(θ) is the group scattering cross section. Henyey-Greenstein phase function χ HG (g, θ) (figure 2b) [8] is calculated from the equation: which is characterized by the asymmetry parameter g: where θ is the scattering angle. Further we use these results for MC numerical simulations in order to find out, whether the interference effects in rigorous χ(θ) have influence on the radiation transfer in two-layered medium.

Monte Carlo simulation
The MC approach for numerical calculation of the intensity radiation scattered within a randomly inhomogeneous turbid medium is based on simulation of the photon trajectories N within the medium between the points of the isotropic source and the detector. MC simulation of light scattering was designed in order to solve RTE for scattering media with parameters from the table 2.
The algorithm consists of three steps [10]: (i) Double-local estimation computation [11] (allows to obtain the scattered radiance in the single point). (ii) Calculation of the collision point. (iii) Calculation of scattering parameters.

Single scattering
Double local estimation in Monte Carlo numerical scheme is used to determine higher order scattered radiation. If it needs to calculate the first-order scattered radiance, we use the following equation [12]: Λε exp(−εr) 4πrsin(θ) π 0 χ(η)exp(−εzsin(θ)(tan(η/2) − tan(θ/2)))dη.  Single scattering albedo, Λ 1.00 1.00 Extinction coefficient, ε ext 20. 9 3.35 Mean free path, l 0.0478 0.2985 Layer thickness, z/l 1.00 1.00 Asymmetry parameter, g 0.5654 0.2978 Number of trajectories 5000 5000 In case of two-layered media it is necessary to change this equation using probabilistic approach. 1. When the single collision is in the second layer, probability is to pass without collisions through the first layer: Value of the single scattering expression 4 decreases proportionally to the probability: 2. When the single collision is in the first layer, probability is to pass without collisions through the second layer: ), The value of the single scattering expression 4 decreases proportionally to the probability: As a result, single-scattered radiance is determined as:

Conclusion
In the present work, the analysis technique of the first-order scattered radiation brightness for medium with 2 layers was shown. Comparative analysis of multiple scattered radiation brightness for medium with two layers shows the following: • interference of the scattered waves inside the group of densely packed particles in stratified medium impacts on the overall radiance distribution, even if the layers have different parameters. • application of model phase function with single parameter for radiation transfer calculation in stratified media leads to significant errors, which can be prevented by using the rigorously obtained phase function.