Issue 
A&A
Volume 573, January 2015



Article Number  A72  
Number of page(s)  20  
Section  Numerical methods and codes  
DOI  https://doi.org/10.1051/00046361/201424042  
Published online  22 December 2014 
Preconditioned backward Monte Carlo solutions to radiative transport in planetary atmospheres
Fundamentals: Sampling of propagation directions in polarising media^{⋆}
^{1}
ESA Fellow, ESA/RSSD, ESTEC, 2201 AZ
Noordwijk, The
Netherlands
email: tonhingm@gmail.com
^{2}
Research School of Physics and Engineering and Fenner School of
Environment and Society, Australian National University, ACT
0200
Canberra,
Australia
email: frank.mills@anu.edu.au
^{3}
Space Science Institute, Boulder, CO
80301,
USA
Received: 21 April 2014
Accepted: 15 August 2014
Context. The interpretation of polarised radiation emerging from a planetary atmosphere must rely on solutions to the vector radiative transport equation (VRTE). Monte Carlo integration of the VRTE is a valuable approach for its flexible treatment of complex viewing and/or illumination geometries, and it can intuitively incorporate elaborate physics.
Aims. We present a novel preconditioned backward Monte Carlo (PBMC) algorithm for solving the VRTE and apply it to planetary atmospheres irradiated from above. As classical BMC methods, our PBMC algorithm builds the solution by simulating the photon trajectories from the detector towards the radiation source, i.e. in the reverse order of the actual photon displacements.
Methods. We show that the neglect of polarisation in the sampling of photon propagation directions in classical BMC algorithms leads to unstable and biased solutions for conservative, opticallythick, strongly polarising media such as Rayleigh atmospheres. The numerical difficulty is avoided by preconditioning the scattering matrix with information from the scattering matrices of prior (in the BMC integration order) photon collisions. Preconditioning introduces a sense of history in the photon polarisation states through the simulated trajectories.
Results. The PBMC algorithm is robust, and its accuracy is extensively demonstrated via comparisons with examples drawn from the literature for scattering in diverse media. Since the convergence rate for MC integration is independent of the integral’s dimension, the scheme is a valuable option for estimating the diskintegrated signal of stellar radiation reflected from planets. Such a tool is relevant in the prospective investigation of exoplanetary phase curves. We lay out two frameworks for disk integration and, as an application, explore the impact of atmospheric stratification on planetary phase curves for large starplanetobserver phase angles. By construction, backward integration provides a better control than forward integration over the planet region contributing to the solution, and this presents a clear advantage when estimating the diskintegrated signal at moderate and large phase angles.
Key words: polarization / radiative transfer
A oneslab, planeparallel version of the PBMC algorithm is available at the CDS via anonymous ftp to cdsarc.ustrasbg.fr (130.79.128.5) or via http://cdsarc.ustrasbg.fr/vizbin/qcat?J/A+A/573/A72
© ESO, 2014
1. Introduction
The gases and aerosols that make up a planetary atmosphere leave characteristic signatures on the radiation emitted and/or reflected from the planet. The technique of polarimetry utilises the polarisation state of emergent radiation to investigate the planet’s atmospheric optical properties. Polarimetry is relevant in the remote sensing of planetary atmospheres both as a standalone technique and in combination with photometry. In the solar system, polarimetric observations made from spaceborne and groundbased telescopes have yielded insight into the gas and aerosol envelopes of Earth (Dollfus 1957; Hansen & Travis 1974), Venus (Coffeen 1969; Hansen & Hovenier 1974), Mars (Santer et al. 1985), Jupiter and Saturn (Morozhenko & Yanovitskii 1973; Schmid et al. 2011; West et al. 1983), Titan (Veverka 1973; West & Smith 1991), and Neptune and Uranus (Joos & Schmid 2007; Michalsky & Stokes 1977; Schmid et al. 2006).
Various spacecraft for Earth (ADEOS I and II, PARASOL) and solar system exploration (e.g. Voyager, Galileo, Cassini) carried instrumentation with (limited) polarimetric capabilities. Most modern groundbased observatories are equipped with polarimeters for either spectroscopy or imaging. Groundbased observations of the outer planets, however, only have partial coverage of the SuntargetEarth phase angle, which limits the possible physical insight from polarimetric investigations. For the above reasons, it is generally agreed that polarimetry’s potential for characterising the atmospheres of Earth and the rest of the solar system planets remains underexploited. Interestingly, the discovery of planets orbiting stars other than our Sun has caused a renewed interest in polarimetry both as a detection and a characterisation technique. The key idea behind this new interest is that stars are typically unpolarised or weakly polarised, whereas planets may be partially polarised, which presents an advantage for separating the planet from the glare of its host star (e.g. Seager et al. 2000; Stam et al. 2004).
The new born field of exoplanet research is prompting significant effort in the development of polarimetric facilities, as demonstrated by proposed space missions such as ESA’s SPICES (Boccaletti et al. 2012), dedicated instrumentation for Gemini (Macintosh et al. 2006), or ESO’s Very Large Telescope and EuropeanExtremely Large Telescope (Beuzit et al. 2008; Kasper et al. 2008). Correspondingly, on the theoretical front, there has been work to investigate polarimetry’s potential for identifying planets’ orbital parameters, as well as for characterising their main atmospheric and surface features (e.g. Bailey 2007; Fluri & Berdyugina 2010; Seager et al. 2000; Stam and collaborators 2004; 2008, but also Karalidi & Stam 2012, Karalidi et al., 2011; 2012; 2013; Williams & Gaidos 2008; Zugger et al. 2010; 2011). Since the number of exoplanets already surpasses the number of solar system planets, theoretical investigations that explore gas, cloud, and surface properties, possibly in the framework of a new generation of general circulation models, will continue to play a key role in the prediction and prospective characterisation of exoplanetary observables.
This paper is devoted to numerical modelling of radiation scattered by planetary atmospheres. Our approach relies on backward Monte Carlo (BMC) integration of the vector radiative transport equation (VRTE). Special attention is paid to the sampling of propagation directions in polarising media. We show that in classical BMC integration, failing to account for polarisation in the sampling of propagation directions may destabilise and bias the numerical solution in conservative, optically thick, strongly polarising media. We propose a preconditioned BMC (PBMC) algorithm and show that preconditioning the scattering matrix with information from prior collisions (in the order of backward integration) eliminates the numerical difficulties. Preconditioning is equivalent to providing information about the history and polarisation state of photons through their simulated trajectories. We describe in detail the algorithm and its performance. Because it consistently delivers precisions of 10^{4} when compared to solutions that are accurate to at least that level, the algorithm may be considered “exact” (in the de Haan et al. 1987 sense) or nearly so. This paper is part of an ongoing effort to build a tool for efficiently simulating the radiation emerging from both diskresolved and diskintegrated realistic planetary atmospheres. In its scalar form, the algorithm has already been used without description (García Muñoz & Pallé 2011; García Muñoz & Mills 2012; García Muñoz et al. 2011; 2012; 2014). The cases investigated here focus on Rayleigh and Mie scattering, for which the scattering matrix is easy to obtain. The theory is more general than that and should also apply to scattering particles with different scattering matrices.
The paper is structured as follows. In Sect. 2, we note some of the differences between forward and backward integration. BMC algorithms are very selective with the planet regions that they probe, and this is a clear advantage, for instance, when producing the diskintegrated signal from a planet at a specified phase angle. We review the fundamentals of BMC algorithms and discuss the sampling of photon propagation directions in classical BMC algorithms and in our PBMC approach. We also present two different schemes for integrating the net radiation reflected from a sphericalshell planet. In Sect. 3, we assess the performance of the classical and preconditioned algorithms with test cases for plane parallel configurations. In Sect. 4, we predict a number of planetary phase curves. The extensive suite of test cases considered will hopefully help guide the decision of potential users of the PBMC algorithm. Finally, in Sect. 5 we summarise the main conclusions and comment on followup work.
2. The BMC algorithm
MC algorithms for radiative transport fall within the general class of Markov chain methods for the statistical simulation of photon collisions in scattering media (Cashwell & Everett 1959; Marchuk et al. 1980). By using appropriate statistical estimators, MC algorithms can estimate the radiation within and emerging from a medium.
MC algorithms are classified as forward or backward (FMC and BMC, respectively), depending on whether the solution is built by simulating the photon trajectories from the radiation source towards the observer or vice versa. FMC algorithms easily account for the photon’s polarisation state in the sampling of the photon propagation direction following a collision (e.g. Bartel & Hielscher 2000; Bianchi et al. 1996; Cornet et al. 2010; Fischer et al. 1994; Hopcraft et al. 2000; Kastner 1966; Schmid 1992; Whitney 2011). That is not immediately possible in BMC algorithms because the scattering events are treated in the reverse order that they actually occur. BMC algorithms generally treat the sampling of propagation directions by omitting the radiation’s polarisation state and correcting subsequently for the bias introduced (Collins et al. 1972 and works thereafter; e.g. Emde et al. 2010; Gay et al. 2010; Oikarinen 2001). As shown below, that approach may fail to render accurate solutions in conditions for which the scattering directions of the photons are strongly influenced by their polarisation states.
MC algorithms are exact in the sense (de Haan et al. 1987) that their accuracy is in principle only limited by the number of photon trajectory simulations. Thus, MC algorithms are often used as standards in the validation of other methods, particularly in cases that involve complex viewing and/or illumination geometries (Loughman et al. 2004; Postylyakov 2004).
BMC algorithms are better suited for problems with small detectors and large radiation sources, and the opposite is true for FMC algorithms (Modest 2003). This important distinction means that BMC integration turns out to be the appropriate choice for numerous applications in investigating planetary atmospheres. By tracing the photon trajectories from the detector towards the planet (or towards a part of the planet that is known to be illuminated), BMC algorithms offer a more efficient approach to achieving the desired accuracy. This is not directly possible in the FMC framework because there is no previous knowledge about the directions the photons will take to exit the medium. In FMC algorithms, moreover, estimating the emergent radiation typically requires averaging over a range of exiting directions. (Alternatively, variance reduction techniques such as the nextevent pointestimator can be utilised, e.g. Kaplan et al. 2001 and Lux & Koblinger 1985; their efficiency however is strongly dependent on the detector’s acceptance angle.) These characteristics penalise the computational efficiency of FMC algorithms, especially when only a specified number of viewing geometries with narrow acceptance angles are of interest.
Additional properties that make FMC/BMC algorithms appealing in their application to planetary atmospheres include:

They are easy to implement and debug. Their description can indeed be accomplished in less than one page (see Appendix A).

Implementing the scattering matrix does not require a series expansion of the matrix elements.

Curvature and twilight effects are naturally accounted for. Limbviewing geometries do not require any special treatment.

It is easy to separate the contributions from the atmosphere and surface, from different atmospheric layers, or from various orders of scattering.

Scattering by large particles, which lead to highly asymmetric scattering phase functions, can be treated without significant computational penalty.

The computational cost for solving the VRTE and its scalar counterpart are comparable.

The accuracy of the solution depends on the number of photon trajectory simulations. Moderately accurate solutions can be obtained at low computational costs.

In BMC algorithms, each photon collision can be utilised to estimate the contribution to the detector from various incident directions of the illuminating source.
Our implementation of the algorithm follows the basic layout by O’Brien (1992; 1998), which we have extended to include polarisation. The implementation makes use of variance reduction techniques, which logically arise from the mathematical elaboration of the integrals that occur in the formal solution to the VRTE. O’Brien (1992; 1998) provides an excellent introduction to these ideas, and we follow the nomenclature in those works to a large extent.
2.1. Fundamentals
Our interest lies in the VRTE for a scattering and absorbing medium without volume or surface emission sources: (1)where x and s are vectors of position and direction, β(x) and γ(x) are the scattering and extinction coefficients of the medium (independent of direction), and dΩ(s′) is the differential solid angle about direction s′. The ratio ϖ(x) = β(x)/γ(x) is the local single scattering albedo of the medium. In terms of the θ and φ angles at the top of Fig. 1, dΩ(s′) = sinθdθdφ. I(x,s) = [ I,Q,U,V ] ^{T} is the Stokes vector that describes the polarisation state of radiation, and ℙ(x,s,s′) is a 4 × 4 matrix for deflection of radiation from the incident direction s′ to the emergent direction s. ℙ(x,s,s′) = (π − i)(x,s,s′)(−i′), and (π − i) and (− i′) are rotation matrices for the conversion of the Stokes vector from the meridional plane (the plane formed by the z axis of a userdefined rest reference frame and the direction of photon propagation) to the scattering plane and vice versa. The rotation matrix is (2)with κ either π − i or −i′, angles i′ and i defined as sketched in (Fig. 1, top), and (x,s,s′) is the scattering matrix, for which we assume (x,s,s′) = (x,s·s′) = (x,cosθ). The normalisation of (x,s·s′) verifies (3)In Mie scattering theory for spherical particles, the matrix is fully prescribed by means of four elements (Mishchenko et al. 2002): (4)In the Rayleigh limit for particle sizes that are much smaller than the radiation wavelength, the four elements take on analytical expressions that, neglecting anisotropy effects, are a_{1} = 3(1 + cos^{2}θ) / 4, b_{1} = 3(−1 + cos^{2}θ) / 4, a_{3} = 3cosθ/ 2, and b_{2} = 0.
Fig. 1 Top: definition of the incident, s′, and emergent, s, photon directions at a scattering event within the atmosphere. The xyz axes form a rest reference frame fixed to the planet. The differential solid angle dΩ(s′) = sinθdθdφ is defined with s serving as polar axis. Angles θ ∈ [ 0,π] and φ ∈[0, 2π]. In the backtracing of photons of BMC algorithms, s is known at each collision and s′ must be sampled from the relevant scattering phase function. Angles i′ and i, both ∈[0, π], are needed for consistent referencing of the Stokes vector throughout the scattering process. Vectors {, and } and {e_{1}, e_{2} and e_{3}} define righthanded coordinate systems at the meridional planes of the incident and emergent photon directions, respectively. Bottom: definition of the incident, s′, and emergent, s, photon directions at a reflection event at the local surface (plane ). Here, n is the inwardpointing normal vector at the surface, and ẑ is oriented along −n. The differential solid angle dΩ(s′) = sinθdθdφ is defined with ẑ serving as polar axis. 
Equation (1) admits the formal solution for the Stokes vector at {x_{k},s_{k}}: (5)On the righthand side, the first term stands for radiation reflected from a point x_{kb} at the boundary of the integration domain into direction s_{k}, whereas the second term represents the radiation scattered within the medium from {x_{ka},s_{ka}} to {x_{k},s_{k}}. Each term may include both diffuse and unscattered radiation components, defined as the contributions from photons that have undergone at least one and zero prior scattering collisions, respectively, and dℓ_{ka} stands for the arclength along the path joining x_{kb} and x_{k}. The transmittance between x_{ka} and x_{k} is (6)and t(x_{k},x_{kb}) is defined analogously.
It is useful to introduce the dimensionless variables (7)a(x_{k},x_{kb}) = 1 − t(x_{k},x_{kb}), and ϖ(x_{ka}) = β(x_{ka}) /γ(x_{ka}), that, by construction, range from 0 to 1, so that the formal solution to Eq. (1) becomes (8)To evaluate Eq. (8), boundary conditions at the top and bottom of the atmosphere are needed. We consider here that the only source of illumination is stellar radiation from direction s_{⊙}, for which the unimpeded, unpolarised irradiance is F_{⊙} = π [ 1,0,0,0 ] ^{T}δ(s′ − s_{⊙}), with ∫m walodΩ(s′)δ(s′ − s_{⊙}) ≡ 1. Tacitly, the given F_{⊙} assumes that the stellar size subtended from the planet is small so that the radiation incident on the planet is oriented in a single direction s_{⊙}. We furthermore assume Lambert reflection with albedo r_{g} at the atmospheric bottom (the planet’s surface) and a transparent atmospheric top for outgoing radiation.
The surface reflectance properties relate I(x_{kb},s_{k}) to the incident Stokes vector at the boundary I(x_{kb},s_{kb}). For Lambert reflection at the atmospheric bottom, (9) Here, n(x_{kb}) is the inwardspointing normal vector at the surface, is the fourbyfour depolarizing matrix with _{1,1} = 1 as the only nonzero entry, and I(x_{kb},s_{kb}) is the Stokes vector for diffuse radiation reaching the surface. The lower panel of Fig. 1 sketches the relevant geometrical parameters for photon collisions at the surface. The two terms of Eq. (9) are the separate contributions to I(x_{kb},s_{k}) from both diffuse radiation and from unscattered stellar radiation reaching the surface. From the definition of the transparent atmospheric top, I(x_{kb},s_{k}) ≡ 0 for x_{kb} at the top of the atmosphere.
Similarly, it is convenient to separate the diffuse and unscattered radiation within the atmospheric medium: (10)In both Eqs. (9) and (10), x_{⊙} is the intersection at the top boundary of the rays traced in the −s_{⊙} direction from x_{kb} and x_{ka}, respectively. Clearly, t(x_{kb},x_{⊙}) and t(x_{ka},x_{⊙}) ≡ 0 if the stellar disk is not visible from either x_{kb} and x_{ka}, respectively.
With the above considerations, Eq. (8) is now expressed as (11)where and In Eq. (11), only the term preceded by a(x_{k},x_{kb}) occurs for x_{kb} at the atmospheric top, but for generality, we retain both of them. Both the ℒ_{ℬ} and terms can be evaluated based on the optical properties of the medium, whereas the ℬ and terms need additional information in the form of the diffuse radiation vectors I(x_{kb},s_{kb}) and I(x_{ka},s_{ka}).
Starting from {x_{0},s_{0}}, which determines the position of and entry direction into the detector, recurrent use of Eq. (11), complemented by Eqs. (12)–(15), leads to an expression for I(x_{0},s_{0}) as an infinite summation series of integrals of increasingly higher dimensions (O’Brien 1992; 1998). Physically, higher dimension integrals account for additional orders of scattering of the simulated photons. Figure 2 shows the definition of the pairs {x_{0a},s_{0a}}, {x_{0b},s_{0b}}, {x_{0aa},s_{0aa}}, {x_{0ab},s_{0ab}}, {x_{0ba},s_{0ba}}, and so forth, which appear in the recurrence law. The series is convergent (and thus amenable to truncation) provided that the optical thickness of the medium is finite and/or the medium is not fully conservative (i.e. either ϖ or r_{g} ≤ 1). The recurrence law builds the solution for I(x_{0},s_{0}) by splitting each summation into a double summation involving new ℬ and integrals at each step.
Appendix A spells out the first few terms in the summation series and summarises the practical implementation in the PBMC algorithm. Rewriting the integrals that appear in the summation in terms of appropriately normalised variables leads to improved convergence rates, an approach that is equivalent to socalled variance reduction techniques (O’Brien 1992). In what follows, we address the integration in solid angle and the significance of polarisation in the sampling of photon propagation directions, which is the feature unique to our PBMC algorithm with respect to other BMC schemes.
Fig. 2 In black, sketch demonstrating the construction of the {x_{k},s_{k}} pairs starting from {x_{0},s_{0}} as the photon is traced back from the detector through the medium. Vectors are pointed in the direction of photon propagation, which is the reverse of the direction of integration in the BMC algorithm. In principle, a {x_{k},s_{k}} pair can lead to two new {x_{ka},s_{ka}} and {x_{kb},s_{kb}} pairs. In the MC implementation of the algorithm, a scheme based on the coefficients of the ℒ_{ℬ}+ℬ and + operators, Eq. (11), determines whether the photon’s next move occurs within the atmospheric medium or whether the photon moves onto the planet’s surface. In blue, one specific photon trajectory within the family of possible trajectories. For this specific trajectory, the red arrows denote the direction of the unscattered stellar photons. 
2.2. Monte Carlo integration
The essence of MC integration is to estimate multidimensional integrals through evaluation of the integrand at properly selected values of the integration variables: (16) Here, each j represents a random draw from the uniform distribution functions u_{k}∈[0, 1]. Importantly, MC integration converges to the exact value at a rate that depends on the number of realisations, N, but not on the dimension of the integral, d.
In a BMC framework, the evaluation of the summation series for I(x_{0},s_{0}) is interpretable in terms of photons whose trajectories are simulated in the backwards direction, i.e. from the detector through the medium, finally reaching the radiation source. Thus, we regularly refer to the determination of the {x_{k},s_{k}} pairs as simulated photon trajectories made up of collision events at x_{0a}, x_{0b}, x_{0aa}, x_{0ab}, x_{0ba}, etc. Ultimately, the solution to the VRTE is built by simulating a number n_{ph} (=N in Eq. (16)) of photon trajectories.
2.3. Integration in solid angle
The summation series for I(x_{0},s_{0}) obtained from recurrent use of Eq. (11) contains multidimensional integrals in solid angles: (17)for collisions within the atmospheric medium. For simplicity in the notation, we removed all references to x_{0k} within the ℙ matrices. The treatment of collisions at the bottom boundary is analogous. In a BMC framework, dΩ(s′) integration at a particular collision event entails selecting an incident s′ direction for a given emergent s direction (see Fig. 1), according to an appropriate probability density function.
2.3.1. The classical sampling scheme
In classical BMC algorithms (Collins et al. 1972, and thereafter), evaluation of Eq. (17) proceeds by separating it into (18)and, subsequently, sampling the θ and φ angles in each integral from the local _{1,1} function and from a uniform distribution between 0 and 2π, respectively. Tacitly, the sampling scheme assumes that the relative orientations between s and s′ must depend on the local properties of the medium but not on the propagation history of the photons, or that any bias introduced by proceeding that way can be subsequently corrected for by dividing by the sampled _{1,1}. The assumption is exact in the treatment of the scalar RTE, but is fundamentally erroneous in polarising media. We refer to the simplified approach based on Eq. (18) as the classical sampling scheme for photon propagation directions.
2.3.2. The preconditioned sampling scheme
A more appropriate approach to the evaluation of Eq. (17) is to sample
(19)By proceeding sequentially, at each step all the involved photon propagation directions but the one being sampled are known. The scheme derives directly from Eq. (17), and preserves the history of the simulated photon trajectories through the ordered arrangement of the products of ℙ matrices. At each collision event, the matrices ℍ(s_{0},s_{0a}) = (≡unity matrix), ℍ(s_{0a},s_{0aa}) = ℙ(s_{0},s_{0a}), ℍ(s_{0aa},s_{0aaa}) = ℙ(s_{0},s_{0a})ℙ(s_{0a},s_{0aa}), ..., effectively precondition the local ℙ matrix, and in turn, the probability that scattering occurs in any of the possible s_{0a}, s_{0aa}, s_{0aaa}, ..., incident directions. The preconditioning matrix evolves as the photon trajectory is being backtraced, and in this way, the photon history is preserved throughout the simulation. Hereafter, we term this approach the preconditioned sampling scheme for photon propagation directions. This scheme is at the core of our PBMC algorithm.
Expanding Eq. (19) yields insight into the preconditioned sampling scheme. For an arbitrary ℍℙ(s,s′)dΩ(s′) = ℍ(π − i)(x,θ)(−i′) dΩ(θ,φ), the (1, 1) entry leads to an expression proportional to f(θ,φ)dθdφ= (20)where we defined q = ℍ_{1,2}/ ℍ_{1,1} and u = ℍ_{1,3}/ ℍ_{1,1}. In the derivation of Eq. (20), we used the geometrical relation between i and φ, for φ locally defined with respect to the meridian plane (see Fig. 1, Top). Angle i′ is evaluated once both θ and φ are determined.
Two important properties apply to f(θ,φ), namely: [1] it is ≥0 for θ∈[0, π] and φ∈[0, 2π]; and [2] its integral over the θ–φ domain is equal to one, which is straightforward for confirming from the normalisation of Eq. (3). Since a_{1}(θ) ≥ 0 and  b_{1}  ≤  a_{1}  (Mishchenko et al. 2002), property [1] requires that  qcos(2φ) − usin(2φ)  ≤ 1. To prove that condition, it suffices to show that the first row of ℍ, [ ℍ_{1,1},ℍ_{1,2},ℍ_{1,3},ℍ_{1,4} ] (=ℍ_{1,1}[ 1,q,u,v ] in our own notation), forms from (21)The vector resulting from Eq. (21) is indeed the transpose of (22)which is the Stokes vector for an associated direct problem of photons propagating onwards from the detector. In this direct problem, the relevant scattering matrix is ^{T}. For Mie scattering, Eq. (4), the matrix satisfies ^{T}(b_{2}) = (−b_{2}), which suggests a connection with the adjoint formulation based on vector Green’s functions proposed by Carter et al. (1978). From the association of the backward problem with its direct counterpart of scattering matrix ^{T}, it becomes apparent that q and u are relative linear polarisations and v is the corresponding relative circular polarisation. As a result,  qcos(2φ) − usin(2φ)  ≤ 1.
Thus, f(θ,φ) is a bivariate probability density function that can be used to sample the propagation directions in the backtracing of photons. Our preconditioned scheme of Eq. (19) is indeed similar in structure to the schemes utilised in some FMC algorithms (e.g. Bartel & Hielscher 2000; Bianchi et al. 1996; Cornet et al. 2010; Fischer et al. 1994; Hopcraft et al. 2000; Kastner 1966; Schmid 1992; Whitney 2011).
In practice, the sampling is facilitated by separating f(θ,φ)= f_{θ}(θ)f_{φ  θ}(φ  θ), with Here, f_{θ} is the conventional θsampling function implemented in most FMC and BMC algorithms, whether treating the scalar or vector RTE. Function f_{φ  θ}(φ  θ) conveys that sampling in φ is constrained by θ and, through q and u, also by the photon polarisation state and history. Figure 3 explores f(θ,φ) for a few combinations of q and u ≡ 0 in the specific case of a Rayleigh medium. The classical sampling scheme is equivalent to drawing the θ and φ from the probability density function f(θ,φ;q ≡ 0). Doing so appears inappropriate in strongly polarising media where b_{1}(θ) /a_{1}(θ), q and u may take absolute values close to one through the photon simulations. The consequences of this are investigated below.
Fig. 3 Probability density function f(θ,φ)=f_{θ}(θ)f_{φ  θ}(φ  θ) in the precondioned sampling scheme of photon propagation directions, Eqs. (23)–(24), for a Rayleighscattering medium. Note the changes in f(θ,φ) with q, especially near the maximum of  b_{1}(θ) /a_{1}(θ)  for θ = 90°. By ignoring polarisation, the classical sampling scheme determines the θ and φ values of the incident propagation direction from f(θ,φ;q ≡ 0). It is apparent that the classical sampling scheme is more likely to fail in strongly polarising media that involve high q values during the backtracing of photons. 
2.4. Diskintegration schemes
We are interested in the radiation emerging from both diskresolved and diskintegrated planetary atmospheres. We derive two diskintegration schemes here and describe their incorporation into the PBMC algorithm.
2.4.1. Integration over the “visible” disk
Horak (1950) laid out the expressions for evaluating the diskintegrated radiation scattered from a planet over its “visible” disk. In this context, “visible” refers to the disk portion that appears illuminated by singlescattered photons as viewed from the observer’s vantage point. We refer to the sketch of Fig. 4, which presents the relevant geometrical parameters. Provided that both the observer and the star are sufficiently far from the planet, Horak (1950) arrives at the expression: (25)which we adapt to the vector case by using the Stokes vector I, in which case F (=[F_{I},F_{Q},F_{U},F_{V}]^{T}) is the irradiance Stokes vector. Here, ρ (=R_{p}+h_{TOA} for planets with a solid core of radius R_{p} and an atmosphere extending up to altitudes of h_{TOA}) and Δ are the radius of the planet’s scattering disk and the observertoplanet distance, respectively. We normalise F by eliminating the (ρ/ Δ)^{2} factor from Eq. (25).
Fig. 4 Geometrical parameters relevant to the integration over the planet’s “visible” disk. For a pair of u_{i} and v_{i} values, N is the location on the disk where the observer’s line of sight intercepts the atmosphere. N is equivalent to x_{0} in our implementation of the PBMC algorithm. 
Rather than working with longitudes, ζ_{d}, and colatitudes, η_{d}, it is convenient to introduce the two auxiliary variables: such that, after some manipulations, Eq. (25) transforms into (28)The premultiplying factor before the double integral is the projected size of the planet’s “visible” disk. The double integral may be seen as an average radiance Stokes vector over that domain.
In the form of Eq. (28), it is straightforward to insert the evaluation of F into the PBMC algorithm as the sum (29)where u_{i} and v_{i} are picked from the random uniform distributions u, v ∈ [ 0, 1]. Each u_{i}, v_{i} yields the location on the planet’s disk where the observer’s line of sight intercepts the planet’s atmosphere or, equivalently, x_{0} in the implementation of the algorithm of Appendix A. For a sufficiently remote observer, s_{0} is, according to Fig. 4, permanently oriented along the x axis. The application of Eq. (29) requires the inversion of Eqs. (26), (27). For u → η_{d}, we interpolate from precalculated tabulations of u = u(η_{d};α). For v → ζ_{d}, the inversion is done analytically. Since in our formulation I is by default referenced to the meridian plane containing the z axis and s_{0}, and s_{0} is fixed in space, there is no need to rotate the emergent ⟨ I(u_{i},v_{i}) ⟩ Stokes vectors, which can be added directly into Eq. (29). In our normalisation, the first of the F elements is A_{g}Φ(α), with A_{g} being the planet’s geometric albedo and Φ(0) ≡ 1.
2.4.2. Integration over the entire disk
Alternatively to the integration over the “visible” disk, one can proceed by integrating over the entire disk. Introducing r and Θ as the polar coordinates that determine the projection of N in Fig. 4 on the yz plane, and the normalised variables u′ = Θ / 2π and v′ = (r/ρ)^{2}, integration over the projected surface element rdrdΘ leads to (30)which, after eliminating the (ρ/ Δ)^{2} factor, translates into (31)in the PBMC algorithm. Again, , are picked from uniform distributions u′, v′∈[0, 1].
Some of the advantages of the latter implementation with respect to the one in Sect. 2.4.1 include [1] it makes no assumption on the extent of the effectivelyscattering disk and, therefore, properly handles the full range of phase angles from superior to inferior conjunctions; [2] each photon trajectory simulation can simultaneously contribute to various specified phase angles. A drawback of the latter implementation (shared with FMC algorithms) is that for a given number of photon realisations n_{ph}, the solution statistics becomes poorer for the larger phase angles because fewer of the simulated photon trajectories actually connect the observer and the direction of illumination. We explore in Sect. 4 some of these issues in application of the two diskintegration schemes to both Rayleigh and Venuslike atmospheres.
3. Comparison of the PBMC algorithm against solutions from other methods
We assessed the performance of our PBMC algorithm against a suite of test cases for which reliable solutions are either available in the literature or can be produced with existing models. The suite includes solutions to both the scalar and vector RTE, different viewing/illumination geometries, and a variety of optical properties for the scattering particles. Here in Sect. 3, we focus on scattering in planeparallel atmospheres.
Figure 5 sketches the relevant angles. In particular, an azimuth of 0 corresponds to the observer facing the Sun, and 180° to the observer looking away from the Sun. In the scalar RTE calculations, polarisation is omitted by zeroeing all entries but _{1,1} in the scattering matrix.
Fig. 5 Illumination and viewing angles for the planeparallel atmosphere test cases discussed in Sect. 3. 
3.1. Nonpolarised Rayleigh scattering
In a first assessment, we compared our PBMC algorithm in its scalar mode against DISORT (Stamnes et al. 1988) solutions in Rayleigh scattering media. The exercise includes 15 876 test cases that explore both optically thin and thick atmospheres with viewing/illumination angles from zenith inclination to nearly horizontal pointing (see Table 1). The comparison, the details of which are given in the appendices, shows an excellent match between the two approaches.
Parameters in the investigation of conservative Rayleigh scattering in planeparallel atmospheres.
3.2. Polarised Rayleigh scattering
Coulson et al. (1960) tabulated solutions for the elements of the Stokes vector in conservative, polarising, Rayleighscattering atmospheres above Lambert reflecting surfaces. More recently, Natraj and collaborators (2009; 2012) have extended the calculations to arbitrarily large optical thicknesses. The newly tabulated Stokes vectors (that we adopt as reference) are claimed to be accurate to within one unit in the eighth decimal place. We computed the 15 876 cases summarised in Table 1 for n_{ph} up to 10^{7} with our PBMC algorithm in its VRTE mode. For comparison, we utilised both the classical and preconditioned sampling schemes introduced in Sect. 2.3.
Figure 6 shows δI (=(I_{BMC}−I_{ref}) /I_{ref} × 100) for the preconditioned (top) and classical sampling schemes (bottom). For the latter, Fig. 7 shows δP (=(P_{BMC}−P_{ref}) /P_{ref} × 100), where .
The δI graphs reveal that the two sampling schemes generally perform well for optical thicknesses ≤4, but that the classical scheme destabilises and/or biases the solutions for larger thicknesses. A similar behaviour also occurs for δP. Median values for  δI  as calculated with the preconditioned scheme are listed in Table C.1. For the PBMC solutions, the convergence rate is comparable to that for the solution of the scalar RTE.
Fig. 6 Differences in intensity, δI, for the solution of the VRTE in conservative, Rayleighscattering atmospheres. The algorithm uses the preconditioned (top) and classical (bottom) sampling schemes for photon propagation directions. 
Fig. 7 Differences in polarisation, δP, for the solution of the VRTE in conservative, Rayleighscattering atmospheres. The algorithm uses the preconditioned sampling scheme for photon propagation directions. 
Figure 8 offers some insight into the stability issue with the classical sampling scheme. It shows the convergence history for the I Stokes element for a cos(OPA) = cos(SPA) = 1 viewing/illumination geometry and varying optical thicknesses above a black surface. (OPA/SPA stands for observer/solar polar angle, Table 1.) The most striking feature of Fig. 8 is that the classical sampling scheme produces abrupt changes in the solution with effects that may not go away even after many photon simulations. The instabilities become more frequent and noticeable for the larger optical thicknesses. Referring to Eq. (24) and Fig. 3, the neglect of polarisation in the classical sampling scheme is likely to favour some propagation directions rather than others and, in turn, erroneously bias the solution. Inspection of some of the abrupt changes indicates that they are associated with a sequence of photon collision events each with scattering angle θ near 90° and therefore likely to be misrepresented by the classical sampling scheme. The disturbance becomes more apparent in optically thick, conservative media because they allow for many more collisions before the photon is lost. Further evidence for the latter comes from the fact that Rayleigh calculations with ϖ ~ 0.95 or less (not shown) show no stability issues for any optical thickness in the range tested. The bottom line is that the primary assumption of the classical sampling scheme, i.e. that the multidimensional integral of Eq. (17) can be approximated by separate integrals as given by Eq. (18) plus a subsequent correction, becomes inappropriate for specific configurations.
Fig. 8 Convergence history of the I Stokes element for a conservative Rayleigh atmosphere over a black surface with cos(SPA) = cos(OPA) = 1 and two different optical thicknesses. Thin and thick curves represent the solutions obtained with the preconditioned and classical sampling schemes, respectively. Abrupt changes in the solution with the classical sampling scheme for optical thickness of 4 occur, but they are not discernible on the scale of the graph. 
The idea is confirmed by investigating the solution to the VRTE in other polarising media. For this purpose, we produced scattering matrices at λ = 0.63 μm for monodisperse droplets of real refractive index equal to 1.53 and a few radii from 1.2 ×10^{1} to 1.7 × 10^{1}μm. Figure 9 shows the corresponding −b_{1}(θ)/a_{1}(θ) ratios, which are properties of the media but also the corresponding degrees of polarisation for photons scattered one single time. When referring to the structure of Eqs. (23), (24), it is apparent that lower b_{1}(θ)/a_{1}(θ) ratios distort the probability density function f(θ,φ;q ≠ 0) less with respect to the case for q = 0. The convergence history for the solutions to the multiple scattering problem in a medium of optical thickness equal to 16 and cos(OPA) = cos(SPA) = 1 are shown in Fig. 10. They reveal that the classical scheme performs poorly in the more strongly polarising media, but performs similarly to the preconditioned sampling algorithm in less polarising conditions.
To the best of our knowledge, there have been no previous reports of difficulties using BMC algorithms with classical sampling, probably because benchmarking solutions for optically thick Rayleigh atmospheres had not been readily available. This example serves to highlight the importance of benchmarking solutions in the literature.
Fig. 9 Polarisation in single scattering for monodisperse droplets of various radii and real refractive index equal to 1.53 at λ = 0.63 μm. 
Fig. 10 Convergence history for scattering by the monodisperse droplets of Fig. 9. Thin and thick curves represent the solutions obtained with the preconditioned and classical sampling schemes, respectively. The classical scheme produces inconsistent solutions for strongly polarising conditions. The calculations assumed optical thickness equal to 16 and cos(OPA) = cos(SPA) = 1. 
Solutions to the Stokes vector in a conservative, hazeL atmosphere of optical thickness equal to 1 and cos(SPA) = 0.5.
3.3. Polarised Mie scattering
We tested our PBMC algorithm against a number of VRTE solutions in Miescattering media. The Stokes vectors for radiation emerging from a conservative atmosphere with socalled hazeL scattering particles have been tabulated by de Haan et al. (1987) from calculations based on the doublingadding method. Tables 2 and B.1 show some of their solutions and the corresponding PBMC calculations. The I Stokes element from both calculation methods generally agrees to the fourth decimal place for n_{ph} = 10^{9}. Typically, solutions accurate to within one per cent in I are obtained for n_{ph} = 10^{5}. Polarisation is low in all cases investigated in Table 2. As a consequence, the convergence of the Q, U and V elements is slower than for I. In the appendices, we extend the comparison by considering the results published by Garcia & Siewert (1986) for scattering within a Venuslike atmosphere. The good match of our PBMC results attests to the capacity of our algorithm to produce accurate solutions to elaborate scattering problems.
4. Planetary phase curves
That the convergence rate of MC integration is independent of the dimension of the integral, Eq. (16), can be used to efficiently estimate the net radiation scattered from the planet. In the solar system, Venus represents a unique demonstration of how diskintegrated polarisation can be used to infer a planet’s cloud composition (Coffeen 1969; Hansen & Hovenier 1974). At remote distances from Earth, exoplanets will not be spatially resolvable in the near future and, thus, their investigation must rely on diskintegrated measurements. Initial attempts to investigate the optical properties of exoplanet atmospheres in reflected light by means of polarisation have been made (Berdyugina et al. 2008; 2011; Wiktorowicz 2009). Foreseeably, a new generation of telescopes and instruments will provide the technical capacity to detect and characterise a variety of exoplanets.
To explore the diskintegration schemes of Sect. 2.4, we utilised a few configurations relevant to both Rayleigh and Venuslike atmospheres. Essentially, the diskintegration scheme selects the entry point of the photon into the atmosphere. The threedimensional photon trajectory is then traced through the medium. The PBMC algorithm is implemented over a spherical shell description of the planet’s atmosphere, which allows us to investigate phenomena related to atmospheric curvature and stratification at large starplanetobserver phase angles.
4.1. Rayleigh phase curves
Buenzli & Schmid (2009) have investigated with an FMC algorithm the phase curves of Rayleighscattering planets. Their study expands on earlier work (e.g. Kattawar & Adams 1971) by systematically exploring the parameter space (optical thickness, atmospheric single scattering albedo, and Lambert surface albedo). Madhusudhan & Burrows (2012) have also produced Rayleigh phase curves on the basis of analytical solutions to the planeparallel problem. Rayleigh scattering may provide a first approximation to the interpretation of a planet’s phase curve. It is, however, of limited usefulness in the general understanding of possibly occurring atmospheres. In such cases, more flexible treatments including Mie and other nonRayleigh forms of scattering are needed. Thus, Stam et al. (2006) have devised an efficient technique for disk integration based on the planeparallel approximation that can deal with arbitrary scattering particles for planets with horizontally uniform atmospheres.
Parameters in the investigation of disk integration for both conservative and nonconservative Rayleigh atmospheres.
We produced Rayleigh phase curves for the configurations listed in Table 3 with the visibledisk integration scheme of Sect. 2.4.1 and compared them to those published by Buenzli & Schmid (2009). Specific properties of the curves such as the geometric or spherical albedo, or the value and position of the polarisation peak have been discussed in that work, and are not discussed here any further. With the PBMC algorithm all properties are evaluated at the specified phase angles without having to bin (and possibly extrapolate) in phase angle.
The agreement between Buenzli & Schmid (2009) and our PBMC calculations is very good. For n_{ph} = 10^{6}, the median of the absolute differences in F_{I} between the two approaches over the α = 7.5–132.5° range is about 0.1%, which is consistent with the accuracy targeted by Buenzli & Schmid (2009). Since the computational time is dictated by the number of photon realisations, it turns out that the computational cost is comparable in both the spatially resolved problems of Sect. 3.2 and in the spatially unresolved problems discussed here. In other words, integrating over the disk involves a computational time comparable to obtaining the solution over a localised region of the planet. With the visibledisk integration scheme, moreover, the convergence properties of the algorithm become independent of phase angle.
Figure E.1 shows phase curves for F_{I} and F_{Q}/F_{I} with ϖ = 1, r_{g} = 1; and τ = 0.1, 0.4, 1, 2, 5, and 10; and n_{ph} = 10^{4}. Because we use the xz meridional plane for referencing the Stokes vector I, the ratio F_{Q}/F_{I} is consistent with positive polarisation in the xz plane perpendicular to the scattering plane. Figure 11 shows the computational times on a 2.6 GHz desktop computer for an average point of a phase curve and n_{ph} = 10^{4}. Computational times depend on the platform, and are not often published in the literature, which prevents a comparison with the performance of other algorithms.
Fig. 11 Computational time per phase curve point for n_{ph}=10^{4} and Rayleighscattering calculations. For a full phase curve with, for instance, 34 evenly separated points from 7.5 to 172.5° (as for Fig. E.1), the computational time is 34 times what is indicated in the plot. The curves correspond to different values of the single scattering albedo of the medium, ϖ, and the Lambert surface albedo, r_{g}. 
4.2. Venus phase curves
Venus is a wellknown example that demonstrates the potential of diskintegrated polarimetry in the remote investigation of clouds (Coffeen 1969; Hansen & Hovenier 1974). The diskintegrated polarisation of Venus is low but sensitive to wavelength and phase angle, facts that were exploited by Hansen & Hovenier (1974) to characterise the droplets that make up the Venus upper clouds. Venus sets a valuable precedent for eventually investigating exoplanetary clouds with polarimetry.
As a further assessment of our PBMC algorithm, we looked into the Venus polarisation phase curves. This analysis has the added value of allowing comparison with real planetary measurements. From the visible through the nearinfrared, the Venus clouds are optically thick and close to fully conservative. The analysis, thus, provides insight into the performance of the PBMC algorithm in conditions that require many photon collisions per simulation.
Fig. 12 Polarisation phase curves for Venus. Black diamonds are the measurements used in Hansen & Hovenier (1974). Colour symbols and curves are our PBMC calculations for n_{ph}=10^{4} and 10^{5}, respectively. This figure can be compared to Figs. 4, 8, 9, and 11, 12 in Hansen & Hovenier (1974). 
Figure 12 shows the digitised data points for the degree of linear polarisation utilised by Hansen & Hovenier (1974) at 0.365, 0.445, 0.55, 0.655, and 0.99 μm. The colour symbols are our PBMC calculations for n_{ph} = 10^{4}, and the underlying solid curves are the calculations for n_{ph} = 10^{5}. For the modelling, we use the prescriptions for particle size distributions (gammadistribution, effective radius r_{eff} = 1.05, effective variance v_{eff} = 0.07), refractive indices, and atmospheric single scattering albedo inferred by Hansen & Hovenier (1974). We assume that the atmosphere is made up of a single slab of optical thickness equal to 30 overlaying a fully reflective Lambert surface. The legends in the panel give additional information about the Rayleighscattering component, f_{R} (see Hansen & Hovenier 1974), which becomes important at UV wavelengths, and various r_{eff} values. For n_{ph} = 10^{4} and the visibledisk scheme of Sect. 2.4.1, the computational time per point in the phase curve is about 8 secs. For the curves of Fig. 12, we took 2°increments in α, which leads to the full phase curve for n_{ph} = 10^{4} being produced in about 12 min. The statistical dispersion of the PBMC calculations is smaller than the dispersion associated with the measurements, and from a practical viewpoint, it seems appropriate to truncate to n_{ph} = 10^{4}. The phase curves of Fig. 12 can be directly compared to the model calculations of Figs. 4, 8, 9, and 11, 12 in Hansen & Hovenier (1974), which confirms the good agreement between both approaches.
In addition, we produced polarisation phase curves at wavelengths from 1.2 to 2.4 μm for the droplets’ size distribution given above and realonly refractive indices based on a 75% H_{2}SO_{4}/H_{2}O solution by mass (Hansen & Hovenier 1974). They are shown in Fig. 13, which further illustrates the sensitivity of polarisation to wavelength.
Fig. 13 Polarisation phase curves for Venus calculated with our PBMC algorithm with =10^{5} at nearinfrared wavelengths for a H_{2}SO_{4}/H_{2}O dilution at 75%. We adopted a single scattering albedo of 1 in all cases. n_{r} is the refractive index in each case. 
As a final exercise, we explored the appropriateness of integrating over the visible disk, Sect. 2.4.1, against the more comprehensive approach of Sect. 2.4.2. For this, we took the Venus atmosphere at 0.55 μm described above as a basis. Differences are expected to arise when the atmosphere is vertically extended. Thus, we stratified the total optical thickness of the atmosphere (=30) with scale heights H (the efolding length for changes in the γ extinction coefficient in the vertical) of 4, 8, 16, and 32 km. Both diskintegration schemes produce nearly identical results (not shown) for the F_{Q}/F_{I} ratio. In contrast, the differences in F_{I}, Fig. 14, can become significant when the planet approaches inferior conjunction, especially for the larger scale heights. Figure 14 provides valuable clues for choosing the appropriate diskintegration scheme for specific applications. The figure also shows that stratification and curvature effects become important for sufficiently large phase angles and H/R_{p} ratios.
Fig. 14 Phase curves F_{I} in the optical for Venuslike planets with atmospheres stratified according to the given scale heights. 
5. Summary and future work
We have presented a novel preconditioned backward Monte Carlo (PBMC) algorithm to solve the vector radiative transport equation (VRTE) in planetary atmospheres. A unique feature of our PBMC algorithm is that it preconditions the scattering matrix before sampling the incident propagation direction at a photon collision. Preconditioning retains some of the information associated with the polarisation state of photons, a feature shown to be critical for correct treatment of conservative, optically thick, strongly polarising media. This is, to the best of our knowledge, the first investigation to report the numerical difficulties that may occur in BMC algorithms (and possibly FMC as well) when polarisation is ignored in sampling propagation directions. We give extensive evidence that our proposed preconditioned sampling scheme ensures the stability of the PBMC algorithm.
We explored the performance of the PBMC algorithm, showing that it consistently produces solutions accurate to better than 0.01% in the Stokes element I when compared to published benchmarks provided that enough photon trajectories are simulated. Our extensive assessment exercise, that includes accuracies and computational times, should help potential users assess the advantages and disadvantages of the method. We believe that similar exercises should become common place in the investigation of VRTE solvers.
In its spherical shell version, our PBMC algorithm is well suited to evaluating the net radiation scattered by a spatially unresolved planet. This feature is particularly interesting in any investigation of the phase curves of solar system planets and exoplanets. We proposed two diskintegration schemes and showed that integration over the “visible” disk incurs a computational cost comparable to solving the VRTE over a localised region of the planet, provided that the spatial details of the emerging radiation can be overlooked. Thus far, we have focused on planetary atmospheres that may be vertically stratified but are otherwise homogeneous in the horizontal direction. Future work will extend the diskintegration scheme to horizontally inhomogeneous planets, thus accounting for a full threedimensional description of the planet. Similar ideas will also be explored for diskintegrated thermal emission and for the simultaneous spectralanddisk integration of both scattered and thermallyemitted radiation.
A oneslab, planeparallel version of our PBMC algorithm is available upon request. Making the algorithm publicly available will hopefully encourage comparative investigations of VRTE solvers.
Acknowledgments
We gratefully acknowledge N. I. Ignatiev for confirming the definition of angles in DISORT, Richard Turner for assistance with the computer cluster “The Grid” at ESA/ESTEC, and the Grupo de Ciencias Planetarias at UPV/EHU, Spain, for access to their computational resources.
References
 Bailey, J. 2007, Astrobiol., 7, 320 [NASA ADS] [CrossRef] [Google Scholar]
 Bartel, S., & Hielscher, A. H. 2000, Appl. Optics, 39, 1580 [NASA ADS] [CrossRef] [Google Scholar]
 Berdyugina, S. V., Berdyugin, A. V., Fluri, D. M., & Piirola, V. 2008, ApJ, 673, L83 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Berdyugina, S. V., Berdyugin, A. V., Fluri, D. M., & Piirola, V. 2011, ApJ, 728, L6 [NASA ADS] [CrossRef] [Google Scholar]
 Beuzit, J.L., Feldt, M., Dohlen, K., et al. 2008, Proc. SPIE, 7014, 18 [Google Scholar]
 Bianchi, S., Ferrara, A., & Giovanardi, C. 1996, ApJ, 465, 127 [NASA ADS] [CrossRef] [Google Scholar]
 Boccaletti, A., Schneider, J., Traub, W., et al. 2012, Exp. Astron., 34, 355 [NASA ADS] [CrossRef] [Google Scholar]
 Buenzli, E., & Schmid, H. M. 2009, A&A, 504, 259 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Carter, L. L., Horak, H. G., & Sandford II, M. T. 1978, J. Comput. Phys. 26, 119 [NASA ADS] [CrossRef] [Google Scholar]
 Cashwell, E. D., & Everett, C. J. 1959, A practical manual on the Monte Carlo method for random walk problems (New York: Pergamon Press) [Google Scholar]
 Coffeen, D. A. 1969, AJ, 74, 446 [NASA ADS] [CrossRef] [Google Scholar]
 Collins, D. G., Blättner, W. G., Wells, M. B., & Horak, H. G. 1972, Appl. Optics, 11, 2684 [NASA ADS] [CrossRef] [Google Scholar]
 Cornet, C., CLabonnote, L., & Szczap, F. 2010, J. Quant. Spectr. Rad. Transf., 111, 174 [NASA ADS] [CrossRef] [Google Scholar]
 Coulson, K. L., Dave, J. V., & Sekera, Z. 1960, Tables related to radiation emerging from a planetary atmosphere with Rayleigh scattering (Berkeley, Los Angeles: University of California Press) [Google Scholar]
 de Haan, J. F., Bosma, P. B., & Hovenier, J. W. 1987, A&A, 183, 371 [NASA ADS] [Google Scholar]
 Dollfus, A. 1957, Supplements aux Annales d’Astrophysique, 4, 3 [Google Scholar]
 Emde, C., Buras, R., Mayer, B., & Blumthaler, M. 2010, Atmos. Chem. Phys., 10, 383 [Google Scholar]
 Fischer, O., Henning, Th., & Yorke, H. W. 1994, A&A, 284, 187 [NASA ADS] [Google Scholar]
 Fluri, D. M., & Berdyugina, S. V. 2010, A&A, 512, A59 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ford, E. B., Seager, S., & Turner, E. L. 2001, Nature, 412, 885 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Garcia, R. D. M., & Siewert, C. E. 1986, J. Quant. Spectr. Rad. Transf., 36, 401 [NASA ADS] [CrossRef] [Google Scholar]
 GarcíaMuñoz, A., & Mills, F. P. 2012, A&A, 547, A22 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 GarcíaMuñoz, A., & Pallé, E. 2011, J. Quant. Spectr. Rad. Transf., 112, 1609 [Google Scholar]
 García Muñoz, A., Pallé, E., Zapatero Osorio, M. R., & Martín, E. L. 2011, Geophys. Res. Lett., 38, L14805 [NASA ADS] [CrossRef] [Google Scholar]
 García Muñoz, A., Zapatero Osorio, M. R., Barrena, R., et al. 2012, ApJ, 755, 103 [NASA ADS] [CrossRef] [Google Scholar]
 García Muñoz, A., PérezHoyos, S., & SánchezLavega, A. 2014, A&A, 566, L1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gay, B., Vaillon, R., & Mengüç, M. P. 2010, J. Quant. Spectr. Rad. Transf., 111, 287 [NASA ADS] [CrossRef] [Google Scholar]
 Hansen, J. E., & Hovenier, J. W. 1974, J. Atmos. Sci., 31, 1137 [NASA ADS] [CrossRef] [Google Scholar]
 Hansen, J. E., & Travis, L. D. 1974, Space Sci. Rev., 16, 527 [NASA ADS] [CrossRef] [Google Scholar]
 Hopcraft, K. I., Chang, P. C. Y., Walker, J. G., & Jakeman, E. 2000, in Light Scattering from Microstructures, eds. F. Moreno, & F. González, Lect. Notes Phys., 534, 135 [Google Scholar]
 Horak, H. G. 1950, ApJ, 112, 445 [NASA ADS] [CrossRef] [Google Scholar]
 Joos, F., & Schmid, H. M. 1980, A&A, 463, 1201 [Google Scholar]
 Kaplan, B., Ledanois, G., & Drévillon, B. 2001, Appl. Optics, 40, 2769 [NASA ADS] [CrossRef] [Google Scholar]
 Karalidi, T., & Stam, D. M. 2012, A&A, 546, A56 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Karalidi, T., Stam, D. M., & Hovenier, J. W. 2011, A&A, 530, A69 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Karalidi, T., Stam, D. M., & Hovenier, J. W. 2012, A&A, 548, A90 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Karalidi, T., Stam, D. M., & Guirado, D. 2013, A&A, 555, A127 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kasper, M. E., Beuzit, J.L., Verinaud, C., et al. 2008, Proc. SPIE, 7015, id. 70151S [Google Scholar]
 Kastner, S. O. 1966, J. Quant. Spectr. Rad. Transf., 6, 317 [NASA ADS] [CrossRef] [Google Scholar]
 Kattawar, G. W., & Adams, C. N. 1971, ApJ, 167, 183 [NASA ADS] [CrossRef] [Google Scholar]
 Loughman, R. P., Griffioen, E., Oikarinen, L., et al. 2004, J. Geophys. Res., 109, 6303 [CrossRef] [Google Scholar]
 Lux, I., & Koblinger, L. 1991, Monte Carlo particle transport methods: Neutron and photon calculations (Boca Raton, Florida: CRC Press Inc.) [Google Scholar]
 Macintosh, B., Graham, J., Palmer, D., et al. 2006, Proc. of the SPIE, 6272, 62720L [Google Scholar]
 Madhusudhan, N., & Burrows, A. 2012, ApJ, 747, 25 [NASA ADS] [CrossRef] [Google Scholar]
 Marchuk, G. I., Mikhailov, G. A., Nazaraliev, M. A., et al. 1980, The Monte Carlo methods in atmospheric optics (Berlin, Heidelberg: SpringerVerlag) [Google Scholar]
 Michalsky, J. J., & Stokes, R. A. 1977, ApJ, 213, L135 [NASA ADS] [CrossRef] [Google Scholar]
 Mishchenko, M. I., Travis, L. D., & Lacis, A. A. 2002, Scattering, absorption and emission of light by small particles (Cambridge: Cambridge University Press) [Google Scholar]
 Modest, M. F. 2003, J. Heat Transfer, 125, 57 [CrossRef] [Google Scholar]
 Morozhenko, A. V., & Yanovitskii, E. G. 1973, Icarus, 18, 583 [NASA ADS] [CrossRef] [Google Scholar]
 Natraj, V., & Hovenier, J. W. 2012, ApJ, 748, 28 [NASA ADS] [CrossRef] [Google Scholar]
 Natraj, V., Li, K.F., & Yung, Y. L. 2009, ApJ, 691, 1909 [NASA ADS] [CrossRef] [Google Scholar]
 O’Brien, D. M. 1992, J. Quant. Spectr. Rad. Transf., 48, 41 [NASA ADS] [CrossRef] [Google Scholar]
 O’Brien, D. M. 1998, J. Quant. Spectr. Rad. Transf., 60, 573 [NASA ADS] [CrossRef] [Google Scholar]
 Oikarinen, L. 2001, J. Geophys. Res., 106, 1533 [NASA ADS] [CrossRef] [Google Scholar]
 Postylyakov, O. V. 2004, J. Quant. Spectr. Rad. Transf., 88, 297 [Google Scholar]
 Santer, R., Deschamps, M., Ksanformaliti, L. V., & Dollfus, A. 1985, A&A, 150, 217 [NASA ADS] [Google Scholar]
 Schmid, H. M. 1992, A&A, 254, 224 [Google Scholar]
 Schmid, H. M., Joos, F., & Tschan, D. 2006, A&A, 452, 657 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Schmid, H. M., Joos, F., Buenzli, E., & Gisler, D. 2011, Icarus, 212, 701 [NASA ADS] [CrossRef] [Google Scholar]
 Seager, S., Whitney, B. A., & Sasselov, D. D. 2000, ApJ, 540, 504 [NASA ADS] [CrossRef] [Google Scholar]
 Stam, D. M. 2008, A&A, 482, 989 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Stam, D. M., Hovenier, J. W., & Waters, L. B. F. M. 2004, A&A, 428, 663 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Stam, D. M., de Rooij, W. A., Cornet, G., & Hovenier, J. W. 2006, A&A, 452, 669 [NASA ADS] [CrossRef] [EDP Sciences] [MathSciNet] [Google Scholar]
 Stamnes, K., Tsay, S.C., Jayaweera, K., & Wiscombe, W. 1988, Appl. Optics, 27, 2502 [Google Scholar]
 Veverka, J. 1973, Icarus, 18, 657 [NASA ADS] [CrossRef] [Google Scholar]
 West, R. A., & Smith, P. H. 1991, Icarus, 90, 330 [NASA ADS] [CrossRef] [Google Scholar]
 West, R. A., Hart, H., Hord, C. W., et al. 1983, J. Geophys. Res., 88, 8679 [NASA ADS] [CrossRef] [Google Scholar]
 Whitney, B. A. 2011, Bull. Astron. Soc. India, 39, 101 [NASA ADS] [Google Scholar]
 Wiktorowicz, S. J. 2009, ApJ, 696, 1116 [NASA ADS] [CrossRef] [Google Scholar]
 Williams, D. M., & Gaidos, E. 2008, Icarus, 195, 927 [NASA ADS] [CrossRef] [Google Scholar]
 Zugger, M. E., Kasting, J. F., Williams, D. M., Kane, T. J., & Philbrick, C. R. 2010, ApJ, 723, 1168 [NASA ADS] [CrossRef] [Google Scholar]
 Zugger, M. E., Kasting, J. F., Williams, D. M., Kane, T. J., & Philbrick, C. R. 2011, ApJ, 739, 12 [NASA ADS] [CrossRef] [Google Scholar]
Appendix A: PBMC algorithm implementation
The implementation of our PBMC algorithm follows the formulation by O’Brien (1992; 1998), which we extend to include polarisation. Other BMC formulations exist, that differ mainly in their definition of the statistical estimator’s kernel (e.g. Collins et al. 1972; Marchuk et al. 1980; Postylyakov 2004). For completeness, we sketch the practical details of the algorithm here.
Starting from {x_{0},s_{0}}, recurrent use of Eq. (11) for the first few pairs {x_{k,sk}} leads to:
A summation series for I(x_{0},s_{0}) is obtained by sequentially inserting the I(x_{kb},s_{kb}) and I(x_{ka},s_{ka}) into the corresponding I(x_{k},s_{k}). In doing so, each I(x_{k},s_{k}) turns into a double summation of increasingly higher dimension integrals. In the PBMC framework, each of those integrals is estimated by their integrands at properly selected values of the integration variables.
The overall process, however, is greatly simplified if at each step only one of the two summations is pursued. The structure of Eq. (11), with coefficients 1−a(x_{k},x_{kb}) and a(x_{k},x_{kb}), suggests the way to proceed. In the more general case, it is convenient to draw a random number ϱ∈[0, 1] and follow the summation if a(x_{k},x_{kb}) ≥ ϱ> 0 and the ℬ summation if a(x_{k},x_{kb}) <ϱ< 1.
The ultimate goal of the PBMC algorithm is to estimate the Stokes vector at the detector from a number n_{ph} of single photon experiments. For a fixed {x_{0},s_{0}}, this is done by evaluating (A.1)where each ⟨ I^{iph}(x_{0},s_{0}) ⟩ is an estimate based on a single photon simulation. The estimate becomes statistically meaningful by repeating the process n_{ph} of times. When the integration is over the planetary disk, Eq. (A.1) is replaced by Eqs. (29) or (31), and the position x_{0} of entry of the simulated photon into the atmosphere is determined with the corresponding sampling scheme.
The process that yields ⟨ I(x_{0},s_{0}) ⟩ (index i_{ph} omitted) starts by tracing the ray from x_{0} in the −s_{0} direction, following the instructions below:

1.
Initialise ⟨ I(x_{0},s_{0}) ⟩ = 0, x_{k} = x_{0}, s_{k} = s_{0}, ℍ_{k} = ℍ_{0} (≡unity matrix) and w_{k} = 1.

2.
Determine x_{kb} and a(x_{k},x_{kb}). Then,

(a)
If r_{g}(x_{kb}) = 0 or if vector −s_{k} does not intersect the planet’s surface:

g = a(x_{k},x_{kb}).

Go for A at step 3.


(b)
Otherwise,

g = 1.

Draw a random number ϱ ∈ [0, 1]. Then,

If a(x_{k},x_{kb}) ≥ ϱ> 0, go for A at step 3.

If a(x_{k},x_{kb}) <ϱ< 1, go for B at step 4.



(a)

3.
Going for A: Collision in between boundaries.

Draw a random number ϵ_{ka} ∈ [0, 1] and displace the photon from x_{k} to x_{ka} along −s_{k} according to Eq. (7).

At x_{ka}, draw a random number ζ_{ka} ∈ [0, 1] and find θ_{ka} from the probability distribution function of Eq. (23). This is done by tabulation and subsequent inversion of ∫f_{θ}(θ)dθ∈ [0, 1].

At x_{ka}, find φ_{ka} from the probability distribution function of Eq. (24). This is done by means of the rejection method and the fact that by construction 2πf_{φ  θ}(φ  θ) ≤ 2.

Update:

⟨ I(x_{0},s_{0}) ⟩ ← ⟨ I(x_{0},s_{0}) ⟩ + w_{k} g ϖ(x_{ka})t(x_{ka},x_{⊙})ℍ_{k} ℙ(x_{ka},s_{k},s_{⊙})F_{⊙}

w_{k} ← w_{k} g ϖ(x_{ka})

ℍ_{k} ← ℍ_{k}ℙ(x_{ka},s_{k},s_{ka})/(ℍ_{k}ℙ(x_{ka},s_{k},s_{ka}))_{1,1}.

x_{k} ← x_{ka} and s_{k} ← s_{ka}.


If w_{k} ≥ ε_{ph}, go to step 2. Otherwise, go to step 5.


4.
Going for B: Collision at the bottom boundary. Draw random numbers ζ_{kb},η_{kb}∈[0, 1].

Displace the photon from x_{k} to x_{kb} along −s_{k}.

At x_{kb}, evaluate φ_{kb} = ζ_{kb}*2π and cos( .

Update:

⟨ I(x_{0},s_{0}) ⟩ ← ⟨ I(x_{0},s_{0}) ⟩ + w _{k} g r_{g}(x_{kb}) (n(x_{kb})·s_{⊙}) /π t(x_{kb},x_{⊙})ℍ_{k}F_{⊙}

w_{k} ← w_{k} g r_{g}(x_{kb})

ℍ_{k} ← ℍ_{k}ℙ(x_{kb},s_{k},s_{kb})/(ℍ_{k}ℙ(x_{kb},s_{k},s_{kb}))_{1,1}.

x_{k} ← x_{kb} and s_{k} ← s_{kb}.


If w_{k} ≥ ε_{ph}, go to step 2. Otherwise, go to step 5.


End of ⟨ I(x_{0},s_{0}) ⟩ loop.
The loop ends when the weight w_{k} reaches a userdefined threshold ε_{ph} that truncates the summation series. The value of ε_{ph} has an impact on both the solution’s accuracy and execution time. Values in the range 10^{4}–10^{5} are adequate for required accuracies of about 0.1% in the I Stokes element.
BMC algorithms with classical sampling schemes for photon propagation directions have a structure similar to the above. In the classical sampling scheme, the incident photon directions s′ are sampled from the local f(θ,φ) for q ≡ 0 (Fig. 3, top panel) rather than from the full f(θ,φ) of Eqs. (23)–(24) (Fig. 3, panels for  q  > 0). The classical BMC algorithm can be seen as a variation to the above algorithm, with the main difference being the definition of ℍ_{k}.
Appendix B: Comparison with de Haan et al. (1987)
Solutions to the Stokes vector in a conservative, hazeL atmosphere of optical thickness equal to 1 and cos(SPA) = 0.1.
Appendix C: Comparison with DISORT
For an assessment of the PBMC algorithm against the problem of radiation emerging from a conservative, nonpolarising Rayleigh atmosphere above a Lambert reflecting surface, we built a battery of solutions with DISORT (Stamnes et al. 1988). DISORT is a welldocumented and thoroughly tested solver of the scalar RTE for monochromatic radiation in multiplescattering media based on the discreteordinate method.
Table 1 summarises the model parameters and their ranges for the comparison exercise. They include the atmospheric optical thickness, surface albedo, and the three angles of Fig. 5. The atmospheric single scattering albedo is taken to be one as corresponds to conservative scattering. In total, the battery comprises 15 876 test cases that explore both optically thin and thick atmospheres with viewing/illumination angles from zenith inclination to nearly horizontal pointing. The PBMC calculations were carried out with n_{ph}=10^{4}, 10^{5}, 10^{6}, and 10^{7} photon realisations. Figure C.1 shows the relative differences, defined as δI=(I_{PBMC}−I_{ref}) /I_{ref}×100, between the computations with DISORT (reference model) and our PBMC algorithm. Median values for  δI  are listed in Table C.1. The convergence rate is consistent with the expected law for MC integration.
Fig. C.1 Relative differences in intensity, δI, between DISORT and our PBMC algorithm for conservative, nonpolarising, Rayleigh atmospheres. The full set of cases is summarised in Table 1 of the main text. Dashed vertical lines separate cases run with different atmospheric optical thickness (τ, in the graph). The PBMC algorithm simultaneously runs 7 × 7 = 49 configurations for the illumination geometry. 
Appendix D: Comparison with Garcia & Siewert (1986)
Scattering in media inspired by the Venus atmosphere has been investigated by Garcia & Siewert (1986) with a generalised spherical harmonics method. We refer to two of their study cases, namely: L = 13 (r_{eff} = 0.2 μm, v_{eff} = 0.07, λ = 0.951 μm,
refractive index of 1.44; akin to the mode1 haze atop the upper clouds of Venus) and L = 60 (r_{eff} = 1.05 μm, v_{eff} = 0.07, λ = 0.782 μm, refractive index of 1.43; akin to the mode2 droplets that make up the upper clouds of Venus). Tables D.1–D.16 show the Garcia & Siewert (1986) solutions and our PBMC calculations for n_{ph} up to 10^{9}.
Stokes vector element I at the top of the atmosphere for an L = 13 atmosphere: r_{eff} = 0.2 μm, v_{eff} = 0.07, λ = 0.951 μm, refractive index of 1.44, with optical thickness of one, single scattering albedo ϖ = 0.99, and surface albedo of 0.1.
Stokes vector element I at the top of the atmosphere for an L = 60 atmosphere: r_{eff} = 1.05 μm, v_{eff} = 0.07, λ = 0.782 μm, refractive index of 1.43, with optical thickness of one, single scattering albedo ϖ = 0.99, and surface albedo of 0.1.
Same as Table D.9 for Stokes vector element Q.
Same as Table D.9 for relative azimuth between the incident and emerging directions equal to π/2.
Same as Table D.11 for Stokes vector element Q.
Same as Table D.11 for Stokes vector element U.
Same as Table D.11 for Stokes vector element V.
Same as Table D.9 for relative azimuth between the incident and emerging directions equal to π.
Same as Table D.15 for Stokes vector element Q.
Appendix E: Rayleigh phase curves
Fig. E.1 Diskintegrated phase curves for F_{I} and F_{Q}/F_{I} for a Rayleighscattering atmosphere with ϖ = 1 and r_{g} = 1 and optical thicknesses from 0.1 to 10. Our PBMC calculations for n_{ph} = 10^{4} (red curves) are compared to the calculations by Buenzli & Schmid (2009) using a forward Monte Carlo algorithm (black curves) and at least 2 × 10^{6} photons per exit direction bin over the 30–120°αrange. Both sets of curves are nearly indistinguishable on the scale of the graph, even for the relatively small number of photon realisations used in our PBMC calculations. 
All Tables
Parameters in the investigation of conservative Rayleigh scattering in planeparallel atmospheres.
Solutions to the Stokes vector in a conservative, hazeL atmosphere of optical thickness equal to 1 and cos(SPA) = 0.5.
Parameters in the investigation of disk integration for both conservative and nonconservative Rayleigh atmospheres.
Solutions to the Stokes vector in a conservative, hazeL atmosphere of optical thickness equal to 1 and cos(SPA) = 0.1.
Stokes vector element I at the top of the atmosphere for an L = 13 atmosphere: r_{eff} = 0.2 μm, v_{eff} = 0.07, λ = 0.951 μm, refractive index of 1.44, with optical thickness of one, single scattering albedo ϖ = 0.99, and surface albedo of 0.1.
Stokes vector element I at the top of the atmosphere for an L = 60 atmosphere: r_{eff} = 1.05 μm, v_{eff} = 0.07, λ = 0.782 μm, refractive index of 1.43, with optical thickness of one, single scattering albedo ϖ = 0.99, and surface albedo of 0.1.
Same as Table D.9 for relative azimuth between the incident and emerging directions equal to π/2.
Same as Table D.9 for relative azimuth between the incident and emerging directions equal to π.
All Figures
Fig. 1 Top: definition of the incident, s′, and emergent, s, photon directions at a scattering event within the atmosphere. The xyz axes form a rest reference frame fixed to the planet. The differential solid angle dΩ(s′) = sinθdθdφ is defined with s serving as polar axis. Angles θ ∈ [ 0,π] and φ ∈[0, 2π]. In the backtracing of photons of BMC algorithms, s is known at each collision and s′ must be sampled from the relevant scattering phase function. Angles i′ and i, both ∈[0, π], are needed for consistent referencing of the Stokes vector throughout the scattering process. Vectors {, and } and {e_{1}, e_{2} and e_{3}} define righthanded coordinate systems at the meridional planes of the incident and emergent photon directions, respectively. Bottom: definition of the incident, s′, and emergent, s, photon directions at a reflection event at the local surface (plane ). Here, n is the inwardpointing normal vector at the surface, and ẑ is oriented along −n. The differential solid angle dΩ(s′) = sinθdθdφ is defined with ẑ serving as polar axis. 

In the text 
Fig. 2 In black, sketch demonstrating the construction of the {x_{k},s_{k}} pairs starting from {x_{0},s_{0}} as the photon is traced back from the detector through the medium. Vectors are pointed in the direction of photon propagation, which is the reverse of the direction of integration in the BMC algorithm. In principle, a {x_{k},s_{k}} pair can lead to two new {x_{ka},s_{ka}} and {x_{kb},s_{kb}} pairs. In the MC implementation of the algorithm, a scheme based on the coefficients of the ℒ_{ℬ}+ℬ and + operators, Eq. (11), determines whether the photon’s next move occurs within the atmospheric medium or whether the photon moves onto the planet’s surface. In blue, one specific photon trajectory within the family of possible trajectories. For this specific trajectory, the red arrows denote the direction of the unscattered stellar photons. 

In the text 
Fig. 3 Probability density function f(θ,φ)=f_{θ}(θ)f_{φ  θ}(φ  θ) in the precondioned sampling scheme of photon propagation directions, Eqs. (23)–(24), for a Rayleighscattering medium. Note the changes in f(θ,φ) with q, especially near the maximum of  b_{1}(θ) /a_{1}(θ)  for θ = 90°. By ignoring polarisation, the classical sampling scheme determines the θ and φ values of the incident propagation direction from f(θ,φ;q ≡ 0). It is apparent that the classical sampling scheme is more likely to fail in strongly polarising media that involve high q values during the backtracing of photons. 

In the text 
Fig. 4 Geometrical parameters relevant to the integration over the planet’s “visible” disk. For a pair of u_{i} and v_{i} values, N is the location on the disk where the observer’s line of sight intercepts the atmosphere. N is equivalent to x_{0} in our implementation of the PBMC algorithm. 

In the text 
Fig. 5 Illumination and viewing angles for the planeparallel atmosphere test cases discussed in Sect. 3. 

In the text 
Fig. 6 Differences in intensity, δI, for the solution of the VRTE in conservative, Rayleighscattering atmospheres. The algorithm uses the preconditioned (top) and classical (bottom) sampling schemes for photon propagation directions. 

In the text 
Fig. 7 Differences in polarisation, δP, for the solution of the VRTE in conservative, Rayleighscattering atmospheres. The algorithm uses the preconditioned sampling scheme for photon propagation directions. 

In the text 
Fig. 8 Convergence history of the I Stokes element for a conservative Rayleigh atmosphere over a black surface with cos(SPA) = cos(OPA) = 1 and two different optical thicknesses. Thin and thick curves represent the solutions obtained with the preconditioned and classical sampling schemes, respectively. Abrupt changes in the solution with the classical sampling scheme for optical thickness of 4 occur, but they are not discernible on the scale of the graph. 

In the text 
Fig. 9 Polarisation in single scattering for monodisperse droplets of various radii and real refractive index equal to 1.53 at λ = 0.63 μm. 

In the text 
Fig. 10 Convergence history for scattering by the monodisperse droplets of Fig. 9. Thin and thick curves represent the solutions obtained with the preconditioned and classical sampling schemes, respectively. The classical scheme produces inconsistent solutions for strongly polarising conditions. The calculations assumed optical thickness equal to 16 and cos(OPA) = cos(SPA) = 1. 

In the text 
Fig. 11 Computational time per phase curve point for n_{ph}=10^{4} and Rayleighscattering calculations. For a full phase curve with, for instance, 34 evenly separated points from 7.5 to 172.5° (as for Fig. E.1), the computational time is 34 times what is indicated in the plot. The curves correspond to different values of the single scattering albedo of the medium, ϖ, and the Lambert surface albedo, r_{g}. 

In the text 
Fig. 12 Polarisation phase curves for Venus. Black diamonds are the measurements used in Hansen & Hovenier (1974). Colour symbols and curves are our PBMC calculations for n_{ph}=10^{4} and 10^{5}, respectively. This figure can be compared to Figs. 4, 8, 9, and 11, 12 in Hansen & Hovenier (1974). 

In the text 
Fig. 13 Polarisation phase curves for Venus calculated with our PBMC algorithm with =10^{5} at nearinfrared wavelengths for a H_{2}SO_{4}/H_{2}O dilution at 75%. We adopted a single scattering albedo of 1 in all cases. n_{r} is the refractive index in each case. 

In the text 
Fig. 14 Phase curves F_{I} in the optical for Venuslike planets with atmospheres stratified according to the given scale heights. 

In the text 
Fig. C.1 Relative differences in intensity, δI, between DISORT and our PBMC algorithm for conservative, nonpolarising, Rayleigh atmospheres. The full set of cases is summarised in Table 1 of the main text. Dashed vertical lines separate cases run with different atmospheric optical thickness (τ, in the graph). The PBMC algorithm simultaneously runs 7 × 7 = 49 configurations for the illumination geometry. 

In the text 
Fig. E.1 Diskintegrated phase curves for F_{I} and F_{Q}/F_{I} for a Rayleighscattering atmosphere with ϖ = 1 and r_{g} = 1 and optical thicknesses from 0.1 to 10. Our PBMC calculations for n_{ph} = 10^{4} (red curves) are compared to the calculations by Buenzli & Schmid (2009) using a forward Monte Carlo algorithm (black curves) and at least 2 × 10^{6} photons per exit direction bin over the 30–120°αrange. Both sets of curves are nearly indistinguishable on the scale of the graph, even for the relatively small number of photon realisations used in our PBMC calculations. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.