Markov chain formalism for polarized light transfer in plane-parallel atmospheres , with numerical comparison to the Monte Carlo method

Building on the Markov chain formalism for scalar (intensity only) radiative transfer, this paper formulates the solution to polarized diffuse reflection from and transmission through a vertically inhomogeneous atmosphere. For verification, numerical results are compared to those obtained by the Monte Carlo method, showing deviations less than 1% when 90 streams are used to compute the radiation from two types of atmospheres, pure Rayleigh and Rayleigh plus aerosol, when they are divided into sublayers of optical thicknesses of less than 0.03. ©2011 Optical Society of America OCIS codes: (030.5620) Radiative transfer; (290.4210) Multiple scattering; (290.5850) Scattering, particles; (290.5855) Scattering, polarization. References and links 1. L. W. Esposito, and L. L. House, “Radiative transfer calculated by a Markov chain formalism,” Astrophys. J. 219, 1058–1067 (1978). 2. L. W. Esposito, “An „adding‟ algorithm for the Markov chain formalism for radiation transfer,” Astrophys. J. 233, 661–663 (1979). 3. J. E. Hansen, and L. D. Travis, “Light scattering in planetary atmospheres,” Space Sci. Rev. 16(4), 527–610 (1974). 4. K. F. Evans, and G. L. Stephens, “A new polarized atmospheric radiative transfer model,” J. Quant. Spectrosc. Radiat. Transf. 46(5), 413–423 (1991). 5. Q. L. Min, and M. Duan, “A successive order of scattering model for solving vector radiative transfer in the atmosphere,” J. Quant. Spectrosc. Radiat. Transf. 87(3-4), 243–259 (2004). 6. S. Y. Kotchenova, E. F. Vermote, R. Matarrese, and F. J. Klemm, Jr., “Validation of a vector version of the 6S radiative transfer code for atmospheric correction of satellite data. Part I: path radiance,” Appl. Opt. 45(26), 6762–6774 (2006). 7. P. W. Zhai, Y. Hu, C. R. Trepte, and P. L. Lucker, “A vector radiative transfer model for coupled atmosphere and ocean systems based on successive order of scattering method,” Opt. Express 17(4), 2057–2079 (2009), http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-17-4-2057. 8. J. H. Hannay, “Radiative transfer: exact Rayleigh scattering series and a daylight formula,” J. Opt. Soc. Am. A 26(3), 669–675 (2009). 9. F. Weng, “A multi-layer discrete-ordinate method for vector radiative transfer in a vertically-inhomogeneous, emitting and scattering atmosphere–I. Theory,” J. Quant. Spectrosc. Radiat. Transf. 47(1), 19–33 (1992). 10. R. J. D. Spurr, “VLIDORT: A linearized pseudo-spherical vector discrete ordinate radiative transfer code for forward model and retrieval studies in multilayer multiple scattering media,” J. Quant. Spectrosc. Radiat. Transf. 102(2), 316–342 (2006). 11. F. M. Schulz, K. Stamnes, and F. Weng, “VDISORT, an improved and generalized discrete ordinate method for polarized (vector) radiative transfer,” J. Quant. Spectrosc. Radiat. Transf. 61(1), 105–122 (1999). 12. J. V. Dave, “Intensity and polarization of the radiation emerging from a plane-parallel atmosphere containing monodispersed aerosols,” Appl. Opt. 9(12), 2673–2684 (1970). 13. J. W. Hovenier, “Symmetry relations for scattering of polarized light in a slab of randomly oriented particles,” J. Atmos. Sci. 26(3), 488–499 (1969). 14. C. E. Siewert, “On the phase matrix basic to the scattering of polarized light,” Astron. Astrophys. 109, 195–200 (1982). #136522 $15.00 USD Received 14 Oct 2010; revised 1 Dec 2010; accepted 20 Dec 2010; published 7 Jan 2011 (C) 2011 OSA 17 January 2011 / Vol. 19, No. 2 / OPTICS EXPRESS 946 15. J. W. Hovenier, and C. V. M. van der Mee, “Basic Relationships for Matrices Describing Scattering by Small Particles”, in Light Scattering by Nonspherical Particles, pp. 61–85, M. Mishchenko, J. W. Hovenier and L. Travis, (eds.), Academic Press, San Diego, 2000. 16. A. B. Davis and F. Xu, “Monte Carlo modeling of polarized light transfer in vertically varying plane-parallel atmospheres, with application to lofted aerosol layer detection using O2 spectroscopy,” J. Quant. Spectrosc. Radiat. Transf. 17. A. A. Kokhanovsky, V. P. Budak, C. Cornet, M. Duan, C. Emde, I. L. Katsev, D. A. Klyukov, S. V. Korkin, L. C-Labonnote, and B. Mayer, “Benchmark results in vector atmospheric radiative transfer,” J. Quant. Spectrosc. Radiat. Transf. 111(12-13), 1931–1946 (2010). 18. G. Marchuk, G. Mikhailov, N. Nazaraliev, R. Darbinjan, B. Kargin, and B. Elepov, The Monte Carlo Methods in Atmospheric Optics, Springer-Verlag, New-York, NY (1980). 19. K. F. Evans, and A. Marshak, “Numerical Methods,” in 3D Radiative Transfer in Cloudy Atmospheres, pp. 243– 281, A. Marshak and A. B. Davis (eds.), Springer, Heidelberg, Germany (2005). 20. A. A. Kokhanovsky, J. L. Deuzé, D. J. Diner, O. Dubovik, F. Ducos, C. Emde, M. J. Garay, R. G. Grainger, A. Heckel, M. Herman, I. L. Katsev, J. Keller, R. Levy, P. R. J. North, A. S. Prikhach, V. V. Rozanov, A. M. Sayer, Y. Ota, D. Tanré, G. E. Thomas, and E. P. ZegeA. A. Kokhanovsky, J. L. Deuzé, D. J. Diner, O. Dubovik, F. Ducos, C. Emde, M. J. Garay, R. G. Grainger, A. Heckel, M. Herman, I. L. Katsev, J. Keller, R. Levy, P. R. J. North, A. S. Prikhach, V. V. Rozanov, A. M. Sayer, Y. Ota, D. Tanré, G. E. Thomas, and E. P. Zege, “The intercomparison of major satellite aerosol retrieval algorithms using simulated intensity and polarization characteristics of reflected light,” Atm. Meas. Techn. 3, 909–932 (2010). 21. D. J. Diner, J. C. Beckert, T. H. Reilly, C. J. Bruegge, J. E. Conel, R. A. Kahn, J. V. Martonchik, T. P. Ackerman, R. Davies, S. A. W. Gerstl, H. R. Gordon, J. Muller, R. B. Myneni, P. J. Sellers, B. Pinty, and M. M. Verstraete, “Multi-angle Imaging Spectro-Radiometer (MISR) instrument description and experiment overview,” IEEE Trans. Geosci. Rem. Sens. 36(4), 1072–1087 (1998). 22. P. Y. Deschamps, F. M. Breon, M. Leroy, A. Podaire, A. Bricaud, J. C. Buriez, and G. Seze, “The POLDER mission: Instrument characteristics and scientific objectives,” IEEE Trans. Geosci. Rem. Sens. 32(3), 598–615 (1994). 23. M. I. Mishchenko, B. Cairns, G. Kopp, C. F. Schueler, B. A. Fafaul, J. E. Hansen, R. J. Hooker, T. Itchkawich, H. B. Maring, and L. D. Travis, “Accurate monitoring of terrestrial aerosols and total solar irradiance: Introducing the Glory mission,” Bull. Am. Meteorol. Soc. 88(5), 677–691 (2007). 24. D. Crisp, R. Atlas, F. Breon, L. Brown, J. Burrows, P. Ciais, B. Connor, S. Doney, I. Fung, and D. Jacob, “The Orbiting Carbon Observatory (OCO) mission,” Adv. Space Res. 34(4), 700–709 (2004). 25. V. Natraj, R. J. D. Spurr, H. Boesch, Y. Jiang, and Y. Yung, “Evaluation of errors in neglecting polarization in the forward modeling of O2 A band measurements from space, with relevance to CO2 column retrieval from polarization sensitive instruments,” J. Quant. Spectrosc. Radiat. Transf. 103(2), 245–259 (2007). 26. V. Natraj, and R. J. D. Spurr, “A fast linearized pseudo-spherical two orders of scattering model to account for polarization in vertically inhomogeneous scattering-absorbing media,” J. Quant. Spectrosc. Radiat. Transf. 107(2), 263–293 (2007).


Introduction and overview
The polarization state of atmospheric radiation potentially contains abundant information about aerosol properties.A computationally efficient and accurate method for vector radiative transfer is therefore important for practical retrieval algorithms.For this purpose, the Markov chain method was developed to calculate the diffusely reflected intensity from a vertically homogeneous and inhomogeneous atmosphere [1][2].Application to Venus" atmosphere indicates its higher computation efficiency than that of adding/doubling method [3][4] by one magnitude [2].Moreover, the Markov chain method retains (i) the advantage of the adding/doubling method of easy physical interpretation and high numerical accuracy in calculating the atmosphere of large optical depth, (ii) the advantage of the successive orders of scattering method in separating the contribution of a specific scattering order numerically [5][6][7] or analytically [8], and (iii) the advantage of the discrete source method [9][10][11] for the internal radiation field calculation.These merits form the motivation of the current research to extend the Markov chain method to polarized radiative transfer in a vertically inhomogeneous atmosphere.
This study is organized as follows.In order to ease the derivation of vector radiative transfer formalism, Section 2 gives an introduction to the Markov chain formalism for the scalar case in a physical way.The formalism differs to that in Ref [1]. in two respects: (1) inclusion of transmission, and (2) consideration of a vertically inhomogeneous atmosphere.Then as a generalization, Section 3 extends the Markov chain formalism to vector radiative transfer in a vertically inhomogeneous plane-parallel atmosphere.A comparison of numerical results calculated by the Markov chain method to the counterparts from the Monte Carlo method is presented in Section 4, followed by a summary and outlook in Section 5.

Scalar radiative transfer
Assume that a vertically inhomogeneous plane-parallel atmosphere of optical depth τ 0 is illuminated by a monochromatic solar flux πF 0 .The atmosphere is divided into N sublayers so that each one has its optical thickness Δτ n = τ n+1 τ n (1  n  N with τ N+1 = τ 0 ), single scattering albedo ω 0 (n) and phase function P(cosΘ,n) which is normalized in such a way that its integral over all directions is 4π.The scattering angle Θ is defined by where u 0 = cosθ 0 and u = cosθ denote the direction of the incident and scattered light, respectively, with respect to the downward increasing optical depth, 0  and  describe the azimuthal plane of the incident and diffuse light.When polarization is neglected, the phase function is even function about ( 0   ) and thus can be expanded into the Fourier series components (P (m) ) for cosine modes only: (0) ( ) 0 0 0 0 0 1 (cos , ) ( , , ; ) ( , ; ) 2 ( , ; ) cos ( ).
m m P n P u u n P u u n P u u n m To further specify the scattering type to be reflection or transmission according to the sign of u and specify the incidence to be upwelling or downwelling according to the sign of u 0 , we define μ = |u| and μ 0 = |u 0 |, getting the following four regimes [3]: * 0 0 0 0 ( , , ; ) ( , , ; ), * 0 0 0 0 ( , , ; ) ( , , ; ), where the subscripts "r" and "t" of the phase function denote the reflection and transmission, respectively for the illumination from above while the asterisk (*) in the superscript denote the illumination from below.Accordingly, the notations for the Fourier series components become ( ) ( ) 00 ( , ; ) ( , ; ), ( ) ( ) 00 ( , ; ) ( , ; ),

Initial scattering in the atmosphere
We begin the calculation of multiple scattering with the intensity of the first order scattered light.This is the initial state for the Markov chain.Taking into account the exponential attenuation of the incident solar light (direction cosine μ 0 ) and the outgoing scattered light (direction cosine μ i ), the contribution of first-order scattering to the diffusely reflected intensity from the upper boundary of the n th layer (τ = τ n ) has the Fourier series () and the single-scattering contribution to the diffused transmitted intensity through the lower boundary of the n th layer (τ = τ n+1 ) has the Fourier series () Note that for integrating the contribution of the incident light in all directions μ i to the scattered intensity in one direction in the next order of scattering, a quadrature weight w i is introduced for μ i on the interval 0  μ i  1.Assuming the emergent single scattering intensity from the boundary to be produced by uniform source distribution in each sublayer, the source intensity can be calculated through dividing Eqs. ( 12)-( 13) by an average attenuation factor c (n,i) along the path Δτ n /μ i , namely, ( ) ( ) ( , ),0, / ( , ),0, / ( , ) U / , where

Intermediate scattering in the atmosphere
To analyze the multiple scattering, we assume the light last scattered from the layer n in the direction μ i travels upward to the layer n' for the next order of scattering.The reduced flux for scattering in n' th layer is: where on the basis of the uniform source distribution in each layer the first term stands for the intensity of the light leaving the upper boundary of n th layer of optical thickness Δτ n , the second term indicates the remaining light intensity reaching the lower boundary of n' th layer after traveling an optical depth (τ n τ n'+1 ), and the last term indicates the attenuated light flux in the n' th layer.Equation ( 16) has the following analytical form: In a similar way, we can calculate the attenuated light flux in the n' th layer for the downward travelling light (n < n') In addition to the upward and downward propagation, partial light flux remains in the same layer (n = n') for scattering: ( , ),( , ) 0 Knowing the reduced flux, the scattered intensity from n' th layer in all directions can be calculated from its Fourier series which satisfies Note that Q represents the transition of light between two states due to one order of scattering.It forms the basis for a later description of all intermediate multiple scattering processes.

Scattering to the boundary
The emergent intensities from the top of the atmosphere (TOA, or τ = 0) and from the bottom of the atmosphere (BOA, or τ = τ 0 ) are contributed by both downward-and upward-traveling light.Assuming they come from layer n' in direction j, and are redirected to emerge in direction μ e through scattering, we then analyze their contribution to the radiation fields at the BOA and TOA.

Downwelling light
Downwelling light from n' th layer partially leaves the TOA through diffuse reflection and partially leaves the BOA through diffuse transmission.The diffuse scattering occurs in all layers n with n  n' so that the Fourier series of the emergent intensities are: N mm e n j n e n j nn RR      (23) for diffuse reflection at the TOA and N mm e n j n e n j nn TT      (24) for diffuse transmission at the BOA, where the subscripts of R and T on the left hand side of the above two equations mean the light leaving the atmosphere in the direction μ e is contributed by the scattering of the light coming from n' th layer in direction μ j , and the subscript (n, e) of each series term on the right hand side of the two equations means the scattering into the direction μ e occurs in n th layer.
The contribution of scattering in n th layer to the emergent intensities is calculated by integrating over τ bounded by τ n and τ n+1 .In terms of the Fourier series we have , where the kernels of the above integrals ( ), ,( , ),( , ) Substitution of Eq. ( 27) and (28) into Eq.( 25) and ( 26 In the particular case that diffuse scattering occurs in the layer n = n', Eq. ( 25) and ( 26) become Substitution of Eq. ( 41) and (42) into Eq.( 39) and ( 40 In the case where diffuse scattering occurs in the layer (n = n'), Eq. ( 39) and (40) become

Radiative transfer through multiple scattering
On the basis of the Π vector and the Q, R, and T matrices, the radiative transfer through different orders of scattering can be achieved in matrix form after invoking the orthogonality relation [12] for the trigonometric functions, namely, where "I", "J", "K", and "L", in the subscripts denote the state of the light described by the layer number (n) and the propagating direction (μ i ); "0" denotes the initial light state in the atmosphere and "e" denotes the light leaving the atmosphere in direction μ e .Change of the state of light from "I" to "J" is represented in the sequence "JI".The first term on the right side of the above equation represents the diffusely emergent light experiencing two scattering events, the second term represents the emergent light experiencing three scattering events, and so on.In the form of an infinite series summation, all orders of multiple scattering are accounted for.If a specific order of scattering is not explicitly expressed and counted, Eqs. ( 49)-( 50) have an equivalent but more compact form: where E is the identity matrix. (

2.5) Total radiation field
The above description has not included the contribution of single-scattered light to the total emergent intensity.For a vertically inhomogeneous atmosphere, the contribution of single scattering to diffuse reflection and transmission can be analytically calculated by summarizing the contribution of scattering of incident solar light in all sublayers, namely, and the contribution to the diffuse transmission by each layer is Knowing the contribution of single and multiple scattering, the total intensities of the diffusely reflected light from the TOA and the transmitted light through the BOA are obtained: Recall that this expression is for the diffuse (single and multiply scattered) light.In the case of transmission through the BOA, the directly transmitted component in the direction (µ 0 ,  0 ), which can be estimated from F 0 , τ 0 , and µ 0 , has to be included for the total field.Our development in this section is parallel to that of Refs [1]- [2], reproduces their results, and extends their formalism to diffuse transmission through a vertically inhomogeneous atmosphere.

Vector radiative transfer
In analogy to the phase function used in scalar radiative transfer, the radiative transfer of polarized light is described by a 4 × 4 Mueller matrix 00 ( , , , )

P
which is defined as the multiplication of the phase matrix P(cosΘ,n) with two rotational matrices about angles i 1 and i 2 [13], namely for diffuse transmission.For m = 0 only the cosine mode is needed in calculation so that

Illustrative numerical calculations
The Markov chain vector 1D radiative transfer model described in the previous sections was programmed in MatLab® for ease of matrix and vector manipulations.Results from this prototype were compared with results from 1D vector Monte Carlo models developed simultaneously in FORTRAN, and a selection of such comparisons are presented in this section.

Implementation of the Monte Carlo method
We use two Monte Carlo models that are described in detail by Davis and Xu [16].The only difference between them is in the representation of vertical variability.One code takes exactly the same input information as the Markov chain model, i.e., N layers with their respective single scattering albedo and phase function; it then tracks histories in optical units, assuming the layers are internally uniform.This version of the Monte Carlo model was benchmarked extensively against the high-precision values provided by Kokhanovsky et al. [17] for a Rayleigh case, an aerosol case, and a cloud case, all single uniform layers over an absorbing surface.
The other code tracks trajectories, where convenient, in physical units and treats spatial variability continuously but in a specific parameterized representation; for the present study, the atmosphere is composed of an exponential distribution of (Rayleigh scattering) molecules and a uniform layer of particles, e.g., polydisperse (Mie scattering) spheres.
In all of the following, 10 8 Monte Carlo histories were traced using the code that ingests discrete layers; this leads to relative uncertainties on the order of 10 4 for well-sampled signals but this precision deteriorates for poorly populated signals (typically, the 3 rd and 4 th components of the Stokes vector).Radiances were computed using the local estimation technique [18,19] at the same ordinates as used by the Markov chain model in the principal and perpendicular planes.

Separate orders of scattering in a pure Rayleigh atmosphere
We start our comparisons of the Markov chain and Monte Carlo methods in predicting various orders of scattering by a homogeneous atmosphere of pure Rayleigh scattering with optical depth τ 0 = 0.5 and single scattering albedo ω 0 = 1.The atmosphere is divided into 20 layers.In Markov chain formalism, contributions of single and higher-order scattering are calculated by Eqs. ( 54 61), respectively.In the Monte Carlo method, tallying the various orders of scattering is trivial to encode.An illustrative calculation of the first two, three, then all, orders of scattering contributing to the diffusely reflected and transmitted fields in the principal plane (- 0 = 0° and 180°) is given in Figs.1-4, showing agreement to be better than 99%.In these figures, only the I and Q components are plotted since the U and V components vanish in the principal plane.Moreover, positive viewing zenith angles are used for the - 0 = 0° plane while their negative counterparts are used for - 0 = 180°.Note that, although the Markov chain formalism is presented for the vertically inhomogeneous atmosphere in the current paper, by setting the Mueller matrix P(cosΘ,n) and single scattering albedo ω 0 (n) to be constant for all sublayers, the formalism reduces to its counterpart for the homogeneous atmosphere [1].
This aerosol model was used in a recent intercomparison of aerosol retrieval algorithms described in Ref [20].and the aerosol optical depth corresponds to their "Case 12." The phase matrix for the aerosol distribution is calculated by integrating the contribution of particles over the size interval 0  a  30 μm, which is discretized to 40000 quadrature points.To highlight molecular scattering, we selected the MISR blue wavelength of 446.4 nm [21], leading to the stated Rayleigh optical depth.
The atmosphere is divided into 42 sublayers, with the optical thickness of each layer as well as the Rayleigh scattering fraction in each layer are listed in Table 1.This layering approximates in optical distance space the combination of an exponential decay of the Rayleigh column with a scale height of 8 km with a uniform aerosol layer between 0 and 2 km.The surface is black and solar illumination is at zenith angle θ 0 = 60°.This atmospheric structure, as well as the lower-boundary and illumination conditions, were prescribed in Ref [20].for 443 nm albeit at a slightly different wavelength (hence Rayleigh optical depth and aerosol phase matrix) and without molecular depolarization effects.The layers in Table 1 were generated automatically from top to bottom under the sole condition that none of them has an optical thickness exceeding 0.03, which is taken in Markov chain method as a reasonable rule-of-thumb to ensure the single scattering domination in each layer.
Figures 5-6 give the I and Q components of the reflection field in the principal plane (- 0 = 0° and 180°) while Figs.7-8 give the I and Q components of the transmission field in the same plane.The U and V components vanish identically in the principal plane.In the perpendicular plane - 0 = ±90°, however, all components of the reflected and transmitted fields appear, as demonstrated in the panels of Figs.9-16.Using 45 Gauss-Legendre ordinates on the interval 0 < μ < 1 (equivalently, 90 streams), Figs.5-8 and Figs.9-16 show the disagreement between the Markov chain and Monte Carlo methods to be less than 1%.Using finer (coarse) sublayer thickness and more (less) Gauss-Legendre ordinates on μ"s interval can further reduce (increase) the disagreement, but costing (saving) in computation time.For example, the maximum error of Q for diffuse reflection and transmission in the principal and perpendicular planes increases to about 4% when 23 Gauss-Legendre ordinates are used but 75% computation time is saved.

Conclusion and outlook
Accurate yet efficient modeling of vector (polarized) polarized radiative transfer is becoming increasingly important, largely to keep up with instrument development in several application areas: medical imaging, ocean and atmospheric optics (including remote sensing), planetary missions, etc. Space-based Earth observation sensors with polarization capability include various incarnations of the POLDER (POLarization and Directionality of Earth Reflectances) sensor, most recently on ESA"s PARASOL satellite [22], and the APS (Aerosol Polarimetry Sensor) on NASA"s upcoming Glory mission [23].Even when polarization is not a focus of the mission, it needs to be accounted for to keep the accuracy of the forward radiative transfer modeling at par with instrument precision.For instance the Orbiting Carbon Observatory (OCO) mission [24] targets global column amounts of CO 2 based on differential optical absorption spectroscopy (DOAS) in the near-IR.This calls for DOAS in both CO 2 and O 2 bands, and the later is significantly affected by aerosol scattering, therefore by polarization [25].To remedy this situation without sacrificing computational efficiency, Natraj and Spurr [26] developed an approximation model based on two orders of scattering that compensates for polarization effects left out by OCO"s operational scalar code.
Since the Markov chain formalism has inherent computational efficiency in handling all orders of scattering from a vertically inhomogeneous atmosphere, it forms another attractive solution to polarized radiative transfer.Comparison of the diffusely reflected and transmitted radiation fields calculated by Markov chain method to those by Monte-Carlo method shows the deviation to be less than 1%, depending on the optical depth and quadrature mesh.On the basis of the current work, we plan to further extend the Markov chain approach to the cases concerning non-uniform and/or bi-directional polarized surface reflectivity, spherical shell geometry, and more structurally complex 3D atmospheres.
of the incident light from the n th layer in direction μ i (for the second order of scattering () m B = [ ( from the n'-th layer and direction μ j , which is contributed by the incident light from n th layer and direction μ i , The scattered light becomes the incident light in the next order of scattering.For integrating the contribution of the incident light in all directions to the next order of scattering, a Gaussian quadrature weight w j is introduced to define a transition factor () ( , ),( , )

1 ,
scattering of the downwelling (denoted by the superscript "d") incident light from n' th layer to the diffusely emergent light through reflection (denoted by the subscript "r") and transmission (denoted by the subscript "t") from/through n th layer, respectively.Denoting ( ), 0 md I to be the Fourier series vector for all states of the incident light intensity from n' th layer, multiplication with matrices ( ), series of the reflected and transmitted fields, respectively, where ( ),

36) 2 . 3 . 2
Upwelling lightUpwelling light from n' th layer leaves the TOA through diffuse transmission and partially leaves the BOA through diffuse reflection.The diffuse scattering occurs in all layers n with n  n' so that the Fourier series of the emergent intensities are, over τ of the n th layer bounded by τ n scattering of the upwelling (denoted by the superscript "u") incident light from n' th layer to the diffusely emergent light through transmission (denoted by the subscript "t") and reflection (denoted by the subscript "r), from/through the n th layer, respectively.Denoting ( ), 0 mu I to be the Fourier series vector for all states of the incident light intensity from n' th layer, multiplication with the matrices ( ), of Eq. (33) into Eq.(46) and Eq.(34) into Eq.( n th term on the right hand side of Eqs.(62)-(65) means the contribution of diffuse light after experiencing (n + 1) scattering events.Difference of the vector and scalar radiative transfer case in calculating the Π vector and the Q, R, and T matrices is expressed here.For calculating the contribution of the downwelling incident light to the diffuse field at the TOA and BOA, each element of the R-matrix in Eqs.(62) and (64) becomes a 4×4 cell which is calculated by Eq. (35) after replacing the phase function series()

Fig. 1 .
Fig. 1.The first two, three, then all, orders of scattering contributing to the diffuse reflection field (I component) in the principal plane for a homogeneous Rayleigh atmosphere of optical depth τ0 = 0.5 and unit single scattering albedo.Error bars are based on the standard deviation among 10 realizations of the Monte Carlo simulation, each with 10 7 histories (cf.Figures 12 and 16 for a situation with degraded precision).

Fig. 2 .
Fig. 2. Same as Fig. 1 but the Q component of the diffuse reflection field is plotted.

Fig. 3 .
Fig. 3. Same as Fig. 1 but the I component of the diffuse transmission field is plotted.

Fig. 4 .
Fig. 4. Same as Fig. 1 but the Q component of the diffuse transmission field is plotted.

4. 3
Vertically stratified aerosol-and-Rayleigh atmosphereAs another example of numerical calculation, we choose an Earth-like atmosphere of total optical depth τ 0 = 1.22934.The atmosphere is composed of Mie aerosols (optical depth 1.0074) of refractive index n = 1.38+10 8 i and Rayleigh molecules (optical depth 0.2286).For the aerosols, a lognormal distribution with parameters s = 1 and a 0 = 0.1 μm is assumed:

Fig. 5 .
Fig. 5. Diffuse reflection field (I component) in the principal plane of the vertically inhomogeneous atmosphere described in the main text under solar illumination at 60° incidence.As in Fig. 1, the error bars are based on the standard deviation among 10 realizations of the Monte Carlo simulation, each with 10 7 histories.

Fig. 6 .
Fig. 6.Same as Fig. 5 but the Q component of the diffuse reflection field is plotted.

Fig. 7 .
Fig. 7. Same as Fig. 5 but the I component of the diffuse transmission field is plotted.

Fig. 8 .
Fig. 8. Same as Fig. 5 but the Q component of the diffuse transmission field is plotted.

Fig. 9 .
Fig. 9. Same as Fig. 5 but the I component of the diffusely reflected field in the perpendicular plane (0 = ±90°) is plotted.

Fig. 10 .
Fig. 10.Same as Fig. 9 but the Q component of the diffusely reflected field is plotted.

Fig. 11 .
Fig. 11.Same as Fig. 9 but the U component of the diffusely reflected field is plotted.

Fig. 12 .
Fig. 12. Same as Fig. 9 but the V component of the diffusely reflected field is plotted.

Fig. 13 .
Fig. 13.Same as Fig. 9 but the I component of the diffusely transmitted field is plotted.

Fig. 14 .
Fig. 14.Same as Fig. 9 but the Q component of the diffusely transmitted field is plotted.

Fig. 15 .
Fig. 15.Same as Fig. 9 but the U component of the diffusely transmitted field is plotted.

Fig. 16 .
Fig. 16.Same as Fig. 9 but the V component of the diffusely transmitted field is plotted.