A Theory of Exoplanet Transits with Light Scattering

Transit spectroscopy (Seager & Sasselov 2000; Brown 2001; Hubbard et al. 2001) is currently the leading technique for studying exoplanet atmospheric composition. Following the discovery of the first exoplanet atmosphere (Charbonneau et al. 2002), transit observations have enabled the characterization of a number of different exoplanets for atmospheric molecular species, clouds, and/or hazes (Pont et al. 2008; Swain et al. 2008; Sing et al. 2009; Bean et al. 2010; Fraine et al. 2014; Knutson et al. 2014b). Recent observational results have shown that hot Jupiter transit spectra demonstrate a complex continuum from clear sky to heavily clouded conditions (Sing et al. 2016; Stevenson 2016), and that cloudiness remains a key factor into the regime of lower mass planets (Line et al. 2013; Knutson et al. 2014a; Kreidberg et al. 2014).

The field of exoplanet transit spectroscopy will be revolutionized with the anticipated launch of NASA's James Webb Space Telescope (JWST) in 2018 (Gardner et al. 2006). Over the course of the five year mission for JWST, the observatory is expected to provide in-depth observations of many tens of transiting exoplanets (Beichman et al. 2014). Some of these observations will probe planets in the poorly understood 2–4 Earth mass regime (Deming et al. 2009; Batalha et al. 2015). Excitingly, JWST may even be capable of characterizing temperate Earth-sized planets (Kaltenegger & Traub 2009; Cowan et al. 2015; Barstow et al. 2016), though the ability to conduct such studies will depend largely on how well systematic noise sources can be constrained (Greene et al. 2016).

As the quality of transit spectrum observations continues to improve, so should models of exoplanet transits. Thus, certain processes initially thought to be of second-order importance should be revisited and possibly added to modeling tools. For example, atmospheric refraction, which was initially shown to be unimportant for the case of certain hot Jupiters (Hubbard et al. 2001), has recently been shown to be critically important for understanding some terrestrial exoplanet spectra (Muñoz et al. 2012; Bétrémieux & Kaltenegger 2014; Misra et al. 2014; Bétrémieux & Kaltenegger 2015) and, possibly, gas giant transit spectra (Dalba et al. 2015; Bétrémieux 2016). Additionally, refraction can affect the shape of an exoplanet transit light curve (Hui & Seager 2002; Sidis & Sari 2010).

Beyond refraction, another process that has seen little study with regards to exoplanet transits is light multiple scattering. Hubbard et al. (2001) used a plane-parallel Monte Carlo scattering model to determine the significance of a molecular Rayleigh scattering "glow" around the limb of a transiting exoplanet. For the case of HD 209458b, these authors found that the contribution of Rayleigh scattering glow to the transit light curve is negligible. Similar conclusions were reached by Brown (2001), who used analytic arguments to deduce that multiple isotropic scatterings could have only a small impact on transit depths.

Since the study of Hubbard et al. (2001), scattering opacity in exoplanet transits has largely been treated as being equivalent to absorption opacity (e.g., Irwin et al. 2008; Line et al. 2013; Waldmann et al. 2015). This equivalence cannot always hold. For example, a very strong forward and conservatively scattering aerosol would only weakly attenuate a beam following a straight-line path from a stellar disk, passing through an exoplanet atmosphere, and traveling toward a distant observer. This issue was plainly demonstrated by de Kok & Stam (2012), who used a Monte Carlo scattering model to study the transmission of a pencil beam through a cloud or haze layer with variable optical thickness and whose scatterers had asymmetry parameters of either 0, 0.9, or 0.98 (their Figure 3).

In the work that follows, we do not seek to outline specific definitions for "hazes" or "clouds." In general, though, we use the former to refer to aerosols that form aloft and grow due to coagulation during the sedimentation process, thereby developing a distribution that may look something like Titan's tholin haze, which has a number density profile described (roughly) by an exponential with a constant scale height in Titan's upper atmosphere (Tomasko et al. 2008). For clouds, we envision an aerosol distribution formed by lifting a gas to its condensation point, thereby developing a distinct cloud base with an overlying cloud deck whose thickness is controlled by mixing processes (see, e.g., Ackerman & Marley 2001). Advances in exoplanet and brown dwarf cloud models have been reviewed recently by Marley et al. (2013, pp. 367–391). The most popular tools include the Eddysed model (Ackerman & Marley 2001), where an upward diffusion of condensible vapor balances a downward sedimentation of aerosol particles, which has been successfully used to study both cloud and haze processes in transiting exoplanet atmospheres (Morley et al. 2013). Additionally, the sophisticated microphysical cloud treatments of Helling et al. (2001) continue to see development and application to transiting exoplanets (Lee et al. 2016).

Thus, the question of the relative importance of scattering, refraction, and absorption remains largely unexplored for a wide portion of exoplanet parameter space. Such studies have been hindered by the lack of available simulation tools, and also by a lack of a single, coherent theory of exoplanet transit spectroscopy. Significant advances can be made by developing such a theory, and by building models that implement these ideas.

In this paper, we construct a theory of exoplanet transit spectroscopy that spans key fundamental integral relations through to its efficient, vectorized implementation. Simplifications of the theory cover the geometric (i.e., straight-line) limit (that is currently used in computationally demanding spectral retrieval models; e.g., Benneke & Seager 2012; Line et al. 2012; Barstow et al. 2013; Lee et al. 2014; Waldmann et al. 2015; Morley et al. 2017) and cases that add refraction. In general, though, we emphasize the "full physics" scenario, where a three-dimensional Monte Carlo model is used to incorporate scattering effects.

Below, we introduce the concept of the "path distribution," which effectively separates the paths of photons (or rays) through a transiting exoplanet atmosphere from gaseous and/or aerosol absorption effects. As the former tends to vary weakly in wavelength, while the latter can vary strongly with wavelength, our approach is computationally efficient. Following the presentation of these ideas, we validate our implementation of the theory against a number of trusted simulation tools. We then give examples of the path distribution for a variety of different model atmospheres. Finally, we use our tools to present the first simulated transit spectra of hot Jupiters with multiply-scattering clouds.

Transit spectra typically probe low-density regions of exoplanetary atmospheres. Here, molecular absorption lines are relatively narrow, as the effects of pressure broadening are not dominant over Doppler broadening (as is the case in the deep atmosphere). Thus, the gas opacity can vary by orders of magnitude over narrow spectral ranges. By contrast, gaseous refractive indexes as well as gas and aerosol scattering properties tend to vary smoothly (and, sometimes, weakly) in wavelength. Since refraction and scattering, which are more computationally intensive processes to simulate, influence the path of a photon (or ray) through an atmosphere, significant advances can be made by outlining a technique that separates the processes which influence photon trajectories from those which influence absorption.

We introduce the concept of the ray (or photon) atmospheric path distribution, ${{ \mathcal P }}_{b}\,(h)$, where b is the impact parameter of the ray and h is altitude in the planetary atmosphere. We define the path distribution such that ${{ \mathcal P }}_{b}\,(h){dh}$ is the linear distance traversed by the ray at altitudes between h and h+dh. While ${{ \mathcal P }}_{b}$ is dimensionless, it can be thought of as having units of km linear distance per km vertical distance. The geometry of these parameters, for a simple case, is shown in Figure 1. The Appendix contains additional discussion on the use of the path distribution. Also, for convenience, Table 1 contains a summary of all symbols and variables used throughout this manuscript.

Table 1.  Symbol Usage

Symbol	Description
${\boldsymbol{A}}$	vector/array of area elements
a	planet-star orbital distance
${a}_{\lambda }$	absorptivity along a ray path
${\alpha }_{\lambda }$	atmospheric extinction coefficient
b	impact parameter
d	projected separation between planet center and star center
g	scattering asymmetry parameter, acceleration due to gravity
$H={R}_{{\rm{g}}}T/{mg}$	pressure scale height
h, ${h}_{{\rm{t}}}$	altitude, altitude at effective top of atmosphere
${I}_{{\rm{s}},\lambda }$	stellar surface brightness
${I}_{{\rm{p}},\lambda }$	surface brightness on planetary disk
${I}_{0,\lambda }$	disk-averaged stellar surface brightness
$\tilde{I}$	background to planetary disk surface brightness mapping
m	atmospheric mean molecular weight
$\mu =\cos (\theta )$	cosine of scattering angle
${\mu }_{{\rm{s}}}$	cosine of angle of incidence on stellar disk
$\hat{{\boldsymbol{\mu }}}$	direction of photon/ray travel
Nr	number of impact parameters in model grid
${N}_{{\rm{p}}}$	number of photons in a Monte Carlo simulation
${N}_{\mathrm{lay}}$	number of layers in atmospheric model
${n}_{\mathrm{ref}}$	atmospheric refractive index
$P(\mu )$	scattering phase function
${{ \mathcal P }}_{b}$, ${\boldsymbol{ \mathcal Q }}$	atmospheric path distribution
p	pressure
ϕ	polar angle when ray tracing, scattering azimuth angle
${R}_{{\rm{E}}}$	Earth radius
${R}_{{\rm{g}}}$	universal gas constant
${R}_{{\rm{J}}}$	Jupiter radius
${R}_{{\rm{p}}}$	fiducial planetary radius
${R}_{{\rm{p}},\lambda }$	wavelength-dependent planetary radius
${R}_{{\rm{s}}}$	stellar radius
r	radial distance from planet center
r0, ${\theta }_{0}$	integration bounds for stellar brightness/area
${r}_{{\rm{c}}}$	ray curvature due to refraction
sb	ray path
${t}_{\lambda }$	transmission along ray path
${\bar{{\boldsymbol{t}}}}_{\lambda }$	Monte Carlo average transmission vector
$\partial {\bar{{\boldsymbol{t}}}}_{\lambda }/\partial {\boldsymbol{\Delta }}{{\boldsymbol{\tau }}}_{\lambda ,a}$	Monte Carlo average transmission Jacobian matrix
${\rm{\Delta }}\tau$	layer vertical optical depth
τ	optical depth along ray path or slant optical depth
${\theta }_{r}$	local zenith angle for ray propagation direction
θ	scattering deflection angle, polar angle in disk integration
Ω, $d{\rm{\Omega }}$	solid angle, differential solid angle
ω	refraction integral/angle
${\tilde{\omega }}_{0}$	single-scattering albedo
ξ	random number between zero and one

Download table as:  ASCIITypeset image

Zoom In Zoom Out Reset image size

In transit spectroscopy, the essential quantity is the attenuation of a ray along its path, as this distinguishes the opaque portions of the atmosphere (that block light from the stellar disk) from the transparent portions. In 1D atmospheric structure models, which are the variety most commonly used for exoplanets, the atmospheric opacity is only a function of altitude (or pressure) and wavelength. Here, then, the optical depth along the path of a ray can be computed using the path distribution via

$$\begin{eqnarray}&&{\tau }_{\lambda ,b}={\int }_{0}^{\infty }{\alpha }_{\lambda }(h){{ \mathcal P }}_{b}\,(h){dh},\end{eqnarray} \tag{ 1 }$$

where ${\alpha }_{\lambda }$ is the atmospheric extinction (in units of inverse distance), and the integral can be taken to infinity as the extinction is zero at very large altitudes. Since the layer vertical differential optical depth is defined as $d{\tau }_{\lambda }(h)={\alpha }_{\lambda }(h){dh}$, this expression demonstrates that the path distribution can also be thought of as the linear optical depth encountered between ${\tau }_{\lambda }$ and ${\tau }_{\lambda }+d{\tau }_{\lambda }$. Similar extensions exist for quantities like the column mass and number densities. Note that Equation (1) is critical, as the strongly wavelength-dependent extinction is separated from the ray path. Using the standard definition of transmission,

$$\begin{eqnarray}&&{t}_{\lambda }(b)={e}^{-{\tau }_{\lambda ,b}},\end{eqnarray} \tag{ 2 }$$

a simple transit spectrum in the pure absorbing limit (i.e., where all optical depth is taken as absorption optical depth) can then be computed by considering the light transmitted through concentric annuli on the planetary disk, each at their own impact parameter

$$\begin{eqnarray}&&{\left(\displaystyle \frac{{R}_{{\rm{p}},\lambda }}{{R}_{{\rm{s}}}}\right)}^{2}=\displaystyle \frac{2}{{R}_{{\rm{s}}}^{2}}{\int }_{0}^{\infty }[1-{t}_{\lambda }(b)]{bdb},\end{eqnarray} \tag{ 3 }$$

where ${R}_{{\rm{p}},\lambda }$ is the wavelength-dependent planetary radius, ${R}_{{\rm{s}}}$ is the stellar radius, and ${({R}_{{\rm{p}},\lambda }/{R}_{{\rm{s}}})}^{2}$ is the transit depth.

In practice, model planetary atmospheres are defined on an altitude (or pressure) grid, and Equation (3) is computed using a sum over a collection of impact parameters (and assuming that the planetary disk is opaque below some fiducial radius, ${R}_{{\rm{p}}}$, which is either at the surface or deep in the atmosphere). In this case, the path distribution is a matrix, ${{ \mathcal P }}_{i,j}$, where ${{ \mathcal P }}_{i,j}{\rm{\Delta }}{h}_{j}$ is the path traversed through atmospheric layer "j" (whose thickness is ${\rm{\Delta }}{h}_{j}$) by a ray whose impact parameter is b_i. At any given wavelength, the transmission is now a sum over ${N}_{\mathrm{lay}}$ atmospheric layers

$$\begin{eqnarray}&&{t}_{\lambda ,i}=\mathrm{EXP}\,\left(-\displaystyle \sum _{j=1}^{{N}_{\mathrm{lay}}}\,{\alpha }_{\lambda ,j}{{ \mathcal P }}_{i,j}{\rm{\Delta }}{h}_{j}\right)=\mathrm{EXP}\,\left(-\displaystyle \sum _{j=1}^{{N}_{\mathrm{lay}}}\,{\rm{\Delta }}{\tau }_{\lambda ,j}{{ \mathcal P }}_{i,j}\right),\end{eqnarray} \tag{ 4 }$$

where we have used the definition of the layer vertical differential optical depth, ${\rm{\Delta }}{\tau }_{\lambda ,j}={\alpha }_{\lambda ,j}{\rm{\Delta }}{h}_{j}$. A similar expression was used by Bétrémieux (2016, their Equations (5) and (7)), where these authors work in numerically computed column number densities and molecular opacities, as compared to our dimensionless quantities. Note that Equation (4) can be easily written in matrix notation as,

$$\begin{eqnarray}&&{{\boldsymbol{t}}}_{\lambda }=1-{{\boldsymbol{a}}}_{\lambda }=\mathrm{EXP}\,(-{\boldsymbol{\Delta }}{{\boldsymbol{\tau }}}_{\lambda }\cdot {\boldsymbol{ \mathcal Q }}),\end{eqnarray} \tag{ 5 }$$

where we have defined the absorptivity vector, ${{\boldsymbol{a}}}_{\lambda }$. The pure absorption transit spectrum derives from a sum over N_r impact parameters,

$$\begin{eqnarray}&&{\left(\displaystyle \frac{{R}_{{\rm{p}},\lambda }}{{R}_{{\rm{s}}}}\right)}^{2}=\displaystyle \frac{1}{{R}_{{\rm{s}}}^{2}}\left({R}_{{\rm{p}}}^{2}+2\displaystyle \sum _{i=1}^{{N}_{r}}[1-{t}_{\lambda ,i}]{b}_{i}{\rm{\Delta }}{b}_{i}\right),\end{eqnarray} \tag{ 6 }$$

where ${\rm{\Delta }}{b}_{i}$ is the thickness of the impact parameter gridpoint. Continuing with the matrix notation, if we define a vector of annulus areas as ${A}_{i}=2\pi {b}_{i}{\rm{\Delta }}{b}_{i}$, then the transit spectrum can be written simply as,

$$\begin{eqnarray}&&{\left(\displaystyle \frac{{R}_{{\rm{p}},\lambda }}{{R}_{{\rm{s}}}}\right)}^{2}=\displaystyle \frac{1}{{R}_{{\rm{s}}}^{2}}\left({R}_{{\rm{p}}}^{2}+\displaystyle \frac{1}{\pi }{{\boldsymbol{a}}}_{\lambda }\cdot {\boldsymbol{A}}\right).\end{eqnarray} \tag{ 7 }$$

We briefly note that the transit spectra expressions given in Equations (3), (6), and (7) are in the pure absorption limit, and assume that all rays trace back to the stellar disk, and that the star has uniform surface brightness. While these are common assumptions when computing transit spectra for model atmospheres, we will now discuss the more general multiple scattering case. The Appendix contains brief details about computing the path distribution in the geometric limit or in the case that refraction is considered.

2.1. Multiple Scattering Path Distributions and Transmissions

Cases that include clouds, especially strongly forward scattering clouds, require three-dimensional radiative transfer treatments to compute accurate transit spectra, and are well suited to Monte Carlo models. Here, the path distributions for a number of photons, ${N}_{{\rm{p}}}$, are used to derive an average transmission for a grid of impact parameters. The individual path distributions are determined using a Monte Carlo approach, and only consider the scattering optical depths (${\tau }_{\lambda ,{\rm{s}}}={\tilde{\omega }}_{0}{\tau }_{\lambda ,{\rm{e}}}$, where ${\tilde{\omega }}_{0}$ is the single-scattering albedo, and subscript "s" indicates "scattering," while "e" indicates "extinction") since only ${\tau }_{\lambda ,{\rm{s}}}$ affects the path of rays through the atmosphere. Given the path distribution, ${{\boldsymbol{ \mathcal Q }}}_{m}$, for photon "m," the transmission for this photon is,

$$\begin{eqnarray}&&{{\boldsymbol{t}}}_{\lambda ,m}=\mathrm{EXP}\,(-{\boldsymbol{\Delta }}{{\boldsymbol{\tau }}}_{\lambda ,a}\cdot {{\boldsymbol{ \mathcal Q }}}_{m}),\end{eqnarray} \tag{ 8 }$$

where ${\boldsymbol{\Delta }}{{\boldsymbol{\tau }}}_{\lambda ,a}$ is a vector of the wavelength-dependent layer differential absorption optical depths, and the transmission averaged over all photons is simply

$$\begin{eqnarray}&&{\bar{{\boldsymbol{t}}}}_{\lambda }=\displaystyle \frac{{\displaystyle \sum }_{m=1}^{{N}_{{\rm{p}}}}\,{{\boldsymbol{t}}}_{\lambda ,m}}{{N}_{{\rm{p}}}}.\end{eqnarray} \tag{ 9 }$$

The execution of our Monte Carlo model follows a standard approach. For a given impact parameter, b_i, a photon is directed into the atmosphere along the $-\hat{{\boldsymbol{x}}}$ direction in Figure 2, and then tracked through the atmosphere. As this photon originated from the direction of the observer, we are performing a so-called "backward" Monte Carlo simulation. The optical distance the photon travels, either initially or following a scattering event, is determined by randomly sampling the scattering optical depth distribution, with,

$$\begin{eqnarray}&&\xi ={\int }_{0}^{{\tau }_{{\rm{s}}}}f(\tau )d\tau ,\end{eqnarray} \tag{ 10 }$$

where ξ is a random number between 0 and 1, and $f(\tau )d\tau$ is the probability that a photon scatters between τ and $\tau +d\tau$. As $f(\tau )d\tau ={e}^{-\tau }d\tau$, we have,

$$\begin{eqnarray}&&{\tau }_{{\rm{s}}}=-\mathrm{ln}(1-\xi ).\end{eqnarray} \tag{ 11 }$$

For whichever layer (centered at h_j) the photon is currently in, the pathlength corresponding to ${\tau }_{{\rm{s}}}$ is,

$$\begin{eqnarray}&&s={\tau }_{{\rm{s}}}/{\alpha }_{{\rm{s}},j},\end{eqnarray} \tag{ 12 }$$

where ${\alpha }_{{\rm{s}},j}$ is the layer scattering extinction coefficient. Given a pathlength, a photon location (${\boldsymbol{r}}=x\hat{{\boldsymbol{x}}}+y\hat{{\boldsymbol{y}}}+z\hat{{\boldsymbol{z}}}$), and trajectory ($\hat{{\boldsymbol{\mu }}}={\mu }_{x}\hat{{\boldsymbol{x}}}+{\mu }_{y}\hat{{\boldsymbol{y}}}+{\mu }_{z}\hat{{\boldsymbol{z}}}$), the photon position is updated according to $x\to x+{\mu }_{x}s$, $y\to y+{\mu }_{y}s$, and $z\to z+{\mu }_{z}s$.

Zoom In Zoom Out Reset image size

Depending on the size of s, the photon either experiences a scattering event within layer j, in which case $s/{\rm{\Delta }}{h}_{j}$ is added to the path distribution for this photon (i.e., ${{\boldsymbol{ \mathcal Q }}}_{m}$), or the photon exits the layer before a scattering event occurs. For the latter, s must be compared to the distance to the layer boundaries. Given the trajectory of the photon, and its current radial position, $r=\sqrt{{x}^{2}+{y}^{2}+{z}^{2}}$, the altitude of the photon along its path is

$$\begin{eqnarray}h(s)+{R}_{{\rm{p}}} & = & r(s)\\ & = & \sqrt{{(x+{\mu }_{x}s)}^{2}+{(y+{\mu }_{y}s)}^{2}+{(z+{\mu }_{z}s)}^{2}}.\end{eqnarray} \tag{ 13 }$$

Thus, the quadratic equations (in s),

$$\begin{eqnarray}&&h({s}_{j})={h}_{j}\pm {\rm{\Delta }}{h}_{j}/2,\end{eqnarray} \tag{ 14 }$$

give the distances to the layer boundaries. When a photon travels to a layer boundary before a scattering event, ${s}_{j}/{\rm{\Delta }}{h}_{j}$ is added to ${{\boldsymbol{ \mathcal Q }}}_{m}$, the cumulative optical depth experienced by the photon is updated as $\tau \to \tau +{\alpha }_{{\rm{s}},j}{s}_{j}$, and the photon is passed to the appropriate layer (either $j-1$ or $j+1$). This process of randomly generating ${\tau }_{{\rm{s}}}$ and passing the photon through sequential layers is repeated until the photon either exits the atmosphere, or reaches the lower boundary of the model and is considered absorbed. (Note that, as with any transit spectrum model of gaseous worlds, the lower boundary must be set deep enough to not influence the simulated spectrum.) For photons that exit the atmosphere, the position and trajectory at exit determine whether or not the photon will intersect the stellar disk along a straight-line trajectory.

When a scattering event occurs, a scattering angle must be sampled from the scattering phase function and the photon trajectory must then be updated. Assuming that the scattering phase function, P, is only a function of the cosine of a single angle, $\mu =\cos (\theta )$, then a randomly sampled scattering deflection angle is determined via,

$$\begin{eqnarray}&&\xi ={\int }_{-1}^{\mu }P(\mu ^{\prime} )d\mu ^{\prime} ,\end{eqnarray} \tag{ 15 }$$

and an azimuthal scattering angle, ϕ, is sampled uniformly, with,

$$\begin{eqnarray}&&\phi =2\pi \xi .\end{eqnarray} \tag{ 16 }$$

As the scattering angles are referenced from the propagation direction, a transformation must be done to convert the initial propagation direction, $\hat{{\boldsymbol{\mu }}}$, and the scattering angles into a new propagation direction, $\hat{{\boldsymbol{\mu }}}^{\prime}$. Standard expressions exist for completing this transformation (e.g., Witt 1977, their Equations (22) and (23)).

Refraction can be included in the Monte Carlo simulation. As the photon moves along s, the path can be sub-divided, and a curvature applied to each smaller ${\rm{\Delta }}s$. At the photon height, the curvature, ${r}_{{\rm{c}}}$, is computed according to Equation (38), and a deflection angle for the photon is determined from the refractive portion of Equation (41) (i.e., ${\rm{\Delta }}{\theta }_{r}={\rm{\Delta }}s/{r}_{{\rm{c}}}$). The change in trajectory is referenced from $\hat{{\boldsymbol{\mu }}}$ and occurs in either the local upward or downward direction (depending on the sign of ${r}_{{\rm{c}}}$). A new set of direction cosines is found by requiring,

$$\begin{eqnarray}&&\hat{{\boldsymbol{\mu }}}\cdot \hat{{\boldsymbol{\mu }}}^{\prime} =\cos ({\rm{\Delta }}{\theta }_{r}),\end{eqnarray} \tag{ 17 }$$

$$\begin{eqnarray}&&(\hat{{\boldsymbol{\mu }}}\times \hat{{\boldsymbol{r}}})\cdot \hat{{\boldsymbol{\mu }}}^{\prime} =0,\end{eqnarray} \tag{ 18 }$$

and

$$\begin{eqnarray}&&| \hat{{\boldsymbol{\mu }}}^{\prime} | =1,\end{eqnarray} \tag{ 19 }$$

where the first equation implies that the trajectory changes through ${\rm{\Delta }}{\theta }_{r}$, the second equation forces the trajectory change to be either locally up or down, and the final equation ensures that the new direction of travel is a unit vector.

The computational efficiency of the Monte Carlo approach to multiple scattering in transit spectra outlined above can be increased in several ways. First, for a given impact parameter, the Monte Carlo routine need only be executed if the straight-line scattering optical depth is a substantial fraction of the absorption optical depth. (We choose a conservative cutoff at ${\tau }_{\lambda ,{\rm{s}}}\lt {10}^{-3}{\tau }_{\lambda ,{\rm{a}}}$.) Also, as suggested by de Kok & Stam (2012), the number of photons used in the Monte Carlo simulation at a given impact parameter should be influenced by the straight-line transmission—these authors recommend the use of 10⁵ photons in transparent conditions and up to 10⁸ photons in opaque conditions. However, de Kok & Stam (2012) used a traditional Monte Carlo radiative transfer approach, where photons are lost to absorption while passing through the atmosphere, thus driving the need for very large numbers of photons in their simulations. As we have separated absorption from scattering in our model, (nearly) all of the photons we simulate in our Monte Carlo exit the atmosphere, so that we only require 10⁴–10⁶ photons. Even more dramatic efficiency gains can be made by noting that, since our Monte Carlo model only considers scattering optical depths, these simulations need only be recomputed if either the layer differential scattering optical depths or asymmetry parameters change significantly. Thus, the individual path distribution matrices for each of the photons in our simulations can be saved and reused over a wide range of wavelengths in a high-resolution spectral grid.

Of course, running ${N}_{{\rm{p}}}$ transmission calculations (i.e., using Equation (8)) at every wavelength in a high-resolution grid becomes computationally unfeasible if the grid is large. Here, runtime reductions of factors of 10–100 can be straightforwardly made by computing the analytic Jacobian of Equation (9), $\partial {\bar{{\boldsymbol{t}}}}_{\lambda }/\partial {\boldsymbol{\Delta }}{{\boldsymbol{\tau }}}_{\lambda ,a}$, with,

$$\begin{eqnarray}\begin{array}{rcl}{\left(\displaystyle \frac{\partial {\bar{{\boldsymbol{t}}}}_{\lambda }}{\partial {\boldsymbol{\Delta }}{{\boldsymbol{\tau }}}_{\lambda ,a}}\right)}_{i,j} & = & \displaystyle \frac{1}{{N}_{{\rm{p}}}}\displaystyle \sum _{m=1}^{{N}_{{\rm{p}}}}\displaystyle \frac{\partial {t}_{\lambda ,m}({b}_{i})}{\partial {\rm{\Delta }}{\tau }_{\lambda ,a,j}}\\ & = & -\,\displaystyle \frac{1}{{N}_{{\rm{p}}}}\displaystyle \sum _{m=1}^{{N}_{{\rm{p}}}}{{ \mathcal P }}_{m,i,j}\cdot {t}_{\lambda ,m}({b}_{i})\end{array}.\end{eqnarray} \tag{ 20 }$$

This Jacobian can be used to rapidly adapt ${\bar{{\boldsymbol{t}}}}_{\lambda }$ to changes in ${\boldsymbol{\Delta }}{{\boldsymbol{\tau }}}_{\lambda ,a}$, and our tests indicate that this approach is accurate with up to 25%–50% variations in level-dependent absorption optical depths. Since a set of path distributions, transmissions, and Jacobians apply to wavelengths where the atmospheric optical property profiles are similar, the Monte Carlo approach outlined above is best implemented within a spectral mapping model (West et al. 1990; Meadows & Crisp 1996).

We note that the many efficiency gains outlined above stem from our separate treatments of scattering (which influences the path distributions) and absorption (which influences the path-derived transmissions and Jacobians). So far as we know, this treatment is a novel approach to Monte Carlo radiative transfer. Additionally, our approach to multiple scattering can be applied beyond exoplanet transits, and, for example, could be used to study scattering effects in occultation observations (i.e., an extension of de Kok & Stam 2012).

2.2. Generalized Transit Spectra

The transit spectra expressions given in Equations (3), (6), and (7) are primarily useful in the geometric limit, since rays in refracting or scattering models are not guaranteed to trace back to the stellar disk. To efficiently model spectra in these latter cases, we generalize the concept of the vector of annulus areas, ${\boldsymbol{A}}$. Key variables used in the discussion below are visualized in Figure 3.

Zoom In Zoom Out Reset image size

Ultimately, a transit spectrum is determined by comparing intensities integrated over the range of solid angles influenced by the planet (${{\rm{\Omega }}}_{{\rm{p}}}$) in the case when the planet is present versus when only the star is considered (with corresponding intensities ${I}_{{\rm{p}},\lambda }$ and ${I}_{{\rm{s}},\lambda }$, respectively). Thus,

$$\begin{eqnarray}&&{\left(\displaystyle \frac{{R}_{{\rm{p}},\lambda }}{{R}_{{\rm{s}}}}\right)}^{2}=\displaystyle \frac{{\int }_{{{\rm{\Omega }}}_{{\rm{p}}}}{I}_{{\rm{s}},\lambda }d{\rm{\Omega }}-{\int }_{{{\rm{\Omega }}}_{{\rm{p}}}}{I}_{{\rm{p}},\lambda }d{\rm{\Omega }}}{{\int }_{{{\rm{\Omega }}}_{{\rm{s}}}}{I}_{{\rm{s}},\lambda }d{\rm{\Omega }}},\end{eqnarray} \tag{ 21 }$$

where ${{\rm{\Omega }}}_{{\rm{s}}}$ is the range of solid angles corresponding to the stellar disk. Given the large distance D between Earth and any exoplanetary system, so that $d{\rm{\Omega }}={dA}/{D}^{2}$, and taking the stellar intensity to be normalized such that,

$$\begin{eqnarray}&&{I}_{0,\lambda }=\displaystyle \frac{{\int }_{{{\rm{\Omega }}}_{{\rm{s}}}}{I}_{{\rm{s}},\lambda }d{\rm{\Omega }}}{\ \ {{\rm{\Omega }}}_{{\rm{s}}}},\end{eqnarray} \tag{ 22 }$$

then,

$$\begin{eqnarray}&&{\left(\displaystyle \frac{{R}_{{\rm{p}},\lambda }}{{R}_{{\rm{s}}}}\right)}^{2}=\displaystyle \frac{{\int }_{{A}_{{\rm{p}}}}{I}_{{\rm{s}},\lambda }/{I}_{0,\lambda }{dA}-{\int }_{{A}_{{\rm{p}}}}{I}_{{\rm{p}},\lambda }/{I}_{0,\lambda }{dA}}{\pi {R}_{{\rm{s}}}^{2}}.\end{eqnarray} \tag{ 23 }$$

The first integral in the numerator of Equation (23) is simply a stellar limb darkening law integrated over the portion of the planetary disk that overlaps the stellar disk. Using a polar coordinate system centered on the planetary disk, and letting d be the distance from the center of the star to the center of the planet, this integral is then,

$$\begin{eqnarray}&&{\displaystyle \int }_{{A}_{{\rm{p}}}}\,\displaystyle \frac{{I}_{{\rm{s}},\lambda }}{{I}_{0,\lambda }}{dA}\\ \quad =\,\left\{\begin{array}{ll}0, & d\geqslant {R}_{{\rm{s}}}+{R}_{{\rm{p}}}+{h}_{{\rm{t}}}\\ {\displaystyle \int }_{{r}_{0}}^{{R}_{{\rm{p}}}+{h}_{{\rm{t}}}}{\displaystyle \int }_{-{\theta }_{0}(r)}^{{\theta }_{0}(r)}\displaystyle \frac{{I}_{{\rm{s}},\lambda }}{{I}_{0,\lambda }}({\mu }_{{\rm{s}}}){rd}\theta {dr}, & d\lt {R}_{{\rm{s}}}+{R}_{{\rm{p}}}+{h}_{{\rm{t}}},\end{array}\right.\end{eqnarray} \tag{ 24 }$$

where ${\mu }_{{\rm{s}}}={\mu }_{{\rm{s}}}(r,\theta )$ is the cosine of the angle of incidence on the stellar disk at coordinates r and θ, and

$$\begin{eqnarray}{r}_{0}=\left\{\begin{array}{ll}0, & d\leqslant {R}_{{\rm{s}}}\\ d-{R}_{{\rm{s}}}, & d\gt {R}_{{\rm{s}}},\end{array}\right.\end{eqnarray} \tag{ 25 }$$

$$\begin{eqnarray}{\theta }_{0}(r)=\left\{\begin{array}{ll}\pi , & d\leqslant {R}_{{\rm{s}}}-{R}_{{\rm{p}}}-{h}_{{\rm{t}}}\\ \tfrac{{d}^{2}+{r}^{2}-{R}_{{\rm{s}}}^{2}}{2{dr}}, & d\gt {R}_{{\rm{s}}}-{R}_{{\rm{p}}}-{h}_{{\rm{t}}}.\end{array}\right.\end{eqnarray} \tag{ 26 }$$

However, as limb darkening is accounted for in standard transit observation data reduction procedures (e.g., Mandel & Agol 2002; Kreidberg 2015), one may wish to omit the limb darkening law from the first integral in the numerator of Equation (23), leaving an analytic integral that can be performed by considering the geometry of two overlapping disks.

The second integral in the numerator of Equation (23) contains the details associated with absorption, scattering, and refraction in the planetary atmosphere. We define $\tilde{I}(r,\theta )$ as the background surface brightness mapped onto by an area element on the planet at r and θ. This function is zero for rays that do not map back to the stellar disk, and, in the geometric limit, is equal to ${I}_{{\rm{s}},\lambda }({\mu }_{{\rm{s}}}(r,\theta ))$ for portions of the planetary disk that overlap the star. Given this definition, we then have,

$$\begin{eqnarray}&&{\int }_{{A}_{{\rm{p}}}}\,\displaystyle \frac{{I}_{{\rm{p}},\lambda }}{{I}_{0,\lambda }}{dA}=\displaystyle \frac{1}{{I}_{0,\lambda }}{\int }_{{R}_{{\rm{p}}}}^{{R}_{{\rm{p}}}+{h}_{{\rm{t}}}}{\int }_{-\pi }^{\pi }{t}_{\lambda }(r)\tilde{I}(r,\theta ){rd}\theta {dr},\end{eqnarray} \tag{ 27 }$$

where we have assumed that ${R}_{{\rm{p}}}$ is set sufficiently deep so that ${I}_{{\rm{p}},\lambda }=0$ for $r\lt {R}_{{\rm{p}}}$. Of course, it is straightforward to consider thermal emission from the planetary disk, which can be important at longer wavelengths (Kipping & Tinetti 2010), by including the planet-to-star flux ratio for the planetary nightside on the right-hand side of Equation (23).

For efficient implementation, we define a grid on the atmospheric portion of the planetary disk in ${N}_{\theta }$ angular and N_r radial points. An ${N}_{r}\times {N}_{\theta }$ area matrix, with ${A}_{i,j}={r}_{i}{\rm{\Delta }}{r}_{i}{\rm{\Delta }}{\theta }_{j}$, need only be computed once, and the ${N}_{r}\times {N}_{\theta }$ background surface brightness mapping matrix, $\tilde{{\boldsymbol{I}}}$, is computed alongside the path distributions and transmissions. Given these, we have

$$\begin{eqnarray}&&{\int }_{{A}_{{\rm{p}}}}\,\displaystyle \frac{{I}_{{\rm{p}},\lambda }}{{I}_{0,\lambda }}{dA}=\displaystyle \frac{1}{{I}_{0,\lambda }}{{\boldsymbol{t}}}_{\lambda }\cdot [(\tilde{{\boldsymbol{I}}}\,\circ \,{\boldsymbol{A}})\cdot 1],\end{eqnarray} \tag{ 28 }$$

where "$\circ$" indicates the Hadamard product, $1$ is a vector of ones with length ${N}_{\theta }$, and, as before, ${{\boldsymbol{t}}}_{\lambda }$ is computed using either Equation (5) or Equations (9) and (20). Critically, this approach still efficiently isolates the parameters that vary rapidly in wavelength, ${{\boldsymbol{t}}}_{\lambda }$, from the parameters influenced by geometry and ray paths—the background surface brightness mapping matrix need only be computed once in the geometric limit, and 1–2 times when refraction is included. For the multiple scattering case, full Monte Carlo simulations need only be run at wavelength scales over which the atmospheric scattering properties vary, and transmissions need only be recomputed at wavelength scales over which their linear corrections via the analytic Jacobians become inaccurate. Integration can also be made more efficient by noting the symmetry about the line connecting the planet center to the star center.

The matrix $\tilde{{\boldsymbol{I}}}$ is straightforward to compute in the geometric limit, as rays trace straight lines back to either the stellar disk or not. When refraction is included, as in the Appendix, a ray tracing is performed once for each gridpoint in N_r ray impact parameters. Then, given the exit location and direction for these rays, and the angular location, ${\theta }_{j}$, for each area element, simple geometry indicates which of the ${N}_{\theta }$ gridpoints have rays that map back to the stellar disk (and what the ${\mu }_{{\rm{s}}}$ value is for each of these rays). In the multiple scattering case, ${N}_{{\rm{p}}}$ Monte Carlo simulations are run, and $\tilde{{\boldsymbol{I}}}$ is built up by averaging over the photons in each instance. For each photon that exits the atmosphere, its trajectory leaving the top of the atmosphere is investigated to see if the stellar disk is intersected. The geometry is then rotated through ${\rm{\Delta }}{\theta }_{j}$, and the intersection is reinvestigated. Thus, only ${N}_{{\rm{p}}}$ Monte Carlo simulations are needed to build up the average ${N}_{r}\times {N}_{\theta }$ background mapping matrix.

Figure 4 shows the surface brightness in a full Monte Carlo simulation prior to integration over solid angle. A haze with vertical optical depth of ${\tau }_{s}=0.01$ is placed above the 1 μbar pressure level following a linear power-law in pressure. The atmosphere is isothermal at 1500 K, the stellar and planetary size are appropriate for HD 189733b, and refraction and Rayleigh scattering are included (at 0.8 μm wavelength). A secondary figure enhances the surface brightness off the stellar disk by a factor of 10. For the portion of the planet that is not overlapping the stellar disk, the sub-figure with enhanced surface brightness shows a thin ring due to forward scattered light. Additional surface brightness structure in this ring is related to the scattering phase function—the region of the planetary limb that is at the greatest angular separation from the stellar disk less strongly samples the haze forward scattering peak.

Zoom In Zoom Out Reset image size

2.3. Model Summary

We validated our transit model, which we implement within the framework of the Spectral Mapping Atmospheric Radiative Transfer (SMART) model (developed by D. Crisp; Meadows & Crisp 1996), using a number of techniques and data sources. First, with refraction and scattering optical depth omitted, we verified that our ray tracing routine (discussed in the Appendix) and the Monte Carlo approach return the path distribution in the geometric limit. Then, with refraction included and scattering optical depth omitted, we showed that the ray tracing and Monte Carlo routines are in agreement. Finally, to check our Monte Carlo approach with scattering included, we made the simple switch to a plane-parallel geometry (where height is only measured relative to the z-axis, not radially) and validated the output radiances against results from the DISORT radiative transfer model (Stamnes et al. 1988) for a wide range of conditions. Figure 5 shows the results from two of these experiments, plotting the top-of-atmosphere intensity (scaled by the incident flux) as a function of observer zenith angle. Clouds with a given scattering optical depth and asymmetry parameter were placed over a perfectly absorbing surface, and a Henyey–Greenstein phase function was assumed.

Zoom In Zoom Out Reset image size

To further validate our model, we performed model inter-comparisons. First, using a standard Earth atmospheric model (McClatchey et al. 1972), we compared our transit spectrum for an Earth–Sun twin system to the refracting model of Misra et al. (2014). Note that the Misra et al. (2014) model has been extensively validated against solar occultation data for Earth (Gunson et al. 1996; Irion et al. 2002) and observations of transmission through Earth's atmosphere during a lunar eclipse (Pallé et al. 2009). Figure 6 shows the result of this inter-comparison, where agreement (in effective transit height, equal to ${R}_{{\rm{p}},\lambda }-{R}_{{\rm{p}}}$) is to within the atmospheric model grid spacing. Agreement improves if we use the coarser integration path lengths adopted by Misra et al. (i.e., 5 km).

Zoom In Zoom Out Reset image size

Second, and finally, we applied our transit spectrum tool to a standardized hot Jupiter-like atmospheric model widely used for inter-comparison purposes. This model has a planetary radius (at the 10 bar pressure level) of 1.16 ${R}_{{\rm{J}}}$, a planetary mass of 1.14 ${M}_{{\rm{J}}}$, a stellar radius of 0.78 R_⊙, and atmospheric volume mixing ratios of H₂, He, and H₂O of 0.85, 0.15, and $4\times {10}^{-4}$, respectively. The 126 model layers are placed evenly in log-pressure between 10 bar and ${10}^{-9}\,\mathrm{bar}$. The atmosphere is isothermal at 1500 K. Figure 7 compares our model, in the geometric limit, to output from the CHIMERA retrieval suite (Line et al. 2013) for this standard case. Agreement is within the 20 ppm scatter seen in comparisons between other (Irwin et al. 2008; Waldmann et al. 2015) transit spectrum tools (M. R. Line 2017, private communication). We note that the offset in the Rayleigh scattering slopes between our model and the CHIMERA calculation is due to differing approaches to computing Rayleigh scattering opacities which results in optical depths that differ at the 5% level or less.

Zoom In Zoom Out Reset image size

We use the formalism outlined above to, first, compute the path distribution for Earth-like and hot Jupiter-like atmospheres for a range of conditions. Following these instructive examples, we explore how scattering influences the transit depth associated with a single aerosol layer over a wide range of parameter space. Finally, we demonstrate the effects of thin, scattering clouds on a typical hot Jupiter transit spectrum.

4.1. Path Distributions

The aforementioned standard Earth and hot Jupiter atmospheric models can be used to demonstrate the path distribution for a range of conditions and assumptions. To show the effects of clouds on the path distribution for the hot Jupiter atmosphere, we use the "cloud" and "haze" differential optical depth profiles shown in Figure 8. The former is based on the shape of the aerosol mixing ratio profiles for a condensational cloud from the Ackerman & Marley (2001) model, while the latter follows a power law in altitude with scale height equal to the pressure scale height. Both models are taken to have integrated vertical optical depth unity, and to be conservatively and moderately forward scattering ($g=0.9$). A Henyey–Greenstein phase function is assumed for cloudy and hazy simulations.

Zoom In Zoom Out Reset image size

Figure 9 demonstrates the path distribution in the geometric limit for the standard hot Jupiter-like model atmosphere. Altitudes are measured relative to the 10 bar radius (at $1.16{R}_{{\rm{J}}}$). This graph is interpreted by selecting an impact parameter, and using the color shading to interpret the path distribution (which, as stated earlier, can be thought of as an enhancement of the vertical optical depth). For example, taking an impact parameter "altitude" ($b-{R}_{{\rm{p}}}$) of 3000 km, the path distribution is zero for all atmospheric layers below this height, since rays with this impact parameter never pass through these deeper atmospheric layers. The path distribution is then largest for layers near an altitude of 3000 km, as rays pass through these layers either horizontally or nearly horizontally. Finally, for atmospheric layers well above 3000 km, rays pass through on much less extreme slant paths, so the path distribution is smaller (typically 5–15, for this particular model atmosphere).

Zoom In Zoom Out Reset image size

The path distribution considering refraction for our cloud-free Earth atmosphere is shown in Figure 10. For comparison, a case in the geometric limit is also presented. For elements of the path distribution where the impact parameter altitude is near the surface, the path distribution deviates substantially from the geometric case due to the refractive bending of ray paths. As an example, take an impact parameter altitude of 5 km. In the geometric case, rays with this impact parameter would never reach altitudes below 5 km. But, due to refraction, these rays are bent "downward" to probe deeper atmospheric layers. Thus, the path distribution for the refracting case is now non-zero for atmospheric layers below 5 km. Furthermore, the atmospheric layer that sees the most enhancement of optical depth (i.e., has the largest path distribution) is now slightly below 5 km, instead of being right at 5 km as in the geometric case.

Zoom In Zoom Out Reset image size

Figures 11 and 12 show an average path distribution from our multiply-scattering Monte Carlo approach with

$$\begin{eqnarray}&&\bar{{\boldsymbol{ \mathcal Q }}}=\displaystyle \frac{1}{{N}_{{\rm{p}}}}\displaystyle \sum _{m=1}^{{N}_{{\rm{p}}}}{{\boldsymbol{ \mathcal Q }}}_{m}.\end{eqnarray} \tag{ 29 }$$

We note that such an average path distribution is not used in our calculation of transit spectra, but is simply meant to indicate the characteristic path that scattered photons take through the atmosphere. In the Figure 11 we only consider molecular Rayleigh scattering opacity (at 0.55 μm wavelength) in our hot Jupiter-like model atmosphere. For Figure 12, we show the effects of our nominal haze and cloud models for the hot Jupiter-like case. While both the haze and cloud have the same integrated optical depth, the haze path distribution has a less pronounced transition into the aerosol-affected atmospheric layers since the condensate cloud is concentrated into a more narrow range of altitudes/pressures. Note the use of a logarithmic color contour scale to emphasize the impacts of scattering.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The path distributions for the scattering cases are less straightforward to interpret. For the Rayleigh scattering case (Figure 11), and taking an impact parameter altitude of 800 km, we see that the path distribution is non-zero for altitudes below 800 km (but would be zero in the geometric case). Here, a small fraction of photons have been scattered down to these deeper layers, and the small values of the path distribution at depth come at the (small) expense of the path distribution aloft. Selecting a deeper impact parameter altitude of 200 km, we see that the path distribution is now greatest at large altitudes (where the path distribution would have been roughly 10 in the geometric case), and is small at altitudes below about 400–600 km. Here, photons with a 200 km impact parameter altitude are scattered at altitudes larger than about 400–600 km, so that relatively few of these photons can probe much deeper than this.

The haze case in Figure 12 can be interpreted similarly to the Rayleigh scattering case, although the range of impact parameters that are influenced by the haze is greater than that of the Rayleigh scattering case. The condensate cloud case in Figure 12 appears distinct due to the well-defined cloud top. Here, rays/photons with impact parameter altitudes larger than 1400 km never encounter the cloud, and the path distribution is the same as in the geometric case. At altitudes just below this, though, the photons encounter the cloud, and are scattered before reaching atmospheric layers below the cloud (whose base is at about 1200 km). For these photons, their path distribution is concentrated at altitudes above the cloud and in the cloud itself.

4.2. Exploration of the Importance of Scattering

As most exoplanet transit spectrum models do not include multiple scattering, a key exercise is to explore the range of conditions for which scattering is expected to influence transit depth. For this exploration, we note that the transit depth attributed to a narrow annulus on the planetary disk (i.e., at a single ray impact parameter) scales with $2\pi b{\rm{\Delta }}b/{R}_{{\rm{s}}}^{2}$. Thus, in the absence of limb darkening, with the atmosphere transparent at radii larger than $b+{\rm{\Delta }}b$, and with the planetary disk centered on the star, only three parameters will influence the annulus transit depth for a scattering scenario—the layer scattering optical depth, the aerosol scattering phase function, and the angular size of the stellar disk from the orbital location of the planet.

Figure 13 shows the relative difference between a transit spectrum model in the geometric and pure absorption limits (where scattering optical depth is treated as absorption optical depth) and a multiple scattering model for a single annulus on the planetary disk. Results are given for three different values of the host star angular size, which is determined by the ratio of the stellar radius to the planetary orbital distance (a). These values range from a host star angular size like that for Mercury and the Sun ($a=0.39\,\mathrm{au}$, so that ${R}_{{\rm{s}}}/a\sim {10}^{-2}$) to that of a hot Jupiter-like planet orbiting at 0.05 au from a Sun-like star (${R}_{{\rm{s}}}/a\sim {10}^{-1}$). Also, this range of host star angular sizes spans those for the Habitable Zones of M and K dwarf stars (Kopparapu et al. 2013).

Zoom In Zoom Out Reset image size

The color contours in Figure 13 are given as a function of layer scattering slant optical depth and the scattering asymmetry parameter. No other opacity source is added to the layer. For simplicity, a Henyey–Greenstein phase function is assumed (Henyey & Greenstein 1941) as using a multi-parameter phase function would introduce additional variables to our phase space exploration. Regardless, these results will still serve to indicate under which conditions scattering can become important.

4.3. Thin Clouds and Hazes in Hot Jupiter Transits

The standard hot Jupiter-like atmospheric model outlined in Section 3 provides a useful test case for exploring some important aspects of scattering in exoplanet transit spectra, especially when considering that planets near to their host stars are those for which scattering can have the largest effects. The parameter space for such an exploration is extremely large, since the wavelength-dependent asymmetry parameter, scattering optical depth, and the cloud/haze vertical distribution, will all influence the transit spectrum. Exploring all of this phase space is certainly beyond the scope of this manuscript, so we adopt a straightforward set of conditions. Specifically, we investigate the impact of general, isolated aerosol layers, placed at different pressure levels in the model atmosphere. Rather than adopt a specific cloud or haze optical depth profile (e.g., Figure 8), we simply distribute the aerosol optical depth uniformly over a single pressure scale height. Scattering optical depths and asymmetry parameters are assumed gray, the single-scattering albedo is taken to be unity (i.e., pure scattering clouds), and, as before, a Henyey–Greenstein scattering phase function is adopted.

Figure 14 shows transit spectra for strongly forward scattering cloud/haze particles ($g=0.95$) over the 1–2 μm wavelength range, which overlaps the accessible wavelength range for Hubble/WFC3 (Kimble et al. 2008). Figure 15 is similar, but for less strongly forward scattering cloud/haze particles ($g=0.90$). Single-layer clouds were placed at pressures of 10⁻⁴, 10⁻³, 10⁻², and 10⁻¹ bar, and spanned a vertical range of $d\mathrm{ln}p=1$. The clearsky case is shown in black, while pure absorbing and scattering cases are shown in gray and purple, respectively. For each cloud pressure, three scenarios with different slant scattering optical depths are shown, with ${\tau }_{{\rm{s}}}$ equal to 1, 10, and 100. This range of scattering optical depths spans the limits from a relatively optically thin aerosol layer to a relatively optically thick layer.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The path distribution approach provides a coherent and computationally efficient method to computing transit spectra. Furthermore, analyzing the path distribution can provide an understanding of which atmospheric layers can contribute information to rays emerging from the planetary disk at a given impact parameter. For example, Figure 10 demonstrates how refraction bends rays to probe deeper atmospheric layers than would be encountered in the geometric limit. However, since the path distribution is intrinsic to just the planetary atmosphere, the other important effect of refraction—bending rays such that they do not strike the stellar disk—is not represented. This is a strength of the path distribution approach, since only a single distribution would need to be computed for identical atmospheric models for planets orbiting different host stars.

Similarly, Figure 11 shows how, on average, Rayleigh scattering will scatter a small fraction of photons at large impact parameters "downward" to probe deeper parts of the atmosphere. Conversely, rays at small impact parameters no longer probe the deep atmosphere, as they are scattered at larger altitudes. Of course, as Rayleigh scattering is not strongly forward scattering, it is unlikely that any significant fraction of the scattered photons will still intersect the stellar disk.

Figure 12 shows similar scattering effects to those in Figure 11. However, the impact parameter where (roughly) deeper atmospheric layers become ineffectively probed depends on the vertical structure of the aerosols and, more specifically, where the slant scattering optical depth unity occurs. We note, again, that these figures display average path distributions for our Monte Carlo simulations, and that transit spectra that include scattering must be computed using using the formalism developed in Section 2.1.

Ultimately, when considering scattering, the brightness of any given annulus (at impact parameter b) on the planetary disk will depend on the transparency of atmospheric layers at radii larger than b, the slant scattering optical depth of the layer at b, the scattering phase function for the aerosols in this layer, and the angular size of the host star as seen from the planet. For the situation where overlaying layers are transparent (implying that we can, thus, see down to a cloud), Figure 13 provides a rough guide to the circumstances where scattering can be important. Here, one can use simple details from a forward model—the angular size of the host star, the scattering optical thickness of a cloud layer, and cloud scattering asymmetry parameter—to determine if scattering can safely be ignored in a simulation.

In the case where the host star angular size is much smaller than it would be at a distance of roughly 0.4 au from a Sun-like star (i.e., ${R}_{{\rm{s}}}/a\sim {10}^{-2}$), Figure 13 shows that essentially any amount of scattering will prevent the path of a ray from tracing back to the stellar disk (which is a very small "target"). With ${R}_{{\rm{s}}}/a$ at roughly (1–5) $\times \,{10}^{-2}$, which, for example, is appropriate for the Habitable Zones of late M dwarf stars, strongly forward scattering aerosols ($g\gt 0.9$) can have significant impacts on the transit depth due to an annulus for clouds whose slant optical depths are not very optically thick (i.e., ${\tau }_{s}$ less than 10). This effect becomes quite dramatic for hot Jupiter-like conditions, where the host star is a relatively large "target." Here, even modestly forward scattering aerosols ($g\sim 0.8$) can cause substantial variations as compared to the commonly used pure absorption assumption. This dependence on the asymmetry parameter explains why both Hubbard et al. (2001) and Brown (2001) found that multiple scattering was unimportant for their presented hot Jupiter transit spectra. As these authors considered only either Rayleigh scattering or isotropic phase functions (which both have g = 0), scattering could effectively be treated as absorption.

Our adoption of the Henyey–Greenstein phase function is driven by the computational simplicity of this parameterization and the need to capture the process of forward scattering while minimizing the number of added parameters to our simulations. Nevertheless, it has been shown that the Henyey–Greenstein phase function, when compared to Mie calculations, can underestimate the power in the forward scattering peak (Toublanc 1996; Boucher 1998), especially at wavelengths comparable to the particle size. Thus, depending on the specifics of a given aerosol in an exoplanetary atmosphere, the results shown in Figure 13 may underestimate the increased transmission due to forward scattered light at a given asymmetry parameter (as compared to the pure absorption limit).

Figures 14 and 15 further explore the importance of scattering in hot Jupiter transit spectra. For high-altitude clouds/hazes, these results show that differences between a pure absorption model versus a scattering model can approach (or exceed) 200 ppm, which is larger than the typical uncertainties achieved in current observations with Hubble (Kreidberg et al. 2014) or that are expected for JWST (Greene et al. 2016). Very optically thick clouds (i.e., with slant scattering optical depths approaching 100) approach the pure absorption limit, and the difference between a scattering model and a pure absorption model for a thin cloud (i.e., scattering optical depth less than about unity) is small since the overall impact of the cloud on the transit spectrum is weak. Similarly, the impact of scattering for a deep cloud or haze layer is small, as few wavelength regions are sensitive to the presence of this aerosol layer.

In general, the impact of aerosol multiple scattering will depend on cloud configurations and optical properties, and this manuscript does not undertake the complex microphysical calculations needed to attempt predictions of cloud/haze distributions and compositions in exoplanet atmospheres—our goal is to highlight the conditions where scattering becomes an important consideration. Nevertheless, the results in Figures 14 and 15 indicate that forward scattering exoplanet clouds analogous to Earth's cirrus clouds (which are at higher altitudes and have optical depths less than roughly unity) may warrant a multiple scattering approach. Additionally, any forward scattering cloud with a scale height comparable to (or greater than) the pressure scale height could require a multiple scattering model to properly simulate a transit spectrum, as transit depth variations due to multiple scattering effects in this clouds will be comparable to those of molecular features (whose scale are typically set by the pressure scale height).

Differences between the scattering and pure absorption models for the hot Jupiter cases will also depend on the asymmetry parameter and single-scattering albedo. For very strongly forward scattering aerosols (i.e., asymmetry parameters larger than our adopted value of 0.95), the transit spectrum will approach the clearsky limit (if the aerosols are weakly or non-absorbing). In essence, such clouds would be invisible, having no effect on the transit depth. On the other hand, asymmetry parameters much smaller than roughly 0.8–0.9 will push the scattering transit spectra toward the pure absorption limit, since these scattered photons are increasingly unlikely to connect to the observer. Similarly, as the single-scattering albedo is decreased (away from our adopted value of unity), the scattering spectra will approach the pure absorption case. The extent of the influence of the single-scattering albedo depends on the average number of scatterings the photons experience along their path. For example, for our ${\tau }_{{\rm{s}}}=10$ case, the photons are scattered of order ten times, implying a single-scattering albedo of less than about 0.9 would reduce the transmission to roughly 50%.

Perhaps the most important message from the scattering studies above is the influence of the scattering asymmetry parameter on transit spectra of close-in exoplanets. Aerosols whose scattering properties vary from strongly forward scattering to weakly forward scattering could impart features in a transit spectrum, as the cloud (or haze layer) will be more transparent where the particles have a larger asymmetry parameter. Such signatures would appear along with absorption bands caused by particle vibrational modes (Wakeford & Sing 2015) and cloud base features (Vahidinia et al. 2014). Additionally, as transit spectra of forward scattering clouds would be reproduced by relatively thinner clouds in the pure absorption limit, transit spectral retrievals using pure absorption models could return cloud thickness (or number densities) that are biased low, where retrieved cloud optical depths would, instead, represent an effective optical depth that incorporates information about the optical properties of the cloud particles and the stellar angular size.

Whether or not aerosol multiple scattering will be important for any given exoplanet transit spectrum will depend on a number of parameters, including the slant optical thickness of any haze/cloud structure and the scattering asymmetry parameter (which will both, in general, depend on wavelength). The planet must also be located near to its host star, which limits the range of relevant aerosols to those which form in warm and hot atmospheric conditions. For example, asymmetry parameters for water droplets in the visible wavelength range are 0.8–0.9 (Kokhanovsky 2004), and span 0.74–0.94 for ice crystals (where crystals with plate-like morphologies are the most strongly forward scattering; Macke et al. 1998). Thus, conditions may be appropriate for multiple scattering to influence the transit spectra of potentially habitable worlds around the coolest stars (e.g., TRAPPIST-1 Gillon et al. 2016). Looking beyond Earth in the Solar System, de Kok & Stam (2012) (their Figure 1) show that a variety of ices, dusts, droplets, and fractal hazes are both non-absorptive and forward scattering (with $g\gtrsim 0.9$) below 2–3 μm. The wavelength range for these Solar System cases is especially relevant to transiting exoplanets around later-type stars that may be investigated by Hubble, JWST, the Fast Infrared Exoplanet Spectroscopy Survey Explorer concept (FINESSE; Deroo et al. 2012), and/or the Atmospheric Remote-Sensing Infrared Exoplanet Large-survey concept (ARIEL; Tinetti et al. 2016).

The composition, size distributions, and optical properties of aerosols in warm and hot exoplanet atmospheres are largely unknown. For the hottest exoplanets, metal and silicate condensates may be expected to form (Lunine et al. 1989; Marley et al. 1999; Burrows et al. 2001). Budaj et al. (2015) investigated the optical properties of a large variety of metallic and silicon-bearing aerosols, and showed that many of these species can be strongly forward scattering over certain wavelength ranges. Especially relevant cases are aluminum oxide, forsterite, enstatite, and pyroxine.

It is interesting to note that our results show that multiple scattering can be an important consideration for exoplanet transits when the angular size of the host star is relatively large, which is the opposite regime from where refraction has been shown to be important (Bétrémieux & Kaltenegger 2014, 2015; Misra et al. 2014). For refraction, a host star with small angular size implies that relatively little refractive bending is required to deflect a ray off the stellar disk, thereby setting a floor in the transit spectrum which is not captured in the commonly used geometric, pure absorption limit. By comparison, only a single scattering (with even a strong forward peak) is required to deflect a ray from the stellar disk, which is then in agreement with the geometric, pure absorption limit. Of course, the opposite of these statements is true for a star with large angular size. In future work we plan to further explore the importance of refraction in transit observations of gas giants.

Finally, since this manuscript primarily emphasizes the development of a light scattering theory for exoplanet transits and the impacts of vertically thin cloud layers, a number of future studies will be undertaken. These will include scattering effects in more vertically extended cloud structures, and could also investigate the role of the assumed (or computed) particle scattering phase function. Finally, the formalism outlined above can be used to study light scattering effects in time-resolved transit light curves.

We have detailed a new theory of exoplanet transit spectroscopy that includes the effects of light multiple scattering. By effectively separating the path that photons take through an exoplanet atmosphere from the absorption processes of that atmosphere, the technique discussed in this manuscript yields models that are both physically rigorous and computationally efficient. This approach, which relies on the so-called "path distribution" defined herein, can be extended to other areas of study, including (most straightforwardly) stellar occultations by planetary atmospheres.

When applying our validated scattering model to isolated cloud layers, models show that multiple scattering is most important for cases where the exoplanet host star is large in angular size as seen from the world (i.e., when the orbital distance is less than 10–20 times the stellar radius). In these cases, multiple scattering by aerosols with asymmetry parameters larger than 0.8–0.9 can have substantial effects on the transmission of the cloud layer. For all cases, differences between a multiple scattering model and a geometric pure absorption diminish for clouds (or hazes) with slant scattering optical depths approaching 100.

In an exploratory case of a conservatively and forward (g = 0.90–0.95) scattering cloud/haze layer in the atmosphere of a hot Jupiter, differences in the transit depth from a multiple scattering model and a model in the pure absorption limit can exceed 200 ppm. These differences are most pronounced when the cloud/haze is well above the level of gas absorption optical depth unity, and when the slant scattering optical depth is of order several to several tens. Future work will further explore the situations where scattering in more extended exoplanet cloud (or haze) layers is key to understanding transit observations.

T.R. gratefully acknowledges support from NASA through the Sagan Fellowship Program executed by the NASA Exoplanet Science Institute. The results reported herein benefited from collaborations and/or information exchange within NASA's Nexus for Exoplanet System Science (NExSS) research coordination network sponsored by NASA's Science Mission Directorate. Certain essential tools used in this work were developed by the NASA Astrobiology Institute's Virtual Planetary Laboratory, supported by NASA under Cooperative Agreement No. NNA13AA93A. The author thanks M. Line and J. Lustig-Yaeger for sharing model outputs for inter-comparison purposes, and J. Fortney and R. de Kok for providing helpful feedback. This manuscript is dedicated to the memory of Prince, who passed away during the completion of this project—"thank you for a funky time."

Software: SMART (Meadows & Crisp 1996), DISORT (Stamnes et al. 1988), CHIMERA (Line et al. 2013).

Given the ray atmospheric path distribution as defined in Section 2, a number of key quantities can be derived. First, the pathlength (in km, for example) traversed between any two altitudes, h₁ and h₂, by a ray incident with impact parameter b is simply

$$\begin{eqnarray}&&{s}_{b}\,({h}_{1},{h}_{2})=\left|{\int }_{{h}_{1}}^{{h}_{2}}{{ \mathcal P }}_{b}\,(h){dh}\right|,\end{eqnarray} \tag{ 30 }$$

where the absolute value forces all distances to be measured non-negative. Then, if h = 0 at the bottom of a terrestrial planetary atmosphere or at some reference pressure (e.g., 10 bar) for a gaseous world, and ${h}_{{\rm{t}}}$ is the effective top of the atmosphere, the total linear distance traversed through the entire atmosphere by the ray is

$$\begin{eqnarray}&&{s}_{b}\,(0,{h}_{{\rm{t}}})={\int }_{0}^{{h}_{{\rm{t}}}}{{ \mathcal P }}_{b}\,(h){dh},\end{eqnarray} \tag{ 31 }$$

where we note that this integral would diverge if taken to infinite altitude.

The pathlength integral can also be written in terms of two pressure coordinates, p₁ and p₂, assuming the ideal gas law and hydrostatic equilibrium

$$\begin{eqnarray}{s}_{b}\,({p}_{1},{p}_{2}) & = & \left|{\displaystyle \int }_{{p}_{2}}^{{p}_{1}}{{ \mathcal P }}_{b}\,(p)\displaystyle \frac{{R}_{{\rm{g}}}T\,(p)}{{mg}}d\mathrm{ln}p\right|\\ & = & \left|{\displaystyle \int }_{{p}_{2}}^{{p}_{1}}{{ \mathcal P }}_{b}\,(p)H\,(p)d\mathrm{ln}p\right|,\end{eqnarray} \tag{ 32 }$$

where ${R}_{{\rm{g}}}$ is the universal gas constant, $T(p)$ is the atmospheric temperature profile, g is the acceleration due to gravity, m is the atmospheric mean molecular weight, and $H(p)$ is the pressure scale height. Also, given the definition of the path distribution, and the differential definition of (path) column mass from the atmospheric mass density profile $[\rho (h)$],

$$\begin{eqnarray}&&d{{ \mathcal M }}_{{\rm{c}}}=\rho (h){ds},\end{eqnarray} \tag{ 33 }$$

the air mass encountered by a ray with impact parameter b between two altitudes, h₁ and h₂, is given by

$$\begin{eqnarray}&&{{ \mathcal M }}_{{\rm{c}},b}\,({h}_{1},{h}_{2})=\left|{\int }_{{h}_{1}}^{{h}_{2}}\rho (h){{ \mathcal P }}_{b}\,(h){dh}\right|.\end{eqnarray} \tag{ 34 }$$

Similarly, a column (path) number density is

$$\begin{eqnarray}&&{{ \mathcal N }}_{{\rm{c}},b}\,({h}_{1},{h}_{2})=\left|{\int }_{{h}_{1}}^{{h}_{2}}n(h){{ \mathcal P }}_{b}\,(h){dh}\right|,\end{eqnarray} \tag{ 35 }$$

where $n(h)$ is the number density profile.

In the geometric limit, rays pass straight through the planetary atmosphere. Here, the path distribution can be determined analytically using geometric arguments. For an atmospheric layer centered at h with width ${\rm{\Delta }}h$, and for a ray incident on the atmosphere with impact parameter b, the path distribution is given by

$$\begin{eqnarray}{{ \mathcal P }}_{b}\,(h)=\,\left\{\begin{array}{ll}0, & b\geqslant r+{\rm{\Delta }}h/2\\ \tfrac{2}{{\rm{\Delta }}h}\sqrt{{(r+{\rm{\Delta }}h/2)}^{2}-{b}^{2}}, & r-{\rm{\Delta }}h/2\lt b\lt r+{\rm{\Delta }}h/2\\ \tfrac{2}{{\rm{\Delta }}h}[\sqrt{{(r+{\rm{\Delta }}h/2)}^{2}-{b}^{2}}-\sqrt{{(r-{\rm{\Delta }}h/2)}^{2}-{b}^{2}}], & b\leqslant r-{\rm{\Delta }}h/2,\end{array}\right.\end{eqnarray} \tag{ 36 }$$

where we have defined $r={R}_{{\rm{p}}}+h$ for conciseness. Note that this expression is only independent of ${\rm{\Delta }}h$ in the limit that this value is small when compared to r. Performing a first-order expansion in ${\rm{\Delta }}h$ yields

$$\begin{eqnarray}&&{{ \mathcal P }}_{b}\,(h)\approx \displaystyle \frac{2r}{\sqrt{{r}^{2}-{b}^{2}}},\end{eqnarray} \tag{ 37 }$$

which is in agreement with the linearized geometry discussed in Fortney (2005) (their Figure 1). The geometric path distribution need only be computed once for a model atmosphere. Thus, when paired with Equation (7), the geometric approach can be executed with great computational efficiency.

When considering refraction, rays no longer travel on straight-line paths through the atmosphere. In this case, the path distribution must be computed numerically. Fortunately, as the refractive indexes for gases likely to be major atmospheric constituents vary weakly in wavelength, the path distribution need only be computed at a small number of spectral points when generating a transit spectrum.

Refraction will deflect the trajectory of a ray upward or downward, depending on the sign of the local atmospheric refractive index profile. Here, the local curvature experienced by the ray, ${r}_{{\rm{c}}}$, is

$$\begin{eqnarray}&&\displaystyle \frac{1}{{r}_{{\rm{c}}}}=\sin ({\theta }_{r})\displaystyle \frac{d\mathrm{ln}{n}_{\mathrm{ref}}}{{dh}},\end{eqnarray} \tag{ 38 }$$

where ${\theta }_{r}$ is the zenith angle for the ray trajectory, and ${n}_{\mathrm{ref}}$ is the atmospheric refractive index. The zenith angle is determined via,

$$\begin{eqnarray}&&\cos ({\theta }_{r})=\hat{{\boldsymbol{r}}}\cdot \hat{{\boldsymbol{\mu }}}.\end{eqnarray} \tag{ 39 }$$

Using pathlength as the integration variable, and according to the geometry shown in Figure 2, we have

$$\begin{eqnarray}&&\displaystyle \frac{{dh}}{{ds}}=\cos ({\theta }_{r}),\end{eqnarray} \tag{ 40 }$$

$$\begin{eqnarray}&&\displaystyle \frac{d{\theta }_{r}}{{ds}}=-\left[\displaystyle \frac{\sin ({\theta }_{r})}{r}+\displaystyle \frac{1}{{r}_{{\rm{c}}}}\right],\end{eqnarray} \tag{ 41 }$$

$$\begin{eqnarray}&&\displaystyle \frac{d\phi }{{ds}}=\displaystyle \frac{\sin ({\theta }_{r})}{r},\end{eqnarray} \tag{ 42 }$$

and

$$\begin{eqnarray}&&\displaystyle \frac{d\omega }{{ds}}=\displaystyle \frac{1}{{r}_{{\rm{c}}}},\end{eqnarray} \tag{ 43 }$$

where ϕ is the polar angle measured from the $+\hat{{\boldsymbol{x}}}$ direction, and ω is the so-called refraction integral (which measures the deflection from the initial straight-line trajectory).

For a ray with initial impact parameter b_i, the path is determined by numerically integrating the expressions above using an algorithm like that described in van der Werf (2008). While the ray is passing through the atmospheric layer at h_j, the increments of ${\rm{\Delta }}s$ used in the path integration are divided by ${\rm{\Delta }}{h}_{j}$ and added to ${{ \mathcal P }}_{i,j}$. The path integration proceeds through all atmospheric layers until either the ray exits the atmosphere or the planetary surface is struck. When the ray exits the atmosphere, the radial coordinate and trajectory are stored for later use in determining which rays map back to the stellar disk. The straight-line trajectory the ray follows after exiting the atmosphere is defined by the total refraction angle at exit, which is equal to the direction cosine of the ray along the x-axis in Figure 2 (i.e., $\hat{{\boldsymbol{\mu }}}\cdot \hat{{\boldsymbol{x}}}={\mu }_{x}=\omega$).