{'fraction_non_alphanumeric': 0.034274586173320354, 'fraction_numerical': 0.013164556962025316, 'mean_word_length': 4.665489849955869, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Monte Carlo method for polarized radiative transfer in gradient-index media

Light transfer in gradient-index media generally follows curved ray trajectories, which will cause light beam to converge or diverge during transfer and induce the rotation of polarization ellipse even when the medium is transparent. Furthermore, the combined process of scattering and transfer along curved ray path makes the problem more complex. In this paper, a Monte Carlo method is presented to simulate polarized radiative transfer in gradient-index media that only support planar ray trajectories. The ray equation is solved to the second order to address the effect induced by curved ray trajectories. Three types of test cases are presented to verify the performance of the method, which include transparent medium, Mie scattering medium with assumed gradient index distribution, and Rayleigh scattering with realistic atmosphere refractive index profile. It is demonstrated that the atmospheric refraction has significant effect for long distance polarized light transfer.


Introduction
The refractive index of a material is dependent on its component concentration, temperature, pressure, density, etc. When the dependent properties are distributed inhomogeneously, the refractive index of the media will become a function of space. Examples of participating media with gradient refractive index distribution include the natural media such as earth's (or other planets') atmosphere, the ocean water, the hot air/gas of a flame, and artificial materials, such as gradient index lens, gradient index optical fiber, etc.
In gradient-index media, light beam will generally travel along curved trajectories determined by Fermat's principle [13][14][15]. The gradient index distribution will induce special effects in radiative transfer.
The variation of refractive index along the ray trajectory will cause light beam to converge or diverge and hence influence the transport of radiative intensity. Furthermore, the curved ray trajectory will induce rotation of polarization ellipse as described by the Rytov's field-vector rotation law [14]. Hence polarized radiative transfer in gradient-index media is much more complicated than in uniform-index media. The vector radiative transfer theory developed by Chandrasekhar [16] is the basis for analyzing polarized radiative transfer in media with uniform refractive index distribution. However, it is not appropriate for gradient-index media, since the aforementioned effects induced by the curved ray trajectory are not considered in the formulation. Up till now, there were very few works on polarized radiative transfer in gradient-index media. To our knowledge, Lau and Watson [17] were the first to present a derivation of the radiative transfer equation for polarized light transport in gradient-index media. Recently, Zhao et. al. [18] presented a new derivation and a more generalized form of the equation of polarized radiative transfer in -3 -gradient-index media was given.
Solution methods for polarized radiative transfer in uniform-index media have been studied for a long time, and many successful methods have been developed and used widely, such as Monte Carlo method [19][20][21][22][23], vector discrete ordinates method [24][25][26], adding-doubling method [27], etc. In recent years, scalar radiative transfer in gradient-index media has attracted the interest of many researchers, and many simulation methods were developed, such as curved raytracing method [28][29][30], Monte Carlo method [31,32], finite volume method [33][34][35], finite element method [36][37][38], DRESOR method [39], to name a few. However, there are still very few works on the solution of polarized radiative transfer in gradient-index media, though there were some works on the polarized radiative transfer in multilayer medium, such as, the work of Garcia [40,41]. Most recently, Ben et al. [42] presented a method to simulate polarized radiative transfer in gradient-index media, in which a multilayer slab model was used to approximate the gradient index distribution and then the traditional vector radiative transfer equation is solved by a MC method in each layer. This approximation is inappropriate when a light beam is launched in a direction parallel to the layer. As shown in Fig. 1, the actual ray path in gradient-index media should bend toward the gradient direction indicated by the ray equation [43] as illustrated in Fig. 1(a), whereas the multilayer model predicts a straight path perpendicular to the gradient direction as illustrated in Fig. 1  In this work, a Monte Carlo (MC) method is presented to simulate polarized radiative transfer in gradient-index media. The ray equation is solved to the second order to address the effect induced by curved -4 -ray path. The gradient-index media that support planar ray trajectories are considered for demonstration .   Three types of test cases are presented to verify the performance of the presented MC method and for   comparison with other models, which include a case with analytical solution, a case with Mie scattering and an assumed gradient index distribution, and a case concerning Rayleigh scattering with realistic atmosphere refractive index profile.

Equation of polarized radiative transfer in gradient-index media
A light beam with arbitrary polarization can be described by the Stokes parameters vector ( , , , ) T I Q U V  I [16,44,45]. Additional geometric effect exists when light transfer in gradient-index medium, the curved ray paths will cause the Stokes parameters changing during transport even when the medium is transparent. Generally, the governing equation of polarized radiative transfer in an absorbing, emitting and scattering gradient-index medium (GVRTE) can be written along the ray trajectory as [18]   where s is the coordinate along the ray trajectory, n is the refractive index, e κ is the extinction matrix, a κ is the emission matrix, b I is the Stokes parameters vector of emission in the refractive medium, Z is the scattering phase matrix and R is a matrix accounting for the rotation of the polarization ellipse determined by the Rytov's law [14] defined as (2) in which  stands for the torsion of the trajectory.
If the gradient-index medium only supports ray trajectories being planar curves, then 0   and the GVRTE can be simplified as -5 - There are two groups of media possess this property, such as the case that the refractive index distribution where 0 t is a fixed vector and the case of radial gradient index distribution with spherical symmetry. For isotropic random media, the apparent extinction and absorption process should be polarization independent, then Eq. (3) can be further written as where  , a  and  are the extinction coefficient, absorption coefficient and single scattering albedo, respectively, and ( , )  Z Ω Ω is the normalized scattering phase matrix satisfying 11 4 Note that 11 Z is also the corresponding phase function for scalar radiative transfer.
In this study, only the gradient-index media that support planar ray trajectories are considered. A curved ray tracing MC method is presented to simulate the transport of polarized photons in gradient-index medium described by Eq. (4).

Monte Carlo simulation
The Stokes parameters for quasi-monochromatic beam are given as [46] ** Here we choose the parallel direction vector p , perpendicular direction vector  p and the beam propagation direction vector   Ω pp to form a left-hand reference frame as was used in Refs. [46,47].
The local direction vectors p and  p are defined as the direction being parallel and perpendicular to the meridian plane, which is defined based on the z-axis of the global coordinates system. The definitions of global coordinates system and the local direction vectors are shown in Fig. 2. The definition of the reference frame for polarization is also used for beam of arbitrary polarization.
Geometry for the definition of polarization reference frame and definition of angles in scattering calculation.
In the MC simulation, the state of a photon bundle is traced in the medium. For steady state polarized radiative transfer, a photon bundle carries both the power and the information of polarization. In the present approach, the power carried by a photon bundle is described by a weight variable W , and the polarization state is described by a Stokes parameters vector with Due to moving, scattering and absorption, the properties of a photon bundle, such as location, propagation direction, polarization state and weight will change when it transfers in the gradient-index media.
In the following, the MC simulation of the transport process is presented.

Moving in gradient-index medium
The moving of photon bundle is traced in a global Cartesian coordinate system. In gradient-index media, photon travels in curved ray path determined by the ray equation [43].
 Ω (7) The The free traveling path length l of a photon bundle can be determined in the same way as in traditional MC method for radiative transfer. Considering the extinction coefficient is homogeneous, the free traveling path length l of the photon bundle can be simulated as [20,48] 1 ln l lR   (11) where l R is a random number of uniform distribution in [0,1] . After l is determined, the photon bundle is then moved along the curved ray trajectory using Eq. (8) step by step until tracing distance of each step. In the implementation, the tracing step size is fixed except for the final step to reach the free traveling path, where a tuned step size is used to make the total tracing distance exactly equal to l .

Interaction: absorption, scattering and reflection
Three kinds of interaction are considered here, namely, absorption and scattering in media, and reflection at interfaces. As indicated by the GVRTE, the Monte Carlo models for absorption and scattering of polarized light in uniform-index media should also apply for gradient-index media. There were a few works on Monte Carlo method for polarized radiative transfer, such as Refs. [19][20][21][22]. Though the basic principle is the same, there are a few variants in the detailed implementation of Monte Carlo methods [20]. In the -9 -following, the implementation of present Monte Carlo simulation is presented.
During the interaction with media, the power of a photon bundle (represent as a weight variable W ) will gradually be reduced due to absorption and finally be terminated when it is smaller than a given threshold. When interaction with media occurs, namely, the free traveling distance is reached, some portion of the power is absorbed and some portion is scattered. The absorbed power is calculated as (1 ) If the weight W  is greater than the threshold, the photon bundle is then scattered. It is assumed that the photon bundle is scattered to a specific direction as a whole. The scattered direction and change of polarization state are two key issues in the simulation of scattering process, which are determined from the scattering phase matrix. By definition, the relation between the scattered Stokes parameters I and the incident Stokes parameters vectors ( , , , ) where  Ω and Ω denotes the incident direction and scattering direction respectively. The phase matrix ( , )  Z Ω Ω is usually expressed based on transformation of the scattering matrix P as [16,44] ( , ) ( (13) in which the rotation matrix The definitions of angles for scattering calculation are shown in Fig. 2.
The first step is to determine the scattered direction. It is convenient to use a new local coordinate system, with  p as the local x direction (direction vector ˆx e ),  Ω as the local z direction (direction vector -10 -ˆz e ), and the y direction (direction vector ˆy e ) is determined to make x-y-z a right hand frame. In this local coordinate system, the scattered direction vector can thus be written as ˆˆˆŝ in cos sin sin cos (15) where  and  are the polar and azimuthal angle of Ω defined based on ˆz e and ˆx e in this local coordinate system, and (16) in which os sin cos 1 Note that a similar expressions of Ω in the local coordinates was presented in Ref. [20]. From Eq. (17), the scattered direction can be calculated if a statistical model of  and  is built based on the scattering matrix P . Here the six element scattering matrix for spherical particle is considered, namely, 12 (18) where cos  Ω Ω . Note that     in the local coordinate system, using the relations Eqs. (12), (13) and (18) Because the form of phase function is usually complicated, here rejection method is used to select the scattering directions that satisfy the probability distribution given by Eq. (19). A valid random scattering -11 - In the following, we consider the reflection at the diffuse reflective boundary and an interface between two media. When a photon bundle reaches a boundary or interface, it is considered to be reflected, transmitted or absorbed as a whole. At a Lambertian reflection boundary with reflectance  , the weight is   (28) Note that a similar form of R M was also presented in Ref. [49]. It is noted that Eq. (27) is used because the first component of the polarization state vector should always be unit by definition.

Photon detection and post-processing
For steady state MC simulation, the weight of a photon bundle is considered to have dimension of power (watt). The amount of power received by a detector is calculated by the summation of all the weights of the photon bundles hit on the detector,

Results and discussion
In order to verify the performance of the MC method presented above and to provide examples for comparison with other models, three types of cases are presented. The first is a case for polarized radiative transfer in a non-participating medium with linear gradient index distribution. Analytical results can be obtained for this case due to its simplicity. The second case considers a Mie scattering medium with -14 -collimated beam irradiation. The third case is about sunlight irradiation in a Rayleigh scattering atmosphere with horizontal ray path, which is the path used in limb remote sensing. It will be demonstrated that gradient index distribution will significantly affect both the angular distribution of intensity and the polarization state.
The importance of gradient index distribution in atmosphere on polarized light transport is illustrated. The simulation time was evaluated for the simulation of different cases using a computer with a 2.5GHz CPU.

Case 1: non-participating slab with linear refractive index distribution
To verify the presented method, polarized radiative transfer in a non-scattering and non-absorbing slab with gradient index distribution is considered. The refractive index distribution in the slab is given as ( ) 1.2 0.6 / n z z H  (32) in which H is the height of the slab. Because the medium is non-participating, the Stokes parameters will follow the generalized Clausius invariant for plane ray path [18], namely, There are two different conditions of ray path, a beam transfer along z-direction or the reverse. The beam will converge in the former condition and will diverge in the latter for the given refractive index distribution. Here light irradiation both from the top boundary ( zH  ) and from the bottom boundary ( 0 z  ) are simulated. Figure 3 shows the simulated angular distributions of radiative intensity at different locations in the slab by the MC method, in which the horizontal axis is cos   ,  is the incident polar angle. The analytical results are also shown for comparison. In the MC simulations, 110 7 photon bundles -15 -were used. Since the mean free path will be infinity for non-participating medium, a small optical thickness is taken as an approximation, e.g., 0.001 H   . Intensity detectors with infinite projection area were applied and only the photon bundles fell within a cone with solid angle of  = 0.01 were recorded. The results presented in Fig. 3(a) and 3(b) is for the case with a unit diffuse source ( ( ) 1 irradiated from the top boundary and from the bottom boundary, respectively. From Fig. 3(a), it is seen that the radiative intensity decreases when the location approaches the bottom boundary, which is due to the decrease of refractive index as indicated by Eq. (33). It is also observed that there are some angular range in 0   where the radiative intensity is non-zero for / 0.5 zH  and 1.0. Since the diffuse source irradiates at 0   , this means some of the photons transport to the backward direction, which is attributed to total internal reflection. From Fig. 3(b), it is seen that the angular distribution range of radiative intensity converges with increasing z . In this case, the diffuse source is irradiated at 0 z  , the beam will converge with the increase of refractive index and total internal reflection will not occur. Meanwhile, the magnitude of -16 -radiative intensity increases. As seen from the figure, the results obtained by the MC method agree well with the analytical results within reasonable statistical errors. The location of discontinuities in the intensity distribution are accurately captured, and the beam converging and diverging effects caused by the gradient refractive index distribution are well predicted.

Case 2: Mie scattering of collimated beam in a gradient-index medium
In this case, the combined effects of absorbing, scattering and gradient index distribution are considered.
The case is similar to the 'L=13 problem' of Garcia and Siewert [50], the difference is that the medium is  44. The Legendre coefficients for the phase matrix were used in the scattering calculation, which were obtained by Evens and Stephen [27]. The lower boundary of the -17 -atmosphere is diffuse reflection with  = 0.1. The emission of the atmosphere and the boundary are omitted.
In order for validation, the original 'L=13 problem', namely, with a uniform refractive index distribution, is simulated first. The MC method is applied to obtain the angular distribution of Stokes parameters at plane  =0. Figure 5 shows the angular distribution of the four Stokes parameters obtained using the MC method at / zH =0.5. The results obtained by Garcia and Siewert [50] using the F N method and by using our previous developed spectral element method (SEM) [51] are also shown for comparison. In each of the MC simulation, 8 3 10  photon bundles were launched. The simulation time is about one hour. All the four components of Stokes parameters obtained by the MC results agree well with the reference results, the relative error is less than 1% at most data points. At the direction close to z-axis ( =1   ), the enhancement is the smallest. This can be understood that the ray geometry is more close to straight line at these directions, hence the effect of gradient index distribution is the smallest. At / zH = 0.5 Fig. 6(b), the location of intensity peak at downward direction moves to bigger value of  , indicating the peak is shifted close to the z-axis direction. Furthermore, the width of the peaks become narrower. Both of these are attributed to the beam converging effect caused by gradient index distribution when the beam transfer along the direction where the refractive index gradually increases.

Case 3: Rayleigh scattering atmosphere with gradient index distribution
In the previous two cases, an assumed refractive index distribution was used for demonstration. In this case, we consider a realistic refractive index distribution in the atmosphere. The 1976 version of U.S. standard atmosphere model is used here, which is also a widely used atmosphere model in literatures. The The profiles of temperature and pressure are two key parameters to estimate the refractive index -22 -distribution, which are shown in Fig. 7(a) for the standard atmosphere. The refractivity of dry air is calculated based on Edlen's semi-empirical law [6,9,52]  Nn  is the refractivity,  is wave number (cm -1 ), P is the pressure of dry air (kPa) and T is the air temperature (K). This model is considered to be valid for dry air at altitudes less than 100 km where the mixing ratio between oxygen and nitrogen is fixed and for wavelength from 0.2 to 2,000 m. The refractivity profile for the standard atmosphere calculated based on Eq. (34) for a wavelength of 0.6328 m is plotted in Fig. 7(b). A fitted formula for the refractivity is obtained (also shown in Fig. 7(b)) to ease the curved ray tracing in the MC simulation, which is in the form of exponential function and given as follows. The geometry of the problem is shown in Fig. 8, a collimated light beam with unit radiative flux density is irradiated horizontally to atmosphere from left to right (in the x-direction). The height and width of the simulated region is H = 30 m and L = 40 km, respectively, and the dimension in perpendicular to the paper (y-direction) is assumed infinite. The curvature of the earth surface is neglected. Rayleigh scattering is considered with a single scattering albedo  = 0.99. The ground is assumed to be water with a refractive index 1.33. The radiative flux density of Stokes parameters along the boundaries are recorded. All the photon bundles moving out the simulation domain are terminated. A raytracing test was performed first in order for choosing an appropriate tracing step size. Figure 8 shows the ray paths traced using step size s  = 510 3 L and 110 4 L. The test demonstrates that s  = 510 3 L is sufficient to give accurate ray paths and used in the following simulations. Note that the rays bend toward the ground due to greater refractive index near the ground and reflected at the water surface.
As for verification, we consider a case that Now we consider the case with width L = 40km. When a ray transfers from the source to the right boundary, it will travel a much long distance. In this case, the gradient index distribution of atmosphere will have significant effect on radiative transfer as indicated by the ray path shown in Fig. 8. Figure 10   there exists a height range, i.e. z < 10m, where the radiative flux density has a high value (as shown in Fig.   10(a)), which is even higher than 3 times of the radiative flux density of the source, e.g. for L  = 0.1. This is attributed to the ray bending toward the ground as shown in Fig. 8. Furthermore, due to the reflection at water surface, a sharp enhancement of radiative flux density is induced close to z = 7m. At z > 10m, nearly no direct photons can reach the detectors (shown in Fig. 8). Hence the detected photon flux is from the scattered photons. As shown in Fig. 10(a), the dimensionless radiative flux at this height range increases with optical thickness, this is due to the enhanced atmosphere scattering at big value of optical thickness. As for the uniform index case shown in Fig. 10 where its magnitude has a big value. The magnitude increases with decrease of optical thickness. The Q component is an indication of linear polarization. It has a negative value, which indicates s-polarization. This is an indication of the polarization induced by water surface reflection. At z > 7.5m, Q component radiative flux density is positive and has a relative weak magnitude (on the order 10 -5 to 10 -4 ) as compared to that of z < 7.5m, the latter is negative and has magnitude (on the order 10 -4 to 10 -3 ) about an order greater than the former. Furthermore, as shown in the inset, the magnitude of Q component radiative flux density increases with the increase of optical thickness. This is an indication of scattering enhancement. As for the uniform index case shown in Fig. 10(d), the Q component radiative flux density is positive at different height location for different optical thickness, and the magnitude is on the order 10 -5 to 10 -4 . With the increase of optical -27 -thickness, the magnitude of Q component radiative flux density increases, indicating the polarization induced by scattering process. By comparison, it is shown that the gradient index distribution in atmosphere will have significant effect on polarized light transport, especially for long distance transfer. It can enhance both the detected radiative flux density and the degree of polarization under certain configurations, e.g., the case studied here.

Conclusions
A Monte Carlo method is presented to simulate polarized radiative transfer in gradient-index media, which is demonstrated to be effective to solve polarized radiative transfer problems. The gradient index distribution will make light beam converge or diverge during transfer and induce rotation of polarization ellipse. Furthermore, totally internal reflection may happen for some incident direction. These effects significantly influence the angular distribution of Stokes parameters, even totally alter the angular distribution. As such, omitting gradient index distribution should be taken very carefully.
The combined process of scattering and transfer along curved ray path is very complex, which makes it difficult to interpret the results. For atmosphere with realistic refractive index distribution, it is shown that the gradient index distribution plays an important role for long distance polarized light transport. The gradient index distribution can enhance both the detected radiative flux density and the degree of polarization.
Further study on the effect of gradient index distribution on polarized radiative transfer in scattering media, developing more efficient and accurate solution methods, and application of the theory to engineering problems are very appealing.