{'fraction_non_alphanumeric': 0.03363553457432541, 'fraction_numerical': 0.02242368971621694, 'mean_word_length': 4.7804651162790694, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A Discrete Spherical Harmonics Method for Radiative Transfer Analysis in Inhomogeneous Polarized Planar Atmosphere

A discrete spherical harmonics method is developed for the radiative transfer problem in inhomogeneous polarized planar atmosphere illuminated at the top by a collimated sunlight while the bottom reflects the radiation. The method expands both the Stokes vector and the phase matrix in a finite series of generalized spherical functions and the resulting vector radiative transfer equation is expressed in a set of polar directions. Hence, the polarized characteristics of the radiance within the atmosphere at any polar direction and azimuthal angle can be determined without linearization and/or interpolations. The spatial dependent of the problem is solved using the spectral Chebyshev method. The emergent and transmitted radiative intensity and the degree of polarization are predicted for both Rayleigh and Mie scattering. The discrete spherical harmonics method predictions using 36 streams are found in good agreement with benchmark literature results. The maximum deviation between the proposed method and literature results and for absolue value of the polar directions greater than 0.1 is less than 0.5% and 0.9% for the Rayleigh and Mie scattering, respectively. These deviations for directions close to zero are about 3% and 10 % for Rayleigh and Mie scattering, respectively.


Introduction
The multiple scattering of any radiation by the earth atmosphere in all directions at all frequency of the incident beam generally gives birth to a light with a different state of polarization to that of the incident beam (Cornet et al., 2010). The solar radiation in traversing the atmosphere is reflected by the atmospheric particles and also at the level of the ground. The reflected radiation and the outgoing long-wave radiation escape to the space. Information about the atmosphere is collected in measuring that emerged solar radiation (De Rooij, 1985). These measurements should be compared to the computation results for checking the assumptions made theoretically (Hovenier, 1969). Knowing the scattering properties of the atmospheric particles the reflected radiance is calculated for a better interpretation of the observations (De Rooij, 1985). The interest of the interpretation of the observations aims to determine the vertical structure and the physical properties of the atmosphere (De Rooij, 1985;Ackerman and Stokes, 2003). The measurement of the emergent and transmitted solar radiation from the atmosphere and its comparison to computation results are very important, for example, to verify approximations done on the size distribution of the atmospheric particles as functions of the radius. The study of the multiple scattering of polarized light is also important in the study of the energy balance of the earth (Wauben et al., 1993), in the calculation of the infrared of planetary atmospheres (Wauben et al., 1993), in the climate study and the weather prediction (Ackerman and Stokes, 2003), in the remote sensing (Ackerman and Stokes, 2003;Marshak and Davis, 2005;Barlakas et al., 2016) and in the atmospheric optics.
The above mentioned applications of the polarization indicate that the investigation of the polarization characteristics of the light propagating through the atmosphere remains an active task for the researchers. The transport of the polarized radiation in the atmosphere is described by the radiative transfer equations (RTE). This equation is a set of four coupled equations of Stokes parameters that can be put in a vector form called vector radiative transfer equation (VRTE). Although this vector equation has been widely studied in literature, accurate and simple techniques are still needed even for the simplest one-dimensional geometry (Spurr, 2008;Budak et al., 2017). Many methods have been presented in the literature for solving the vector equation of the propagation of the polarization light in the atmosphere including the Monte Carlo method (Cornet et al., 2010;Muñoz and Mills, 2015;Barlakas et al., 2016), ), the "Facile" or FN method (Garcia and Siewert, 1989), the spherical harmonics and variant (Benassi et al., 1985;Garcia and Siewert, 1986;Budak and Korkin, 2006) and the discrete ordinates method and variant (Schulz et al., 1999;Rozanov and Kokhanovsky, 2006;Garcia and Siewert, 2011;Korkin, 2013;Collin et al., 2014),. Each of these methods presents the advantages and the drawbacks. The Monte Carlo method, generally adopted for benchmark solutions, requires a large computer memory and a high computational time. The FN method, which is also generally used for benchmark solutions, converges slowly at the domain boundaries leading to a high approximate order for accurate results. The spherical harmonics method can be used for predictions at any polar and zenith angles within the atmosphere. However, this method requires very high order of angular approximations to converge. The discrete ordinates method is easy to implement for polarized radiative transfer. However, it requires interpolation and/or linearization of the Stokes vector at an arbitrary polar angle.
The present work presents a simple and efficient method for radiative transfer in an inhomogeneous polarized planar atmosphere illuminated at the top by a collimated sunlight while the bottom reflects the radiation. The method called the discrete spherical harmonics method (DSHM), combines principles of the spherical harmonics method and the discrete ordinates method. First, both the Stokes vector and the phase matrix are expanded in a finite series of the generalized spherical functions. Second, the resulting nonlinear coupled differential equations of spatial variable are expressed in a set of polar angular directions. These equations are then put in the matrix form and solved numerically using the spectral Chebyshev method (STM). The emergent/transmitted radiative intensity and the degree of polarization predictions at any polar direction for both Rayleigh and Mie scattering of homogeneous/inhomogeneous polarized planar atmosphere are predicted and compared with available benchmark literature results.

Radiative transfer equation
In a scattering, absorbing, and emitting plane-parallel atmosphere with polarization effects the transport of the energy is described by the vector radiative transfer equation (VRTE) which can be derived by applying the energy conservation to an element of volume to yield (Chandrasekhar, 1960;Wauben et al., 1993) (Wauben et al., 1993). The degree of polarization of the scattering light is calculated at any point within the atmosphere and at any angular direction by (Stam et al., 2006) ( , , ) = √ ( , , ) 2 + ( , , ) 2 + ( , , ) 2 ( , , ) . (2) The degree of polarization can be viewed as a sum of the degree of linear polarization (Bohren and Clothiaux, 2006), and the degree of circular polarization superposed to give the total scattered wave. As consequence of the illumination, the total radiation in the atmosphere is composed of two parts: the collimated or the direct component (Benassi et al., 1985;Mishchenko, 1990) and the diffuse part, ( , , ), so that one can express the total polarized radiation as Inserting Eq. (6) into Eq. (1) yields to the following equation of the diffuse component In Eq. (7) and (8) the subscript has been dropped down and ( , , ) denotes the diffuse part of the Stokes vector. For > 0, the solution of Eq. (7) is sought under the boundary conditions (Garcia and Siewert, 1986) (0, , ) = 0 (9) and where 0 is the coefficient for Lambertian reflection or the surface albedo and the matrix of reflectance.

The scattering matrix
The phase matrix contains the scattering properties of the atmosphere and describes how the intensity and the state of polarization change when a radiation propagating in an incident direction ( ′ , ′ ), after scattering propagates in the direction ( , ). It is related to the scattering matrix (Θ) by (Hovenier, 1969;Sievert, 1982;Collin et al., 2014) where the scattering angle Θ is the angle between the incoming and the outgoing beam defined as (Korkin, 2013: Collin et al., 2014 cos Θ = ′ + √1 − 2 √1 − ′2 cos( − ′ ).
The matrices ( − ) and (− ) are that needed to rotate the meridian plane before and after the scattering into a local scattering plane. Thus, scattering matrix relates the Stokes parameters of incident beam and scattered beam, defined with respect to their scattering plane. The matrix of rotation is defined in terms of the rotation angle, , by (Hovenier, 1969;Korkin, 2013) For a distribution of particles that scatter independently the phase matrix depends on the radius , the scattering angle and may be obtained by the Mie theory. In general, this matrix has 16 nonzero elements. These elements are expressed in terms of the amplitude of the scattered light and can be found in detail in (Mishchenko et al., 2004). For some assumptions about the distribution of particles, it is possible to reduce the number of nonzero elements of the matrix.
Assuming that in the atmosphere particles are of the same kind, are randomly oriented and have at least one plane of symmetry or particles and their mirror plane in the same number ( Van de Hulst, 1981) this matrix takes the form ( Van de Hulst, 1981;Mishchenko, 1990) . (14) The calculation of these matrix elements are found in (Hansen and Travis, 1974;van de Hulst, 1981). In the independent scattering the transformation matrix is integrated over all particles in the distribution. More generally, in the radiative transfer analysis, it is preferable to use the normalized scattering matrix than the non-normalized one. The normalization constant is determined so that (Muñoz and Mills, 2015) 1 4 ∫ 1 Ω 4 = 1.
Denoted here as P, the dimensionless elements of the normalized scattering matrix for a unit volume are calculated as (Hansen and Travis, 1974) where = 1,2, … 6, the wavelength of the illumination, the radius of the particle, ( ) the particle size distribution per unit volume, 1 and 2 the smallest and the largest size of the distribution, = 2 ⁄ the size parameter of the particle and the scattering coefficient For the Mie scattering and in the case of isotropic and homogeneous spherical or polydispersed spherical particles, 1 = 5 , 3 = 6 (Evans and Stephens, 1991); the scattering matrix does not depend on the azimuthal angle ( Van de Hulst, 1981). For the molecular scattering and the scattering by particles small compared to the wavelength of the illumination, Eq. (15) reads (Hovenier, 1971;Mischenco et al., 2006) (Θ) = 3 4 [ 1 + cos 2 Θ −1 + cos 2 Θ 0 0 −1 + cos 2 Θ 1 + cos 2 Θ 0 The phase matrix defined by Eq. (11) is expanded in a Fourier series as (Garcia and Siewert, 1986) In Eq. (18) the matrix, ( , ′ ), can be expanded in series of generalized spherical functions as (Garcia and Siewert, 1986) where ℓ is the normalized associated Legendre polynomials, ℓ and ℓ are defined in terms of the normalized spherical functions. The normalized Legendre polynomials are defined in terms of the associated Legendre polynomials ℓ ( )as (Dave, 1975) The upper limit of Eq. (22) It should be noted that the properties of the generalized spherical functions, the recursion formulae for computing them and the relations for the computation of the Greek constants are given in detail in (Siewert, 1982;De Rooij, 1985).

Discrete spherical harmonics transform of the VRTE
Following the expansion of the phase matrix given in Eq. (18), the Stokes vector is expanded in Fourier series of the azimuthal angle as (Garcia and Siewert, 1986) The coefficients of expansion define the vector radiance moment. Inserting Eq. (18) and (26) into Eq. (7) and (8) The upper limit of Eq. (29) is chosen according to Evans (1993). The odd number L refers to the order of the discrete spherical harmonics method.
In the DSHM hypothesis Eq. (22) and (29)  In calculating the integral term in the right hand side of Eq. (27) the following orthogonal relation has been used (Siewert and McCormick, 1993) and l = ℓ ℓ /(2ℓ + 1) .
The second main point of the DSHM transform of the VRTE consists to consider Eq. (30) for each polar direction as where − ≤ ≤ , and integer. In order to have a system of ( + 1) coupled differential equations in ( + 1) unknowns in vector moments, the value of must be chosen in such a way that 2 = + 1.
In the same manner, the boundary conditions for the resolution of Eq. (46) read and The solution of Eq. (46), (48) and (49) allows to determine the Chebyshev coefficients and therefore the radiance moments vector. Consequently, the total intensity and the Stokes vector characterizing the polarization state in the atmosphere can easily be deduced.

Numerical results
Two test problems are considered in this paper for the DSHM predictions. were also used for the spatial discretization of the problems. The homogeneous atmosphere corresponds to the case = 0. This later case is used in this work to show that the DSHM is well implemented and can be used for the radiative transfer predictions in polarized atmosphere.
The first problem considered is that of molecular scattering by the atmosphere. The scattering is conservative and the Rayleigh scattering law is given in table 1. The atmosphere of optical thickness = 0.5 is bounded at the bottom by a non-reflecting surface. Figure 1 presents at the boundaries of the atmosphere the scattering intensity of the polarized radiation for different values of the azimuthal angle. Since the fourth Stokes parameter is zero for molecular scattering, calculation indicates that at azimuthal angles of 0° and 180° the third component vanishes. In this situation, the axes of the ellipse defining the plane of polarization are along the meridian plane. Figure 1 shows that the scattering polarized intensity is greatly affected by the value of the constant . The lowest value of the radiance is observed for = 1.
As approaches zero this value increases, the corresponding curve is closed to the curve of = 0 and the atmosphere behaves as a homogeneous slab medium. This aspect is fundamental to demonstrate that the DSHM presented in this work was well implemented. Moreover, a comparison between the DSHM predictions with the results obtained by Natraj et al. (2009) using a method based on the X and Y functions for homogeneous atmosphere also shows the validity of the method. Another aspect to verify the accuracy of the DSHM is the calculation of the relative difference in intensity which is less than 0.5% for | | ≥ 0.1 except for the directions close to zero where a maximum relative error of less than 3% is found. The last graph at the left of figure 1 informs one on how the radiative intensity is affected by the azimuthal angle.
The second example used in this paper is the scattering by polydispersed spherical particles. The atmosphere of optical thickness = 1 is illuminated at the top by a solar radiation at = 0.951 . In addition, the atmosphere is bounded at the bottom by the ground of surface albedo 0 = 0.1. The size of the particles is described by a gamma distribution with the effective radius = 0.2 , the effective variance = 0.07 and the refractive index of = 1.44 (Deirmendjian, 1969). The Mie theory is used to compute the Greeks constant for this problem and the results are listed in table 2. These coefficients agree to ten decimal places to those reported by Vestrucci and Siewert (1984). The constant scattering albedo is slightly modified to consider the absorption by the atmosphere. Thus the scattering albedo used for computation of the results is 0 = 0.99. Figure 2 presents at the boundaries of the atmosphere the scattering intensity of the polarized radiation for different values of the azimuthal angle. In figure 3 the degree of polarization of the scattered light is also presented. As in the Rayleigh scattering, figure 2 shows that the scattered intensity is greatly affected by the value of the constant . The lowest values of the intensity are obtained for = 1, the value at which the degree of polarization shown in figure 3 is highest. As approaches zero the intensity increases, the degree of polarization decreases, the corresponding curve is closed to the curve of = 0 and the atmosphere behaves as a homogeneous slab medium. This result is a first point indicating that the DSHM can be applicable for radiative transfer. Note that at the azimuthal angle = 0 and = , the fourth Stokes component vanishes and the degree of polarization is reduced to the degree of linear polarization. Garcia and Siewert (1989) used the method for solving this problem for the homogeneous atmosphere ( = 0). The DSHM predictions are compared to their results and a good agreement is found. Once more, the relative difference is calculated to check how far the DSHM results are deviated from the FN results. The absolute value of the maximum deviation is less than 0.9% for | | ≥ 0.1 in the radiative intensity and less than 0.8% in the degree of polarization. For directions → 0 a maximum deviation less than 10% and 7% is observed respectively in radiative intensity and the degree of polarization. This result is not surprising since the SHM needs high order to converge for the directions close to zero.

Conclusion
The radiative transfer through a plane-parallel inhomogeneous polarized atmosphere is studied using a discrete spherical harmonics method. In the development of the method, the Stokes vector and the phase matrix are expanded in a finite series of the generalized spherical functions. Then, the resulting first order differential equations of vector moments are expressed in a set of discrete polar directions. The combination of the spherical harmonics and the discrete ordinates principles make the proposed method useful for radiative intensity and the degree of polarization predictions at any polar direction without using an interpolation function or linearization as it is generally encountered in discrete ordinates method. The resulting equations from the discrete spherical harmonics ordinates were put in a matrix form and the Chebyshev spectral method was adopted for its spatial discretization. Accuracy of the method was tested in the case of homogeneous atmosphere for both Rayleigh and Mie scattering. It was found that for a low order of the discrete spherical harmonics method, the emergent/transmitted radiative intensity and the degree of polarization predictions agree well with literature benchmark results except for polar directions close to zero where high order of angular discretization is required for accurate predictions.