Exoplanetary Monte Carlo Radiative Transfer with Correlated-k I. Benchmarking Transit and Emission Observables

TL;DR

This paper benchmarks a 3D Monte Carlo radiative transfer model for exoplanet atmospheres, demonstrating its accuracy and applicability for detailed atmospheric characterization including clouds and complex geometries.

Contribution

It introduces and validates the CMCRT 3D Monte Carlo radiative transfer code, comparing its results to existing models and observational data, highlighting its suitability for complex exoplanet atmosphere simulations.

Findings

CMCRT matches transmission spectra benchmarks within tens of ppm.

Emission spectra benchmarks agree within 10% of 1D models.

The method effectively models inhomogeneous atmospheres with clouds and multiple scattering.

Abstract

Current observational data of exoplanets are providing increasing detail of their 3D atmospheric structures. As characterisation efforts expand in scope, the need to develop consistent 3D radiative-transfer methods becomes more pertinent as the complex atmospheric properties of exoplanets are required to be modelled together consistently. We aim to compare the transmission and emission spectra results of a 3D Monte Carlo Radiative Transfer (MCRT) model to contemporary radiative-transfer suites. We perform several benchmarking tests of a MCRT code, Cloudy Monte Carlo Radiative Transfer (CMCRT), to transmission and emission spectra model output. We add flexibility to the model through the use of k-distribution tables as input opacities. We present a hybrid MCRT and ray tracing methodology for the calculation of transmission spectra with a multiple scattering component. CMCRT compares well to the transmission spectra benchmarks at the 10s of ppm level. Emission spectra benchmarks are consistent to within 10% of the 1D models. We suggest that differences in the benchmark results are likely caused by geometric effects between plane-parallel and spherical models. In a practical application, we post-process a cloudy 3DHD 189733b GCM model and compare to available observational data. Our results suggest the core methodology and algorithms of CMCRT produce consistent results to contemporary radiative transfer suites. 3D MCRT methods are highly suitable for detailed post-processing of cloudy and non-cloudy 1D and 3D exoplanet atmosphere simulations in instances where atmospheric inhomogeneities, significant limb effects/geometry or multiple scattering components are important considerations.

Tables

Table 1. Table 1: Gas phase absorbers and Rayleigh scattering species used as opacity sources in this study. The current HELIOS-K tables do not include Na and K opacities.

Species	Reference: NEMESIS - HELIOS-K
H2O	Barber et al. (2006) - Polyansky et al. (2018)
CH4	Yurchenko & Tennyson (2014)
CO	Rothman et al. (2010) - Li et al. (2015)
CO2	Tashkun & Perevalov (2011) - Rothman et al. (2010)
NH3	Yurchenko et al. (2011)
PH3	Sousa-Silva et al. (2015)
Na	Heiter et al. (2008) - N/A
K	Heiter et al. (2008) - N/A
H2-H2 CIA	Richard et al. (2012)
H2-He CIA	Borysow et al. (2001); Borysow (2002)
Scattering
H2	Irwin (2009)
He	Irwin (2009)

Equations

$$\psi (x_0)=\frac{{\int }_a^{x_0}p(x)dx}{{\int }_a^{b}p(x)dx},$$

$$g_{samp}=\frac{w_g}{{\sum }_g{w}_g}.$$

$${\left(\frac{R_p}{R_{\star }}\right)}^2=\frac{1}{R_{\star }^2}\left(R_{p,0}^2+2\sum _{i=1}^{N_i}[1-⟨T_i⟩]b_i\Delta b_i\right),$$

$$T_{i,ph}=\exp (-{\tau }_{i,ph}),$$

$$T_{i,ph}=\omega W_{ph}\Phi (\alpha )\exp (-{\tau }_{i,ph}),$$

$$⟨T_i⟩=\frac{\sum T_{i,ph}}{N_{i,ph}}.$$

$$L_{i,\lambda }=4\pi {\rho }_i{V}_i{\kappa }_{i,\lambda }^{abs}{B}_{\lambda }(T_i),$$

$$L_{i,g}=4\pi {\rho }_i{V}_i{\kappa }_{i,g}^{abs}{B}_{\lambda }(T_i),$$

$$L_{i,\Delta \lambda }=\sum _g{w}_g{L}_{i,g}.$$

$$g_{samp}=\frac{w_g{L}_{i,g}}{{\sum }_g{w}_g{L}_{i,g}},$$

$$L_{p,\Delta \lambda }=\sum _i{f}_i{L}_{i,\Delta \lambda }.$$

Full text

We aim to compare the transmission and emission spectra results of a 3D Monte Carlo Radiative Transfer (MCRT) model to contemporary radiative-transfer suites.

We perform several benchmarking tests of a MCRT code, Cloudy Monte Carlo Radiative Transfer (CMCRT), to transmission and emission spectra model output. We add flexibility to the model through the use of k-distribution tables as input opacities. We present a hybrid MCRT and ray tracing methodology for the calculation of transmission spectra with a multiple scattering component.

CMCRT compares well to the transmission spectra benchmarks at the 10s of ppm level. Emission spectra benchmarks are consistent to within 10% of the 1D models. We suggest that differences in the benchmark results are likely caused by geometric effects between plane-parallel and spherical models. In a practical application, we post-process a cloudy 3D HD 189733b GCM model and compare to available observational data.

Our results suggest the core methodology and algorithms of CMCRT produce consistent results to contemporary radiative transfer suites. 3D MCRT methods are highly suitable for detailed post-processing of cloudy and non-cloudy 1D and 3D exoplanet atmosphere simulations in instances where atmospheric inhomogeneities, significant limb effects/geometry or multiple scattering components are important considerations.

††pubyear: 2019††pagerange: Exoplanetary Monte Carlo Radiative Transfer with Correlated-k I. Benchmarking Transit and Emission Observables–Exoplanetary Monte Carlo Radiative Transfer with Correlated-k I. Benchmarking Transit and Emission Observables

1 Introduction

Observing and characterising the three-dimensional atmospheric structure of exoplanets is a continuing endeavour for the exoplanetary community. A keystone tool for investigating exoplanet atmospheres in forward and retrieval efforts is radiative-transfer modelling, examining how radiation interacts with the atmosphere and gives rise to the observable properties of each individual exoplanet.

Exoplanet atmospheres are continuing to be observed in ever greater detail. Transmission spectroscopy from the ground (e.g. Chen et al., 2017; Gibson et al., 2017; Kirk et al., 2017) and space (e.g. Sing et al., 2016) has shown many hot and warm Jupiters and Neptunes to contain a variety of molecular and atomic species, for example, Na (Charbonneau et al., 2002), K (Sing et al., 2011), H2O (Swain et al., 2009), CH4 (Swain et al., 2008) and CO (Snellen et al., 2010), revealing the atmospheric composition at the transmission limbs of the exoplanet atmosphere (e.g. Barstow et al., 2017; Tsiaras et al., 2018; Fisher & Heng, 2018). Many exhibit evidence of cloud coverage at the terminator limbs due to an observed optical wavelength Rayleigh-like slope (e.g. Pont et al., 2013), muted IR water vapour features (e.g. Deming et al., 2013) or consistent with a flat, featureless (e.g. Kreidberg et al., 2014; Wakeford et al., 2017) grey spectra. Emission (e.g. Knutson et al., 2012) and reflection (e.g. Evans et al., 2013) observations have provided information on the dayside temperature structure, and suggest the presence in some exoplanetary atmospheres of a significant optical wavelength scattering component larger than the gas phase Rayleigh scattering component.

3D modelling efforts using Global Circulation Models (GCMs) (e.g. Showman et al., 2009; Rauscher & Menou, 2012; Dobbs-Dixon & Agol, 2013; Mayne et al., 2014; Mendonça et al., 2016) of hot Jupiters and Neptunes suggest atmospheric inhomogeneities in temperature, velocity fields, and vertical mixing rates (e.g. Parmentier et al., 2013; Charnay et al., 2015a; Kataria et al., 2016). Gas phase chemistry modelling in 1D (e.g. Drummond et al., 2016; Tsai et al., 2017; Blumenthal et al., 2018) and 2D/3D (e.g. Agúndez et al., 2012; Drummond et al., 2018a, b; Mendonça et al., 2018; Steinrueck et al., 2018) suggest that considering kinetic non-equilibrium chemistry on the 3D chemical composition of exoplanet atmospheres is important to the resulting temperature structures and observational properties of each individual exoplanet.

Cloud modelling efforts in 1D (e.g. Morley et al., 2013; Helling et al., 2016; Lavvas & Koskinen, 2017; Ohno & Okuzumi, 2017; Gao et al., 2018; Gao & Benneke, 2018; Powell et al., 2018) and 3D (e.g. Charnay et al., 2015b; Lee et al., 2016; Parmentier et al., 2016; Lewis et al., 2017; Roman & Rauscher, 2017; Lines et al., 2018b; Roman & Rauscher, 2019) of a variety of sophistications also suggest inhomogeneous cloud coverage in latitude, longitude and depth across the globe of their atmospheres.

With the advent of space missions with dedicated exoplanet atmospheric characterisation payloads, such as the James Webb Space Telescope (JWST) (e.g. Bean et al., 2018), the Wide Field Infra-Red Space Telescope (WFIRST) (e.g. Robinson et al., 2016), and the Atmospheric Remote-sensing Infrared Exoplanet Large-survey (ARIEL) (Tinetti et al., 2016), the need for accurate radiative-transfer modelling for these 3D inhomogeneous objects will become increasingly pertinent for the physical interpretation of exoplanetary observables. Recently, Feng et al. (2016) and Blecic et al. (2017) have shown the impacts of considering 3D temperature structures when interpreting the retrieval results of emission spectra data.

Through these observational and modelling efforts it is clear that exoplanet atmospheres are inherently 3D and inhomogeneous in nature, with unique local dynamical, temperature, chemical and cloud properties. Due to the stochastic and microphysical nature of the Monte Carlo Radiative Transfer (MCRT) technique, it is well suited to modelling this complex 3D environment. In the astrophysical community, MCRT modelling techniques have been proven to be highly successful and constitute the majority of in use 3D radiative-transfer models (Steinacker et al., 2013). For exoplanets, MCRT has been applied for calculations of transmission spectra (de Kok & Stam, 2012; Robinson, 2017), geometric albedo and albedo spectra (Hood et al., 2008; Garcia Munoz & Isaak, 2015; García Muñoz et al., 2017) plus emission spectra (Lee et al., 2017; Stolker et al., 2017). In addition, modelling photon processes such as polarization (Stolker et al., 2017) and atmospheric refraction (Robinson, 2017) can be readily included.

In this study, we systematically compare and benchmark the results of the 3D MCRT code Cloudy Monte Carlo Radiative Transfer (CMCRT) presented in Lee et al. (2017) to several contemporary radiative transfer codes, in particular the NEMESIS radiative transfer suite (Irwin et al., 2008), the 3D MCRT model ARTES (Stolker et al., 2017) and the Baudino et al. (2017) benchmark protocols. In Section 2, we outline updates to the Lee et al. (2017) model and describe aspects of using k-distribution opacities with MCRT. Sections 3 and 4 detail our approaches for calculating transmission spectra and emission spectra using MCRT respectively. In Sect. 5 we present a direct comparison of the MCRT results to the NEMESIS suite. In Sect. 5.2 we present the benchmarking to the emission spectra results of Stolker et al. (2017). Section 5.3 presents the benchmarking to the Baudino et al. (2017) protocols. Section 6 applies our new methodologies to post-process output of the cloudy HD 189733b GCM simulation of Lee et al. (2016). Section 7 contains the discussion and Sect. 8 contains the conclusions.

2 MCRT Methods Using K-distributions

We update and expand the MCRT model described in Lee et al. (2017), based on the original work by Hood et al. (2008). In Lee et al. (2017) pre-mixed Sharp & Burrows (2007) mean absorption coefficient tables were used. The MCRT methods are generalised to include the use of individual and mixes of gas species with k-distribution tables for application in exoplanet atmosphere radiative transfer problems.

In the k-distribution method, a high fidelity wavelength or wavenumber absorption coefficient table is re-ordered by its cumulative distribution counterpart for a pre-defined bin width. A number of points are then sampled from the cumulative distribution function (k-coefficients), given by the values of a Gaussian ordinance (g-ordinance). Each point is then assigned a weight given by the g-ordinance method used. However, wavelength information is scrambled in the process, leading to the assumption, for an inhomogeneous media, that the k-coefficients at each Gaussian ordinance are correlated. This assumption leads to an error when modelling an inhomogeneous atmosphere, as the opacity distribution shifts temperature and pressure dependently. This error has been shown to retain accuracy to within 10% when compared to line-by-line tests (e.g. Amundsen et al., 2014).

The advantage of the k-distribution method is that it significantly reduces the computational burden by orders of magnitude compared to line-by-line calculations (Heng, 2017), and is a standard methodology in the planetary science community. More in-depth descriptions of k-distributions and correlated-k in exoplanet contexts can be found in Amundsen et al. (2014); Grimm & Heng (2015); Heng (2017) and Amundsen et al. (2017). When required, we apply the random overlap method (e.g. Lacis & Oinas, 1991; Amundsen et al., 2017) to combine individual gas k-tables weighted by their relative abundances.

Due to the statistical nature of both MCRT and k-coefficients, they have very complimentary properties. In MCRT, each packet is evolved in the simulation by sampling a normalised cumulative distribution function (CDF), $\psi$, of a probability function, $p(x)$, from the general equation (e.g. Stolker et al., 2017)

[TABLE]

where $x_{0}$ is a randomly sampled variable, with $a$ and $b$ the lower and upper limit of the distribution.

MCRT can make use of k-tables by sampling the weighted cumulative distribution properties of the k-distribution method. After a packet is spawned in the simulation for a given wavelength bin, a g-ordinance can be randomly selected for that packet by sampling the cumulative distribution function of the Gaussian quadrature weights, wg, where the probability of sampling a specific g-ordinate, $g$, is

[TABLE]

In the MCRT simulation, a uniformly sampled random number, $\zeta$ $\in$ [0,1], is drawn for each packet which corresponds to a single g-ordinance value. This g-ordinance value is then retained throughout the packet’s lifetime, and used for all further calculations involving the packet. The correlated-k approximation (e.g. Heng, 2017) is therefore inherent since we are assuming that the g-ordinance sampled for the packet is correlated across the entire 3D atmosphere.

To visualise the sampling of the g-ordinance, Fig. 1 shows examples of the g values and weights from numerous correlated-k approaches in the literature:

• •

The uniform 20 ordinance points typically used by the NEMESIS (Irwin et al., 2008) radiative-transfer suite .

• •

The 4+4 used in the SPARC/MITgcm (Showman et al., 2009) GCM model.

• •

The 8+8 used in the Exo-REM (Baudino et al., 2015, 2017) radiative-convective equilibrium model.

In typical use, the additional higher order ordinance points are used in GCM and radiative-convective modelling due to their larger bin sizes (e.g. Showman et al., 2009; Amundsen et al., 2014). This is performed in an attempt to capture the opacity contribution of numerous line centers in each bin, while keeping the efficiency of the radiative-transfer scheme reasonable. For retrieval and post-processing efforts with smaller bin sizes (hence less numbers of line centers), a uniform spacing with a larger number of g-ordinances is sufficient to capture the opacity distribution well. In traditional correlated-k radiative transfer codes, the Gaussian quadrature procedure is applied to calculate integrated mean quantities for a given wavelength bin.

In the context of MCRT, Eq. 2 shows that larger weighted g-ordinance values are more likely to be sampled by the scheme. Since the weights are normally distributed, central g-ordinance points are more likely to be sampled compared to those at the distribution wings. The MCRT scheme therefore acts as a numerical integrator by sampling the weighted distribution and summing the results of each sample at the simulation end. As with any Monte Carlo integration calculation, less variance in the final answer is obtained by increasing the number of samples. In principle, any number or sampling scheme of g-ordinances with weights can be used in the MCRT scheme. Using less g-ordinances will not necessarily decrease the computational run-time in this model, since the total number of sampled packets primarily controls the run-time. Individual wavelength calculations are retained by using a single g-ordinate with unity weight.

3 Transmission Spectra

Stellar light travelling through a planetary atmosphere principally interacts through the extinction of photons in the line of sight direction. This extinction reduces the average transmission of stellar light through the atmosphere and imprints atomic and molecular features in the transmission spectra. In addition, photons may also undergo multiple scattering interactions in the atmosphere, and subsequently escape towards the line of sight. Multiple scattered photons exiting the atmosphere towards the line of sight will increase the average transmission compared to extinction alone, resulting in the planet appearing smaller than it would without a multiple scattering component. This brightening effect from multiple scattered photons was first studied in the context of exoplanet transmission spectra by Hubbard et al. (2001).

The transmission spectrum formula in discretised form is given by (e.g. Dobbs-Dixon & Agol, 2013; Robinson, 2017)

[TABLE]

where $R_{p}$ [cm] is the apparent radius of the planet, $R_{\star}$ [cm] the radius of the host star, $R_{p,0}$ [cm] defined as the radius below which the planet can be considered an opaque solid body, $\langle T_{i}\rangle$ the mean transmission at impact parameter index $i$, $b_{i}$ [cm] the height of the impact parameter chord and $\Delta b_{i}$ [cm] the radial width of transit chord $i$.

Central to the CMCRT transmission spectra model is calculating the transmission through the atmosphere at each transit chord including the effects of multiple scattering. The sequence of events for a packet in the transmission spectrum calculation is as follows:

A packet is spawned at a random impact parameter at the transmission annulus of the 3D spherical grid. 2. 2.

The packet is assigned a g-ordinate value, by randomly sampling the cumulative k-distribution weights (Eq. 2). 3. 3.

An initial transmission calculation towards the observational direction is performed. 4. 4.

The packet is evolved through the simulation grid until termination or escape. 5. 5.

The packet’s contribution to the transmission from scattering events is recorded throughout its lifetime through the use of the next event estimation method.

Figure 2 (LHS) shows a diagrammatic representation of the path of a photon packet and the key steps considered in the simulation.

In our hybrid scheme, we combine the MCRT with a ray tracing method, the next event estimation technique (Yusef-Zadeh et al., 1984; Wood & Reynolds, 1999), to calculate the contribution of each packet to the transmission through each impact parameter. The initial transmission contribution of a single packet, $T_{i,ph}$, initially at a random location across the impact parameter width of vertical index $i$ is

[TABLE]

where $\tau_{i,ph}$ is the total optical depth towards the simulation boundary in the observational direction. At each scattering location for the packet during its lifetime, the scattering event also contributes an additional transmission to the packet’s current impact parameter index $i$ given by

[TABLE]

where $\omega$ is the local single scattering albedo, $W_{ph}$ is the current weight of the packet and $\Phi(\alpha)$ the normalised scattering phase function probability towards observation direction $\alpha$. It is important that this scattering contribution is calculated at the impact parameter of the current 3D location of the packet, which may be different from its originally initialised impact parameter after many scattering events. The mean transmission through impact parameter index $i$, $\langle T_{i}\rangle$, is then the total transmission from all contributing packets (initial plus scattering contributions), normalised by the number of packets that was originally randomly initialised across the impact parameter index, $N_{i,ph}$,

[TABLE]

This mean transmission at each impact parameter index is then used in Eq. (3) to calculate the radius ratio.

In summary, all simulated packets contribute a transmission through an impact parameter vertical index $i$, which is tracked for all packets throughout the simulation runtime. Since we are using the next event estimation method, the end points of the packets do not determine the transmission spectrum, and every packet contributes to the final solution. This increases the overall efficiency, and helps lower the variance of the model (e.g. see discussion in Lee et al. (2017)). We also apply the survival biasing with Russian Roulette scheme (e.g. Dupree & Fraley, 2002; Lee et al., 2017) in order to more accurately capture the effects of multiple scattering on the transmission spectra.

A feature of the model is that should the atmosphere have zero scattering opacity, or packet scattering be turned off, the scheme reduces to a 3D extinction limit ray tracing model. The impact of individual scattering components on the end spectra can be examined by running the simulation with and without packet scattering.

4 Emission Spectra

For emission spectra calculations we use the same scheme as in Lee et al. (2017) but with some additional properties when using k-distribution tables. The monochromatic luminosity, $L_{i,\lambda}$ [erg s*-1* cm*-1*], of cell $i$ is given by (e.g. Pinte et al., 2006; Stolker et al., 2017)

[TABLE]

where $\rho_{i}$ [g cm*-3*] is the cell gas density, $V_{i}$ [cm3] the cell volume, $\kappa_{i,\lambda}^{abs}$ [cm2 g*-1*] the mass absorption opacity (possibly including cloud particle absorption opacity) and $B_{\lambda}(T_{i})$ [erg s*-1* sr*-1* cm*-3*] the Planck function at temperature $T_{i}$ [K].

Using k-distribution coefficients, the luminosity contribution from each g-ordinate in the wavelength bin is given by

[TABLE]

where $\lambda$ [cm] is taken as the bin center wavelength, with the total pseudo-spectral luminosity in a band $\Delta\lambda$ given by

[TABLE]

The probability of sampling a particular g-ordinate for the photon packet in each cell is now given by

[TABLE]

which has the important property that the emission spectra in each wavelength bin is more strongly determined by the higher order g-ordinances, corresponding to greater sampling of (near) line-centers in the wavelength bin. The total luminosity of the planet at a given wavelength bin, $L_{p,\Delta\lambda}$ [erg s*-1*], is the sum of the luminosity of each cell, multiplied by the fraction of energy that escaped from that cell, $f_{i}$, toward the observational direction through the next event estimation method (Lee et al., 2017),

[TABLE]

The sequence of events for an emission spectra packet is as follows:

A packet is spawned at a random starting position within a cell volume. 2. 2.

A g-ordinance is sampled for the packet given by Eq. (10). 3. 3.

An initial transmission calculation towards the observational direction is performed. 4. 4.

The packet is evolved through the simulation until termination or escape. 5. 5.

The fraction of energy escaping towards the observational direction through scattering events is tracked throughout the packet’s lifetime through the next event estimation method.

Figure 2 (RHS) shows a diagrammatic representation of the emission scheme in CMCRT.

An assumption in this approach is that Eq. (8) implies that the value of $B_{\lambda}(T_{i})$ is constant across the wavelength bin range, and equal to the bin center wavelength. This approximation is reasonable when considering the small bin sizes ($<$ 0.05 $\mu$m) of the current post-processing efforts. However, this approximation may be invalid for large bin sizes where the Planck function would vary significantly across the bin edges.

Lastly, we note that our emission spectra model accounts for the wavelength dependent photospheric radius effect discussed in Fortney et al. (2019). The wavelength dependence of the photospheric radius is self-consistency accounted for in the 3D grid, since the fraction of a cell’s luminosity, directly weighted by the available emitting material in the cell (Eq. 7), escaping toward the observational direction is calculated. Additionally, any variation in the emitting area as a function of planetary orbital phase is readily accounted for.

5 Benchmarking Tests

We perform all simulations on a 3D spherical grid with longitude and latitude sizes ($\theta$, $\phi$) = (141, 61) corresponding to a typical exoplanet GCM simulation grid resolution (e.g. Dobbs-Dixon & Agol, 2013; Mayne et al., 2014; Kataria et al., 2016). The radial grid size is variable depending on the benchmark protocols. All temperature-pressure profiles used in this study are presented in Fig. 3.

The MCRT model uses a radial height based spherical coordinate grid for all calculations. We apply the height output from the NEMESIS code directly as the grid for the NEMESIS benchmark tests. For all other tests, a second order hydrostatic calculator (benchmarked to the NEMESIS output) is used to calculate the altitude given a temperature-pressure structure.

For input gas phase opacities we use k-distribution tables produced by the NEMESIS (Irwin et al., 2008) and HELIOS-K (Grimm & Heng, 2015) opacity calculators, dependent on the benchmark. All benchmarks presented here were repeated using both sets of k-tables, with only minor differences found between the output, mainly stemming from the different wavelength/wavenumber resolutions used for each opacity set. The input line-lists are presented in Table 1.

We use a custom made opacity tool which reads in individual k-tables, performs the T-p grid interpolation and random overlap k-table mixing, and is responsible for generating the opacity structure of the model atmosphere to be read in by CMCRT. This code also calculates any required gas phase Rayleigh scattering opacities, continuum opacities and cloud optical properties.

When applying the opacities produced by NEMESIS, we use an evenly spaced wavelength resolution of \Delta$$\lambda = 0.005 $\mu$m for 5929 wavelength bins between 0.305 $\mu$m and 29.945 $\mu$m. Each bin is sampled using 20 uniform g-ordinance points (Fig. 1). Several tests were conducted with 5$\times$ (\Delta$$\lambda = 0.001 $\mu$m) and a tenth (\Delta$$\lambda = 0.05 $\mu$m) this wavelength resolution for both the emission and transmission spectrum modes, all producing consistent results. This suggests a low-sensitivity of the results for these wavelength bin widths. For the HELIOS-K opacities, we use an evenly spaced grid of \Delta$$\nu = 10 cm*-1* between 0-30000 cm*-1* wavenumber. We use the Chebyshev polynomial opacity distribution fits from HELIOS-K (Grimm & Heng, 2015), sampling 20 uniform g-ordinance points for each bin, the same as the NEMESIS tables. We currently do not include Na and K gas phase opacities in the HELIOS-K benchmarks, as they are not yet included in the HELIOS-K database.

5.1 NEMESIS Benchmarking

In our first benchmark test, we perform a direct comparison between CMCRT and the NEMESIS radiative-transfer suite of Irwin et al. (2008). We use the same input values (temperatures, pressures, heights, opacities etc) as to allow a direct comparison between the core algorithms. For input values, the commonly used HD 189733b-like benchmark from Robinson (2017), with bulk planetary parameters $R_{p,0}$ (10 bar) = 1.16 $R_{j}$ (1 Rj = 6.9911 $\cdot$ 109 cm), Mp = 1.14 Mj (1 Mj = 1.89813 $\cdot$ 1030 g) is used. The atmospheric properties are 156 layers between 10*-9*-10 bar in log-space gas pressure, 1500 K isothermal gas temperatures and constant mean molecular weight $\bar{\mu}$ = 2.316 g mol*-1*. A constant volume mixing ratio of 0.85, 0.15 and 4 $\cdot$ 10*-4* is applied for H2, He and H2O respectively. The input opacity sources are H2-H2 and H2-He collisional induced absorption, H2 and He Rayleigh scattering and H2O absorption.

Figure 4 presents the results of the transmission spectra benchmark. CMCRT compares well to the NEMESIS output, with the differences between the models remaining below 30 ppm for the equivalent extinction case. A systematic positive offset of $\approx$10-30 ppm is seen across the infrared wavelengths. When multiple scattering is allowed in CMCRT, the increased transmission drops the optical Rayleigh slope by $\approx$ 20 ppm. Since NEMESIS includes multiple-scattering in its formalisms (Barstow et al., 2014), the better agreement with CMCRT in the multiple-scattering case is a promising indication that the scattering contributions are being accurately captured in CMCRT.

As a test of a more realistic atmospheric environment, we produced transmission and emission spectra of the Teff = 1500 K benchmark case from Sect. 5.3 (Fig. 5). The transmission spectra agree to within $\approx$10 ppm, with another small $\approx$5 ppm systematic positive offset in the infrared regions. The emission spectra also agrees well, with a relative difference to within 10% across the wavelength range. Slightly deeper Na, K, H2O and CO absorption features are produced by CMCRT compared to the NEMESIS spectra, suggesting a different spectral line formation region present in the 3D CMCRT grid.

The high frequency variations between the models are attributed to the Monte Carlo noise error of sampling the g-ordinance weights, evidenced by the direct imprint of absorption opacity features in the lower panels of Fig. 4 and Fig. 5. However, it is clear a small systematic offset is seen between the models for both the transmission and emission spectra results. In Sect. 7 we discuss possible geometric differences as an explanation for the systematic offsets.

5.2 Stolker et al. benchmark

In Stolker et al. (2017) several emission spectra benchmarks were performed using the 3D MCRT model ARTES, primarily to investigate emergent polarisation signatures from directly imaged exoplanets. We repeat the cloud free, self-luminous planet emission spectra protocol from Appendix B in Stolker et al. (2017), denoted by Teff = 400 K, 800 K and 1200 K respectively (Fig. 3). We run two simulations for each profile, one with the mean opacity tables directly used in Stolker et al. (2017) (T. Stolker priv. corr.) in individual wavelength mode and one using the HELIOS-K opacities in correlated-k mode.

Figure 6 presents a comparison between our results and Stolker et al. (2017). The simulations using the Stolker et al. (2017) opacity tables are in excellent agreement with the ARTES output, except at optical wavelengths in the Teff = 1200 K test case. The results using the HELIOS-k opacities show differences in the optical regime due to the non-inclusion of Na and K opacity, but also show a general positive offset with shallower absorption features in the infrared. We suggest that the differences in the results produced by CMCRT when using the Stolker et al. (2017) opacity tables and HELIOS-K tables stem from the choice of Voigt line wing cutoff (Sect. 7) used for the input opacities between the models. However, it is encouraging that CMCRT produces consistent results to the Stolker et al. (2017) spectra when using the ARTES opacities directly, as these two models share similar 3D Monte Carlo methodologies.

5.3 Baudino et al. Benchmarks

In Baudino et al. (2017), the radiative-convective models ATMO (Tremblin et al., 2015; Drummond et al., 2016; Goyal et al., 2018), Exo-REM (Baudino et al., 2015) and petitCODE (Mollière et al., 2015, 2017) were benchmarked across a variety of 1D temperature and pressure conditions suitable for exoplanet atmospheres. We benchmark to the transmission and emission spectra of the prescribed Guillot (2010) T-p profiles, denoted as Teff = 500 K, 1000 K, 1500 K, 2000 K and 2500 K (Fig. 3). Figures 7, 8, 9, 10 and 11 present the Teff = 500 K, 1000 K, 1500 K, 2000 K and 2500 K benchmark results respectively. For these tests we apply the NEMESIS k-distribution tables as the input opacities. All transmission and emission results from CMCRT are convolved to the benchmark wavelength resolution using the SpectRes package (Carnall, 2017).

5.3.1 Molecular abundances

In Baudino et al. (2017) the chemical equilibrium (CE) schemes of each model were compared and used as input gas phase abundances. We apply the publicly available CE with condensation code of Woitke et al. (2018), GGchem, to each T-p profile and reduce the full species database (Worters et al., 2017) to the gas and solid/liquid phase species listed in Baudino et al. (2017). Input elemental ratios at solar metallicity are taken from Asplund et al. (2009). As per the benchmark protocols, H3PO4[l] (Phosphoric acid) is added to the GGchem condensate list, and we additionally include the calculation of gas phase SiO since the vapour pressure expression from Wetzel et al. (2013) is used to calculate the SiO[s] supersaturation ratio. The corrected thermochemical data for PH3 from Lodders (1999) is also used.

The CE calculations compare well to the other codes, with differences occurring for PH3 mole fractions at high pressures in the $T_{\rm eff}$ = 1500 K, 2000 K and 2500 K profiles. A significant deviation between the calculations is seen in the $T_{\rm eff}$ = 2000 K and 2500 K profiles for $p_{\rm gas}$ $<$ 10*-3* bar where the H2O, CO, CO2, Na and K abundances drop off. We attribute this to the thermal disassociation and ionisation of these molecular and atomic species at lower pressure and higher ($T_{\rm gas}$ $>$ 2000 K) gas temperature atmospheric conditions. We note that considering thermal dissociation of molecules was not explicitly included in the Baudino et al. (2017) benchmark protocols.

5.3.2 Transmission spectra

The transmission spectra for each profile is typically consistent to the benchmark spectra to within the $\approx$ 10s of ppm level. CMCRT compares best with the ATMO results, but are generally slightly positively offset from the ATMO results. We suggest a 3D grid effect in our methodology is responsible for this positive offset (Sect. 7).

5.3.3 Emission spectra

The emission spectra are generally consistent with the benchmark results to within the $\approx$ 10% level. A small positive offset between CMCRT and the benchmark results is seen in the Teff = 500 K and 1000 K profile at the peak emission and the Rayleigh-Jeans tail wavelengths. We suggest a difference between the 3D spherical and plane-parallel models produces this offset (Sect. 7). As in Sect. 5.1, for the Teff = 1500 K profile a CO absorption feature is present in the CMCRT model between 4-5 $\mu$m not seen in the benchmark results. Since all the models in this test used the Rothman et al. (2010) CO line-list, this suggests that the CO line forming regions in the deep atmosphere are different in the 3D grid compared to the 1D codes for this profile.

6 Post-processing GCM output

As a 3D spherical grid model, CMCRT is well placed to accurately post-process 3D GCM output. To test this capability, we post-process output of the cloudy 3D GCM of HD 189733b presented in Lee et al. (2016) using our new methodologies. We examine the effect of the cloud structures by performing MCRT simulations using the same 3D GCM output with and without cloud opacity. For input to CMCRT, we use the individual cell temperatures, pressures, cloud properties (mean sizes, number density) and the gas phase element abundances involved in the cloud formation process. This is not a self-consistent comparison between a cloud-free and cloudy GCM simulation however, since the GCM thermal structures were produced with cloud opacity effects included. To increase the near-IR wavelength resolution, we use the currently available HELIOS-K produced gas phase opacities in this section, Na and K gas phase opacities are therefore neglected for this model. In order to focus on the thermal emission of the planet only, we do not include the contribution from reflected stellar light.

For gas phase abundances, the simulation is post-processed assuming chemical equilibrium using GGChem (Woitke et al., 2018), with the local depletion of elements from the cloud formation processes taken from the GCM results. For the cloud opacity, the simulation is post-processed using effective medium theory and Mie theory the same way as in Lee et al. (2017). To simplify the calculation, we assume all cloud opacity and scattering properties at the mean particle radius. We note this approximation has been investigated by Powell et al. (2018) who showed differences in cloud opacity and scattering properties when considering the integrated effect of a full cloud particle size distribution. Examining differences in the transmission and emission properties with a full particle size distribution is left to future efforts.

Figure 12 (top left) presents the transmission spectra results of the GCM post-processing with and without cloud opacity. The model transmission spectra is compared to published HST (Pont et al., 2013; McCullough et al., 2014) and Spitzer (Knutson et al., 2012, and references within) data. In the infrared, a broad absorption feature from 8 - 30 $\mu$m is produced due to the cloud particles dominant silicate composition at the limbs of the simulated atmosphere (e.g. Wakeford & Sing, 2015). The results suggest strong CH4 absorption features at $\approx$ 3.3 $\mu$m and 7.8 $\mu$m, however, these features may be modified by non-equilibrium chemistry such as quenching (e.g. Steinrueck et al., 2018), potentially altering the CH4 abundance at higher altitudes compared to the local chemical equilibrium assumption in this study. The GCM transmission spectra produces a flatter optical slope and very muted water features in the near-IR compared to the observational data for the cloudy case. This suggests that the cloud particles in the upper atmosphere of the GCM are too large compared to what fitting and retrieval to the HD 189733b observational data suggest (e.g. Wakeford & Sing, 2015; Barstow et al., 2017). A similar conclusions was discussed in Lines et al. (2018b, a) for microphysical cloud GCM modelling of HD 209458b.

Figure 12 (top right, bottom left) presents the dayside emission spectra results of the GCM post-processing with and without cloud opacity. The flux ratio results are compared to published HST (Crouzet et al., 2014; Barstow et al., 2014) and Spitzer (Knutson et al., 2012; Todorov et al., 2014, and references within) data. For the parent star luminosities we use a Castelli & Kurucz (2004) ATLAS9 stellar atmosphere model with parameters; Teff = 5000 K, $\log$ g = 4.5, [M/H] = 0.0). The results show that cloud coverage significantly dampens the optical and near-IR emission, generally lowering the emitted flux and increasing the strength of the absorption features across infrared wavelengths. Our results compare well to the HST data and photometric Spitzer data, however there is an offset between the Spitzer IRS analysis from Todorov et al. (2014). This suggests that the photospheric temperature as modelled on the dayside of the GCM may be slightly too high, and a cooler dayside temperature preferred, lowering the peak planetary emission in the Spitzer IRS wavelength range. A strong emission feature is present at $\approx$ 9.7 $\mu$m, attributed to the Si-O stretching mode from the dominant silicate cloud composition present in the GCM model (Lee et al., 2016). Due to the strength of this feature, it is quite possible it would have been detectable in the Spitzer IRS observations (Grillmair et al., 2008; Todorov et al., 2014). This suggests that the silicate emission feature is either muted or not present in the real object. However, this feature arises from the presence of silicates clouds at the temperature inversion regions at higher latitudes in the modelled atmosphere (Lee et al., 2016). The observation of a cloudy emission feature in an object in the future (e.g. JWST) may be evidence of a temperature inversion in the atmosphere and reveal the cloud composition.

Figure 12 (top right, bottom right) shows the nightside emission spectra. The nightside spectra show markedly different features compared to the dayside. The spectra is largely characterised by the strong CH4 absorption at $\approx$ 3.3 $\mu$m and 7.8 $\mu$m. In contrast to the dayside spectrum, the cloudy nightside produces an absorption feature at $\approx$ 9.7 $\mu$m and 15-20 $\mu$m corresponding to the Si-O stretching and bending modes of the silicates. This is reminiscent of the absorption feature observed by Spitzer IRS in some L dwarf spectra (Cushing et al., 2006). The nightside is generally consistent with the 3.6 $\mu$m and 4.5 $\mu$m flux ratio presented in Knutson et al. (2012), but is not able to reproduce the 8.0 $\mu$m and 24 $\mu$m points. This suggests the upper atmospheric temperatures in GCM on the nightside may be an underestimate.

7 Discussion

7.1 Line wing cutoff

A possible explanation for the offsets between the CMCRT and the benchmark emission spectra results in Sect. 5.2 and 5.3 is the treatment of the input opacity calculation. The NEMESIS k-tables used in this study applied a Voigt line wing cutoff of 25 cm*-1*, the HELIOS-K tables a cutoff of 100 cm*-1* and the Baudino et al. (2017) benchmark models use an adaptive cutoff (ATMO) or a sub-Lorentzian lineshape model (petitCODE and Exo-REM). The ARTES opacities used in Stolker et al. (2017) applied an infinite cutoff (Min (2017), M. Min priv. corr.). As shown in Baudino et al. (2017) the choice of line wing cutoff when calculating the input opacities can have significant effects on the end emission spectra results in the forward model.

In order to test the sensitivity of the model output to the line wing cutoff, we produce five H2O k-tables with the HELIOS-k code, using the Barber et al. (2006) line list with a variety of line wing cutoff prescriptions found in the literature:

25, 100 and infinite cm*-1* absolute cutoff. 2. 2.

Hybrid min[25 P (atm), 100] cm*-1* (Sharp & Burrows, 2007). 3. 3.

500 Lorentz widths (Grimm & Heng, 2015).

We use the same T-p profile and input quantities as the Teff = 1500 K Baudino et al. (2017) benchmark case, but with only H2O as the gas phase absorption opacity. We also process this T-p profile using the Stolker et al. (2017) opacity table.

Figure 9 presents the results of this test. The choice of line wing cutoff between 25 cm*-1*, 100 cm*-1*, Hybrid, 500 Lorentz and infinite widths has a negligible effect on the end emission spectra and produce similar results, which suggests that the H2O absorption features are formed at lower pressure pgas $\lesssim$ 1 bar atmospheric layers for this profile, where the pressure broadening effect is less sensitive to the cutoff length. We can therefore safely conclude that the input opacity treatment is not responsible for the small offsets seen in Sect. 5.3. However, an offset is seen for the results using the Stolker et al. (2017) opacities. Deeper absorption features are present, as well as a lower peak emission, similar to the benchmark comparison in Sect. 5.2. Possibly, the differences between the models for the Stolker et al. (2017) benchmarks stem from the combination of assuming an infinite cutoff for all gas species rather than just H2O considered here, such as that examined in Baudino et al. (2017).

The issue of the line wing cutoffs and pressure broadening for exoplanet and Brown Dwarf atmospheric opacities are well discussed in the literature (e.g. Sharp & Burrows, 2007; Amundsen et al., 2014; Grimm & Heng, 2015; Baudino et al., 2017) and extensively interrogated in Hedges & Madhusudhan (2016). Future detailed examination of the effects from the choice of line wing cutoff on forward modelling is warranted.

7.2 3D Geometric Effects

Since the aim of this paper is to compare a 3D radiative transfer model to contemporary 1D models, it is not surprising that differences would occur in the calculated spectra. Below we examine our methodologies, and suggest some differences that can lead to variations between our 3D approach and the 1D models.

The comparison to the benchmark transmission spectra produced by NEMESIS (Sect. 5.1) and Baudino et al. (2017) (Sect. 5.3) showed a systematic $\approx$5-15 ppm positive offset between the model output. In this work, a random starting impact parameter for the ray tracing is chosen on the transit annulus of the planet, in contrast to traditional 1D transmission codes where the ray tracing path is typically taken at the center of the radial cells. Assuming that an impact parameter width is evenly sampled, simply from spherical geometry, packets randomly starting nearer to the base of the impact parameter width pass through more atmosphere, and hence contribute a slightly higher optical depth in the transmission calculation. Since the transmission function is given by $T=\exp(-\tau)$, this biases the average transmission for a specific impact parameter width towards higher optical depths, resulting in a larger transit radius compared to the 1D approach.

In some of the emission spectra benchmarks in Sect. 5.3 a $\approx$5% additional flux output is seen in CMCRT compared to the 1D models. This offset is also wavelength dependent, occurring near the peak of the emission flux and the Rayleigh-Jeans tail. As a model that uses a 3D spherical grid, a fundamental difference to the 1D plane-parallel code is the geometry of the ray tracing through the atmosphere. Two main approximations break down for plane-parallel methods occurring near the limb of the planet:

The infinite-plane approximation becomes increasingly invalid as the emission angle approaches 90*∘*, as the sphericity of the planet becomes more important. 2. 2.

Near the limbs of the planet, radiation does not emerge from a vertical column of a $\tau$ $\sim$ 1 surface, but a combination of horizontal, low optical depth locations.

As a result of these approximations, the calculated slant path length of a ray to the top of the atmosphere becomes increasingly divergent as the emission angle approaches 90*∘*.

This effect gives rise to the well known differences in limb darkening behaviour for plane-parallel and spherical models in stellar atmospheres (e.g. Neilson, 2012). Since the T-p profiles used in Sect. 5.3 are highly isothermal in the upper atmosphere (Fig. 3), no limb-darkening effects are seen in the CMCRT output images or expected from the plane-parallel model either. This suggests that the likely origin of the offsets is directly the difference in path length between the spherical and plane-parallel models, with the 3D model allowing slightly more blackbody emission from the upper atmosphere isothermal regions to escape compared to the 1D approximations. However, since the offset is prominent in the Teff = 500 K and 1000 K benchmark tests and negligible in the higher temperature tests, it suggests that the magnitude of this effect is dependent on the modelled T-p profile.

Modelling emission near the limbs of tidally locked giant exoplanets may be an important consideration as it has been shown atmospheric jet structures can shift energy away from the sub-stellar point and towards the eastward limb (e.g. Heng & Showman, 2015). Extensive cloud structures confined to the westward limb are also theoretically expected to form (e.g. Parmentier et al., 2016; Roman & Rauscher, 2019), suggesting photon scattering effects near the limb of these exoplanets may also be important to factor into the radiative-transfer problem. This contrast between the hotter east, and colder and cloudy west terminator regions is well captured by the 3D MCRT model, and will be important for accurate determination of 3D feedback effects of radiative heating and cooling rates inside the atmosphere.

Overall, a rigorous and detailed quantitive comparison between 3D and 1D geometries and the manifestation of these effects in the observable spectra is beyond the scope of this study. However, such a study is warranted to ensure that radiative-transfer modelling using a variety of different methodologies are well calibrated.

8 Summary and Conclusions

We have expanded our 3D Monte Carlo radiative transfer code (Cloudy Monte Carlo Radiative Transfer: CMCRT), originally presented in Lee et al. (2017), to make use of k-distribution tables and the correlated-k approximation by statistically sampling the g-ordinance weights. We present a hybrid Monte Carlo and ray tracing method for the calculation of transmission spectra which include a multiple scattering component such as clouds and Rayleigh scattering. Our method reduces to an equivalent extinction, randomised transit chord algorithm in the absence of a scattering component, allowing the absorption and scattering contributions to the transmission spectra to be individually examined. We present a MCRT sampling method with k-distributions in emission spectra calculations. Our application highlights the synergy between the stochastic MCRT model, k-distribution properties and ray tracing methods.

We performed several benchmarking tests in CMCRT for transmission and emission spectra taken from the literature. Firstly, a direct comparison using identical inputs with the NEMESIS radiative-transfer suite (Irwin et al., 2008) compared highly favourably, with differences within 30 ppm for the transmission spectra and 10% for the emission spectra. Secondly, we benchmarked CMCRT to the emission spectra results in Stolker et al. (2017), producing consistent results when using the Stolker et al. (2017) opacity table directly, and offsets between the results when using the k-distribution opacities. Thirdly, we benchmarked CMCRT to the transmission and emission spectra for prescribed Guillot (2010) T-p profiles to output from the ATMO (Tremblin et al., 2015), Exo-REM (Baudino et al., 2015) and petitCODE (Mollière et al., 2015, 2017) as presented in Baudino et al. (2017). Our chemical equilibrium abundances compared well, with differences arising from thermal dissociation of molecules for the $T_{\rm eff}$ = 2000 K and 2500 K T-p profiles. Our transmission spectra results generally compare well to the 1D code to within 10s of ppm. The emission spectra results generally agree, except with a 5-10% continuum offset at infrared wavelengths in the CMCRT results for the $T_{\rm eff}$ = 500 K and 1000 K T-p profiles.

Lastly, we applied our new methodologies to post-process transmission spectra and emission spectra of the cloudy HD 189733b 3D GCM simulation from Lee et al. (2016). Our transmission spectra results with cloud produced too flat an optical slope and near-IR water features, suggesting that the modelled cloud particles in the GCM are too large and at too high altitude. Our emission spectra results are consistent with the near-IR HST (Crouzet et al., 2014; Barstow et al., 2014) and Spitzer (Knutson et al., 2012, and references within) photometric data. Comparing to the Todorov et al. (2014) Spitzer IRS data suggests that the peak emission on the dayside of the GCM may be too high, indicating a cooler dayside photospheric temperature for the real object. The presence of an emission feature at the cloud particles stretching and bending modes may be suggestive of an atmospheric temperature inversion where clouds are present, and also reveal the composition of the particles.

We examined the effect of the far line wing cutoff parameter on the emission spectra results, finding negligible differences in our H2O tests. A possible candidate for the differences found between the Stolker et al. (2017) and the CMCRT results is the combination of gas species assuming an infinite cutoff as discussed in Baudino et al. (2017).

We suggest a 3D spherical geometry effect from randomly sampling impact parameters in CMCRT, which leads to biasing towards a lower transmission through each impact parameter layer, resulting in a slightly higher (ppm level) planetary radius compared to the 1D benchmarks. We suggest it is likely due to the differences in path lengths between 3D spherical and 1D plane-parallel models at the limb of the simulated planet, a small additional blackbody flux component from the upper atmospheric regions is present in the 3D model compared to the 1D models, the magnitude of which is dependent on the simulated T-p profile.

Despite differences between individual comparisons of output, our benchmarking efforts as a whole suggest our methods produce highly consistent transmission and emission spectra results of comparable quality to the 1D model outputs. We also conclude from the GCM post-processing results that CMCRT is extremely well suited for accurate simulation of the radiative environment of the complex, 3D inhomogeneous nature of exoplanet atmospheres. The importance of considering a 3D atmospheric structure on the interpretation of observational and model results is beginning to be explored in detail (e.g. Feng et al., 2016; Blecic et al., 2017). As observational data on these objects increases in quality and quantity, 1D radiative-transfer techniques may quickly become intractable to modelling the 3D spatial complexity of the atmosphere and 3D radiative-transfer techniques more appropriate for the interpretation of the observational data.

This study provides a cautionary tale when comparing between forward model outputs, it can be challenging to diagnose the origin of any systematic difference between models, be it 3D geometrical effects or the individual treatment and parameters used by different groups when calculating input opacities. Our results also hopefully guide the future of cloud modelling in GCMs, consistently comparing cloudy GCM output to observational data is invaluable for further refinement of GCM modelling techniques and development and testing of cloud formation theories.

Acknowledgements

G.K.H. Lee acknowledges support from the University of Oxford and CSH Bern through the Bernoulli fellowship and support from the European community through the ERC advanced grant EXOCONDENSE (PI: R.T. Pierrehumbert). The participants of the OWL 2018 summer program are thanked for numerous discussions and encouragement on the development of the model. J. Taylor is a Penrose Scholar and would like to thank the Oxford Physics Endowment for Graduates (OXPEG) for funding this research. J-L Baudino acknowledges the support of the UK Science and Technology Facilities Council. T. Stolker and M. Min are thanked for providing data on the ARTES opacities and model output. V. Parmentier and R.T. Pierrehumbert are thanked for discussions and suggestions on the manuscript content. Most plots were produced using the community open-source Python packages Matplotlib (Hunter, 2007), SciPy (Jones et al., 2001), and AstroPy (The Astropy Collaboration et al., 2018). The authors’ local HPC support at Oxford and CSH Bern is highly acknowledged.