© 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 stand-alone technique and in combination with photometry. In the solar system, polarimetric observations made from space-borne and ground-based 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 ground-based observatories are equipped with polarimeters for either spectroscopy or imaging. Ground-based observations of the outer planets, however, only have partial coverage of the Sun-target-Earth 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 European-Extremely 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 pre-conditioned BMC (PBMC) algorithm and show that pre-conditioning the scattering matrix with information from prior collisions (in the order of backward integration) eliminates the numerical difficulties. Pre-conditioning 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 disk-resolved and disk-integrated 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 disk-integrated 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 spherical-shell planet. In Sect. 3, we assess the performance of the classical and pre-conditioned 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 follow-up 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 next-event point-estimator 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. Limb-viewing 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 user-defined 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₁ = 3(1 + cos²θ) / 4, b₁ = 3(−1 + cos²θ) / 4, a₃ = 3cosθ/ 2, and b₂ = 0.

Equation (1) admits the formal solution for the Stokes vector at {x_k,s_k}: (5)On the right-hand 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 arc-length 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 inwards-pointing normal vector at the surface, is the four-by-four depolarizing matrix with _1,1 = 1 as the only non-zero 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₀,s₀}, 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₀,s₀) 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₀,s₀) 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 so-called 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.

2.2. Monte Carlo integration

The essence of MC integration is to estimate multi-dimensional 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₀,s₀) 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₀,s₀) obtained from recurrent use of Eq. (11) contains multi-dimensional 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 pre-conditioned 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₀,s_0a) = (≡unity matrix), ℍ(s_0a,s_0aa) = ℙ(s₀,s_0a), ℍ(s_0aa,s_0aaa) = ℙ(s₀,s_0a)ℙ(s_0a,s_0aa), ..., effectively pre-condition 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 pre-conditioning 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 pre-conditioned sampling scheme for photon propagation directions. This scheme is at the core of our PBMC algorithm.

Expanding Eq. (19) yields insight into the pre-conditioned 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₁(θ) ≥ 0 and | b₁ | ≤ | a₁ | (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₂) = (−b₂), 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 pre-conditioned 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₁(θ) /a₁(θ), q and u may take absolute values close to one through the photon simulations. The consequences of this are investigated below.

2.4. Disk-integration schemes

We are interested in the radiation emerging from both disk-resolved and disk-integrated planetary atmospheres. We derive two disk-integration 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 disk-integrated radiation scattered from a planet over its “visible” disk. In this context, “visible” refers to the disk portion that appears illuminated by single-scattered 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 observer-to-planet distance, respectively. We normalise F by eliminating the (ρ/ Δ)² factor from Eq. (25).

Rather than working with longitudes, ζ_d, and co-latitudes, η_d, it is convenient to introduce the two auxiliary variables: such that, after some manipulations, Eq. (25) transforms into (28)The pre-multiplying 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₀ in the implementation of the algorithm of Appendix A. For a sufficiently remote observer, s₀ 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 pre-calculated 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₀, and s₀ 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/ρ)², integration over the projected surface element rdrdΘ leads to (30)which, after eliminating the (ρ/ Δ)² 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 effectively-scattering 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 disk-integration schemes to both Rayleigh and Venus-like 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 plane-parallel 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.

3.1. Non-polarised 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.

3.2. Polarised Rayleigh scattering

Coulson et al. (1960) tabulated solutions for the elements of the Stokes vector in conservative, polarising, Rayleigh-scattering 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⁷ with our PBMC algorithm in its VRTE mode. For comparison, we utilised both the classical and pre-conditioned sampling schemes introduced in Sect. 2.3.

Figure 6 shows δI (=(I_BMC−I_ref) /I_ref × 100) for the pre-conditioned (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 pre-conditioned 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.

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 multi-dimensional integral of Eq. (17) can be approximated by separate integrals as given by Eq. (18) plus a subsequent correction, becomes inappropriate for specific configurations.

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₁(θ)/a₁(θ) 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₁(θ)/a₁(θ)| 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 pre-conditioned 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.

3.3. Polarised Mie scattering

We tested our PBMC algorithm against a number of VRTE solutions in Mie-scattering media. The Stokes vectors for radiation emerging from a conservative atmosphere with so-called haze-L scattering particles have been tabulated by de Haan et al. (1987) from calculations based on the doubling-adding 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⁹. Typically, solutions accurate to within one per cent in I are obtained for n_ph = 10⁵. 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 Venus-like 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 disk-integrated 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 disk-integrated 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 disk-integration schemes of Sect. 2.4, we utilised a few configurations relevant to both Rayleigh and Venus-like atmospheres. Essentially, the disk-integration scheme selects the entry point of the photon into the atmosphere. The three-dimensional 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 star-planet-observer phase angles.

4.1. Rayleigh phase curves

Buenzli & Schmid (2009) have investigated with an FMC algorithm the phase curves of Rayleigh-scattering 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 plane-parallel 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 non-Rayleigh forms of scattering are needed. Thus, Stam et al. (2006) have devised an efficient technique for disk integration based on the plane-parallel approximation that can deal with arbitrary scattering particles for planets with horizontally uniform atmospheres.

We produced Rayleigh phase curves for the configurations listed in Table 3 with the visible-disk 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⁶, 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 visible-disk 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⁴. 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⁴. Computational times depend on the platform, and are not often published in the literature, which prevents a comparison with the performance of other algorithms.

4.2. Venus phase curves

Venus is a well-known example that demonstrates the potential of disk-integrated polarimetry in the remote investigation of clouds (Coffeen 1969; Hansen & Hovenier 1974). The disk-integrated 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 near-infrared, 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.

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⁴, and the underlying solid curves are the calculations for n_ph = 10⁵. For the modelling, we use the prescriptions for particle size distributions (gamma-distribution, 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 Rayleigh-scattering component, f_R (see Hansen & Hovenier 1974), which becomes important at UV wavelengths, and various r_eff values. For n_ph = 10⁴ and the visible-disk 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⁴ 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⁴. 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 real-only refractive indices based on a 75% H₂SO₄/H₂O solution by mass (Hansen & Hovenier 1974). They are shown in Fig. 13, which further illustrates the sensitivity of polarisation to wavelength.

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 e-folding length for changes in the γ extinction coefficient in the vertical) of 4, 8, 16, and 32 km. Both disk-integration 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 disk-integration 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.

5. Summary and future work

We have presented a novel pre-conditioned 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 pre-conditions the scattering matrix before sampling the incident propagation direction at a photon collision. Pre-conditioning 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 pre-conditioned 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 disk-integration 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 disk-integration scheme to horizontally inhomogeneous planets, thus accounting for a full three-dimensional description of the planet. Similar ideas will also be explored for disk-integrated thermal emission and for the simultaneous spectral-and-disk integration of both scattered and thermally-emitted radiation.

A one-slab, plane-parallel 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

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₀,s₀}, recurrent use of Eq. (11) for the first few pairs {x_k,sk} leads to:

A summation series for I(x₀,s₀) 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₀,s₀}, this is done by evaluating (A.1)where each ⟨ I^i_ph(x₀,s₀) ⟩ 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₀ of entry of the simulated photon into the atmosphere is determined with the corresponding sampling scheme.

The process that yields ⟨ I(x₀,s₀) ⟩ (index i_ph omitted) starts by tracing the ray from x₀ in the −s₀ direction, following the instructions below:

• 1.

  Initialise ⟨ I(x₀,s₀) ⟩ = 0, x_k = x₀, s_k = s₀, ℍ_k = ℍ₀ (≡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.

• 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₀,s₀) ⟩ ← ⟨ I(x₀,s₀) ⟩ + 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₀,s₀) ⟩ ← ⟨ I(x₀,s₀) ⟩ + w _k g r_g(x_kb) (n(x_kb)·s_⊙) /π t(x_kb,x_⊙)ℍ_kF_⊙

    • 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₀,s₀) ⟩ loop.

The loop ends when the weight w_k reaches a user-defined 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)

Appendix C: Comparison with DISORT

For an assessment of the PBMC algorithm against the problem of radiation emerging from a conservative, non-polarising Rayleigh atmosphere above a Lambert reflecting surface, we built a battery of solutions with DISORT (Stamnes et al. 1988). DISORT is a well-documented and thoroughly tested solver of the scalar RTE for monochromatic radiation in multiple-scattering media based on the discrete-ordinate 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⁴, 10⁵, 10⁶, and 10⁷ 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.

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 mode-1 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 mode-2 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⁹.

Appendix E: Rayleigh phase curves

All Tables

All Figures