3D Radiative Transfer for Exoplanet Atmospheres. gCMCRT: A GPU-accelerated MCRT Code

The volume of observational data of exoplanet atmospheres from telescopes continues to grow exponentially, uncovering a wide variety of atmospheric properties across the exoplanet population. With the upcoming launch of JWST (Bean et al. 2018) and the Ariel mission (Tinetti et al. 2016), precise spectral data across the infrared wavelength regime will become available for the community to analyze. Kepler (Howell et al. 2014), TESS (Ricker et al. 2014), and CHEOPS (Broeg et al. 2013) have provided optical-band photometric data of exoplanet phase curves, allowing examination of their albedo and thermal properties (e.g., Demory et al. 2013; Wong et al. 2020; Morris et al. 2021). The unique 3D nature of each exoplanet atmosphere manifests itself in the observational properties of that specific object. Therefore, a 3D understanding of their atmospheres is key to a holistic investigation of their atmospheric structure and composition. For phase-curve observations, different hemispheres of the planet come in and out of view of the observational direction, leading to an intrinsically 3D mapping of the atmosphere. For the dayside of atmospheres that are highly irradiated and puffed up, for example, ultrahot Jupiters (UHJs), the radiation from the dayside can be blended with the cold nightside emission when both the dayside and nightside hemispheres are in view, for example, when taking phase-curve observations (e.g., Irwin et al. 2020; Taylor et al. 2021). These close-in planets are also expected to be tidally locked to their host star, showing a permanent dayside and nightside hemisphere. In addition, as a result of the 3D dynamical properties of these UHJ planets and their short dayside radiative timescales, their day-to-night energy transport can be weak (e.g., Parmentier et al. 2018). Accounting for all these effects accurately requires 3D radiative transfer (RT) modeling techniques, able to produce spectra taking into account the contribution from the entire 3D atmosphere to the end result.

In the past decade, the number and precision of high-resolution observations have greatly increased, revealing the rich chemical composition of hot Jupiter atmospheres (e.g., Seidel et al. 2019; Ehrenreich et al. 2020; Merritt et al. 2020; Pino et al. 2020). Current high-resolution spectral instruments, for example, ESPRESSO (Pepe et al. 2021), CARMENES (Quirrenbach et al. 2016), GIANO-B (Oliva et al. 2012), and CRIRES+ (Follert et al. 2014), provide detailed mapping of the chemical structure and wind profiles in the atmosphere. The fidelity and wavelength range of such high-resolution data will also increase as the next generation of telescopes and instruments come online in the next decade, such as HIRES on the EELT (Marconi et al. 2016). The number of wavelength points for these high-resolution studies is large. For example, there are approximately 170,000 wavelength points in the CRIRES+ wavelength range (0.95 μm < λ < 5.2 μm) for a resolution of R = 100,000. Performing RT for 3D models with this many wavelength points is a challenging prospect, not withstanding the fact that modeling at high resolution is typically performed at many times the resolution of the spectrograph. Should phase information be required, this RT modeling needs to be repeated for each phase of interest. For modelers, there is a need for fast, accurate, and flexible 3D model post-processing capabilities that are able to meet the challenges of modeling planets at high resolution.

As the number of characterized planets increases in both low and high resolution, with thousands of exoplanets scheduled to be characterized in the coming decades, interpreting physical trends from the observational data may become possible. However, this requires large parameter grids of models to be run, which can be a time-consuming and expensive undertaking. Several modeling groups have performed parameter investigations using global circulation models (GCMs), from examining effects of equilibrium temperature, atmospheric drag, and chemistry (Komacek et al. 2017, 2019; Baeyens et al. 2021) to grids of models used to examine cloud properties across the hot Jupiter equilibrium regime (Parmentier et al. 2016, 2021; Roman et al. 2021). Post-processing these large grids is a significant project in itself, typically requiring processing simulations manually one by one. Therefore, typically modelers use available data from the GCM model itself, such as the temperature structures and outgoing longwave radiation data, to make large-scale diagnostic prescriptions. For future large GCM grid projects, there is a need for a uniform and fast spectral 3D RT processing method, able to process the raw GCM output formats and perform the RT modeling in a simple and automated manner.

Recently, several studies have investigated biases that occur when the 3D structure of the atmosphere is not taken into account. In emission, Feng et al. (2016) and Taylor et al. (2021) investigate the biases that can occur when retrieving emission spectra assuming a single T–p profile across a hemisphere that is inhomogeneous in temperature. They find in some cases that assuming a single global T–p profile can produce spurious molecular detections from the blending of cooler and hotter profiles of the atmosphere. Blecic et al. (2017) retrieved 1D T–p profiles of 3D GCM output from Dobbs-Dixon & Agol (2013). They found that it was challenging to represent the 3D T–p structure dayside well using 1D retrieved T–p profiles. Irwin et al. (2020) and Feng et al. (2020) present 2.5D retrieval methods that are able to take into account the variable T–p structures of the atmosphere, increasing accuracy when retrieving phase-curve data. Dobbs-Dixon & Cowan (2017) show how the wavelength-dependent photosphere changes across the planet owing to changes in the temperature and chemical composition, showing a complex contribution function structure that varies across the globe. Fortney et al. (2019) investigate this photospheric radius effect, suggesting that not taking into account the 3D emitting area of the planet can lead to errors of 5% or more in the flux.

In transmission, Caldas et al. (2019) and Pluriel et al. (2020) investigated biases that occur when not assuming 3D transmission geometry, showing significant differences between the full 3D model and a mean profile and that retrievals may produce biased results if the 3D chemical and temperature structure of the atmosphere is not taken into account. Line & Parmentier (2016) show that nonuniform, patchy cloud structures between the east and west terminators can be degenerate with high molecular weight. MacDonald et al. (2020) and Lacy & Burrows (2020) have also shown that chemical gradients in the 3D structure across both east and west limbs can also produce biased results when retrieved in 1D.

To date, modeling high-resolution spectra in 3D has mostly stemmed from the Miller-Ricci Kempton & Rauscher (2012) model, with the exception of the Showman et al. (2013) study. Flowers et al. (2019) post-process 3D hot Jupiter GCMs to investigate the high-resolution transmission spectra and their cross-correlation properties. Beltz et al. (2021) show that post-processing using the 3D structure of GCM models improves the detection significance of molecules at high resolution in emission. Both of these studies used output from a Rauscher & Menou (2012) GCM model. Savel et al. (2021) use MITgcm models of WASP-76b from May et al. (2021) to produce high-resolution spectra, finding that high-altitude cloud particles can reproduce well the phase variations observed by Ehrenreich et al. (2020). Wardenier et al. (2021) use SPARC/MITgcm (Showman et al. 2009) GCM output of a WASP-76b model, finding that Fe depletion at the limbs of the planet can also provide a reasonable explanation for the Ehrenreich et al. (2020) data.

From the above studies, it is highly warranted to develop the accurate and flexible 3D RT techniques required to tackle the challenges associated with the upcoming high-fidelity observational data on exoplanet atmospheres. These requirements make the Monte Carlo Radiative Transfer (MCRT) method a highly suitable methodology to meet these demands, as it is able to directly take into account the 3D geometry of the atmosphere, including multiple scattering effects and other photon microphysics. In transmission, there is a need for a model able to take into account both the chemical and cloud gradients in a 3D, inhomogeneous manner, weighting correctly the contribution of each transmission limb to the transmission spectra. In emission, a model must now take into account the inhomogeneous nature of the 3D temperature structure and correctly weight each emitting area of the exoplanet.

MCRT is a standard 3D RT method in the astrophysical community (Steinacker et al. 2013; Noebauer & Sim 2019), but it has seen limited use in planetary contexts (e.g., Hood et al. 2008; de Kok & Stam 2012; Garcia Munoz & Isaak 2015; Robinson 2017; Stolker et al. 2017). MCRT methods have also been used for simulations of RT through Earth clouds (e.g., Mayer 2009). In astrophysics, MCRT models have been used to post-process hydrodynamical models to produce synthetic spectra and compare to observation data (e.g., Robitaille 2011). Additionally, directly using MCRT as the RT solver inside the hydrodynamical model itself has been performed, for example, in the protoplanetary disk modeling of Harries et al. (2019). 3D GCM model output is typically transformed into latitude–longitude grids, making 3D spherical geometry a natural choice for the MCRT model. Of primary interest is producing transmission, emission, and albedo spectra, as well as phase curves. Enabling the calculation of 3D contribution functions is also important to establish what sections of the planet are contributing to the observable features.

In this study, we develop gCMCRT, an open-source GPU version of the Lee et al. (2017, 2019) CMCRT model ⁸ for general use in planetary atmosphere RT. CMCRT has been used to post-process GCM simulations from Lee et al. (2016, 2020) and Lee et al. (2021), showing its usefulness in 3D contexts. In Wardenier et al. (2021), a high wavelength resolution version was developed, including the Doppler shift of the local opacity from rotation and winds, to investigate the cross-correlation signal of Fe in transmission of WASP-76b (Ehrenreich et al. 2020; Kesseli & Snellen 2021), demonstrating the importance of considering the 3D structure of the atmosphere for high-resolution interpretation as well. This used SPARC/MITgcm model output (Showman et al. 2009) as input to CMCRT.

By applying GPU technology to CMCRT, it is now easier and faster to perform 3D modeling without the need of high-performance CPU environments. Due to the high level of parallelism of MCRT methods, the GPU offers a speedup of 3–10× over the older CPU version. This current study was performed solely on a desktop computer with an Nvidia RTX 3080 GPU card.

In this paper, we develop several standard modes for gCMCRT:

• 1.  
  Albedo spectra
• 2.  
  Transmission spectra
• 3.  
  Emission spectra
• 4.  
  Hi-res transmission spectra
• 5.  
  Hi-res emission spectra.

The albedo and emission modes are able to model phase curves by performing the experiment across multiple viewing angles and at different planetary phases. We outline the physics and mechanisms used in gCMCRT to calculate spectral properties of model atmospheres. We post-process hot Jupiter GCMs using standard output from the Rauscher & Menou (2010, 2012, 2013) model, SPARC/MITgcm (Showman et al. 2009), UK Met Office UM (Mayne et al. 2014), ExoRad (Carone et al. 2020), THOR (Mendonça et al. 2016; Deitrick et al. 2020), and Exo-FMS (Lee et al. 2021). Each model is used to showcase a mode of gCMCRT.

In Sections 2, 3, and 4 we detail how the MCRT techniques and physics are performed in gCMCRT. Section 5 describes how observable quantities are derived from the gCMCRT output. Section 6 shows how 3D contribution functions are calculated for the albedo and emission modes. Section 7 details the opacity sources that gCMCRT can use. Section 8 details how chemical abundances are implemented in gCMCRT. Section 9 shows the results of our 3D albedo spectra mode post-processing. Section 10 shows the results of the emission spectra post-processing. Section 11 shows the results of the transmission spectra post-processing. Section 12 shows the results of the high-resolution emission spectra post-processing. Section 13 shows the results of the high-resolution transmission spectra post-processing. Sections 14 and 15 present the discussion and conclusion of our results, respectively.

gCMCRT is an MCRT code that simulates the path of photon packets through a planetary atmosphere in 3D. It is a hybrid MCRT and ray-tracing code, using the "peeloff" ray-tracing method (Section 4.6) to produce images and spectra of the simulated planet. Both the basic MCRT and ray-tracing capabilities can be used separately or together, depending on the desired output. Basic MCRT methods have been extensively detailed in many sources (e.g., Dupree & Fraley 2002; Whitney 2011; Noebauer & Sim 2019). gCMCRT primarily uses a spherical geometry grid, with packets moving and interacting inside the 3D computational domain. As such, it is a subgrid method, meaning that packets evolve through 3D space and the ray-tracing is performed within the spherical geometry, not requiring special placement or selection of vertical layers.

gCMCRT is originally based on the model published by Hood et al. (2008), which has its origins in MCRT codes used to investigate RT and photoionization for astrophysical applications (e.g., Wood et al. 2004). This model was expanded and specialized for exoplanetary atmosphere science in Lee et al. (2017) and benchmarked to contemporary 1D RT codes with correlated-k opacities in Lee et al. (2019). In this work, in addition to adapting to GPU architecture, we further expand the capabilities of the model, making it suitable for a wide variety of RT situations found in planetary science.

The user interacts with the MCRT routines through a custom "main experiment" routine, which the user can then use to call the core routines to perform a series of operations for the photon packets for their desired setup. We include standard setups for each mode to provide a baseline usability of the model.

MCRT is generally described as "embarrassingly parallel," with the program able to be relatively simply parallelized for multiple processors. An advantage of utilizing GPU cards is the large number (thousands) of cores on each GPU card compared to CPUs (tens). This makes them ideal for simulating MCRT models where a large number of independent, simple calculations are required to be performed. This allows the model to be run on desktop machines with GPU cards cheaply, whereas the CPU version required specialized HPC equipment to function efficiently. GPU nodes have also become a popular addition to compute servers in recent years, providing magnitudes of computing scalability for models able to take advantage of the large number of cores.

gCMCRT is written in CUDA FORTRAN, chosen because it is generally a more familiar coding language to scientists in the astrophysics and (exo)planetary science fields, allowing the user a simpler designing of experiments to interface with the core routines. CUDA FORTRAN retains the FORTRAN syntax and program layout, with additional special syntax that enables communication with the GPU architecture. This also allows future development of gCMCRT to take advantage of other well-tested MCRT methods such as photoionization (Wood et al. 2013).

We have written gCMCRT in a highly modular way, with the core routines modifying packet properties and performing the ray-tracing calculations. Grid and opacity data are stored in global GPU modules, able to be accessed by routines that require them. The host (CPU) side of the code reads and calculates the opacity structure of the atmosphere, which is then passed to the GPU memory. Image data are also stored in global GPU memory, updated by each packet as it contributes to the image. After the iteration is finished, the GPU data are given back to the host for output.

In testing, the main bottleneck was found to be the reading of 3D opacity data from the hard disk by the CPU. We have optimized this section of code as best as possible, such as using asynchronous I/O and loading the maximum amount of data into memory as possible for each read statement.

We have specifically avoided the use of additional packages inside the FORTRAN code (such as HDF5) to simplify the usability of the model. Instead, output is produced in FORTRAN binary files, which can then be converted to more flexible storage formats later.

The main methods used in gCMCRT are already detailed in Lee et al. (2017) and Lee et al. (2019). However, in the following sections we provide a brief review of the main processes and describe new capabilities included in gCMCRT.

gCMCRT includes the ability to treat multiple or single scattering using a variety of different phase functions. Due to the flexibility of the method, mixing different scattering phase functions for varying scattering sources can also be performed.

3.1. Isotropic Scattering

For isotropic scattering the phase function, $P(\cos \theta )$, is given by

$$\begin{eqnarray}&&{P}_{\mathrm{Iso}}(\cos \theta )=1.\end{eqnarray} \tag{ 1 }$$

Sampling a cosine angle for isotropic scattering is simply given by

$$\begin{eqnarray}&&\cos \theta =2\zeta -1,\end{eqnarray} \tag{ 2 }$$

where ζ is a uniformly sampled random number between 0 and 1. This equally samples a direction across 4π sr.

3.2. Lambertian Surface Scattering

For scattering of photon packets off a surface, a common and simple phase function to use is the Lambertian phase function

$$\begin{eqnarray}&&{P}_{\mathrm{Lam}}(\cos \theta )=\cos \theta ,\end{eqnarray} \tag{ 3 }$$

which equally samples a hemisphere of π sr. Sampling a cosine angle for Lambertian scattering is given by

$$\begin{eqnarray}&&\cos \theta =\sqrt{\zeta }.\end{eqnarray} \tag{ 4 }$$

3.3. Rayleigh Scattering

For scattering off gas particles or small size parameter cloud particles, the Rayleigh scattering phase function is used, given by

$$\begin{eqnarray}&&{P}_{\mathrm{Ray}}(\cos \theta )=\displaystyle \frac{3}{4}(1+{\cos }^{2}\theta ),\end{eqnarray} \tag{ 5 }$$

and a cosine angle for unpolarized Rayleigh scattering is sampled through (e.g., Mayer 2009; Frisvad 2011)

$$\begin{eqnarray}&&\cos \theta =\sqrt[3]{-q+\sqrt{1+{q}^{2}}},q=4\zeta -1.\end{eqnarray} \tag{ 6 }$$

3.4. Aerosol Scattering

A key characteristic of planetary atmospheres is the presence of aerosol material, either photochemically produced hazes or cloud particles formed by condensation. Accurate treatment of light scattering by aerosols is a vital part of RT, and gCMCRT includes several options to treat aerosol scattering.

The Henyey–Greenstein (HG) phase function (Henyey & Greenstein 1941) has been extensively used as a one-parameter fit to the Mie scattering phase function. The HG function is given by

$$\begin{eqnarray}&&{P}_{\mathrm{HG}}(\cos \theta )=\displaystyle \frac{1-{g}^{2}}{{\left(1+{g}^{2}-2g\cos \theta \right)}^{3/2}},\end{eqnarray} \tag{ 7 }$$

where $g=\left\langle \cos \theta \right\rangle$ is the asymmetry factor defined as the mean cosine scattering angle. Sampling a cosine angle from the HG function is given by (e.g., Sharma 2015)

$$\begin{eqnarray}&&\cos \theta =\displaystyle \frac{1}{2g}\left[1+{g}^{2}-{\left(\displaystyle \frac{1-{g}^{2}}{1-g+2g\zeta }\right)}^{2}\right].\end{eqnarray} \tag{ 8 }$$

One of the main limitations of the HG function is that as g approaches zero, the phase function becomes isotropic rather than Rayleigh. This generally occurs for IR wavelengths and for small particles at optical wavelengths. Such particles make a significant proportion of the cloud size distributions in exosolar atmospheres (e.g., Powell et al. 2018). To address this, Draine (2003) proposes a hybrid HG and Rayleigh phase function with the form

$$\begin{eqnarray}&&{P}_{{\rm{D}}03}(\cos \theta )=\displaystyle \frac{1-{G}^{2}}{1+\alpha (1+2{G}^{2})/3}\displaystyle \frac{1+\alpha {\cos }^{2}\theta }{{\left(1+{G}^{2}-2G\cos \theta \right)}^{3/2}},\end{eqnarray} \tag{ 9 }$$

where G(g) is now a function of the asymmetry parameter g, assuming that α is given (Draine 2003; Sharma 2015). When α = 1, this distribution is equivalent to the Cornette & Shanks (1992) function. We note that in Draine (2003) a better fit to the Mie scattering phase function was produced for interstellar dust when α < 1; therefore, we assume α = 0.5 as a default in CMCRT. Sampling a cosine angle from Equation (9) is performed following the Gibbs sampling approach found in Zhang (2019).

A variant of the HG function is the two-term HG function (TTHG), a linear combination of two HG functions, given by

$$\begin{eqnarray}&&{P}_{\mathrm{TTHG}}(\cos \theta )=\alpha {P}_{\mathrm{HG}}({g}_{1})+(1-\alpha ){P}_{\mathrm{HG}}({g}_{2}),\end{eqnarray} \tag{ 10 }$$

with α representing the forward-scattered fraction, g₁ ≥ 0 the forward-scattered asymmetry parameter, and g₂ ≤ 0 the backward-scattered asymmetry parameter. This attempts to better capture the backward scattering fraction produced by Mie calculations, especially important at optical wavelengths. Sampling a cosine angle for the TTHG function is performed by sampling a random number. Should ζ < α, the new direction is sampled using Equation (8) with g = g₁; otherwise, Equation (8) is sampled with g = g₂ (Pfeiffer & Chapman 2008).

For the TTHG function, there is no known general solution to find the parameters α, g₁, and g₂ from the asymmetry factor g directly, though fitting to individual aerosol materials using higher moments of the TTHG function has been performed (e.g., Kattawar 1975). As the default in CMCRT, we use the simple relationship used in Cahoy et al. (2010) and Batalha et al. (2019), which has been found useful for reproducing Jupiter reflectance spectra; these are given as g₁ = g, g₂ = −g/2, and $\alpha =1-{g}_{2}^{2}$.

An interesting consequence of the MCRT method is that the weighted asymmetry factor for the mixture of gas and cloud particles does not need to be calculated. Instead, the probability of a scattering event off a gas particle, P_gas, is sampled directly as

$$\begin{eqnarray}&&{P}_{\mathrm{gas}}=\displaystyle \frac{{\kappa }_{\mathrm{Ray}}}{{\kappa }_{\mathrm{Ray}}+{\kappa }_{\mathrm{sca},\ \mathrm{cld}}},\end{eqnarray} \tag{ 11 }$$

where κ_Ray [cm² g⁻¹] is the Rayleigh scattering opacity and κ_{sca, cld} [cm² g⁻¹] is the scattering opacity of the cloud. A random number can then be drawn every scattering event to determine the appropriate scattering phase function to sample.

The act of including variance reduction schemes, also known as "biasing," attempts to address both the signal properties and variance issues present in a default MCRT scheme. For 3D applications, some form of variance reduction is usually required beyond simple MCRT setups and is key to retaining high accuracy and reducing noise in the model. Below we briefly outline the variance reduction techniques currently implemented as options in gCMCRT.

4.1. Continuous Absorption

In the default, nonbiased MCRT methodology, the weight of the packet (representing the fraction of energy carried by the packet) is usually kept constant as it travels through the simulation, until it interacts with the medium in some way. In the continuous absorption method, however, the weight of the packet degrades smoothly as it travels through the simulation, given by (e.g., Noebauer & Sim 2019)

$$\begin{eqnarray}&&{W}_{\mathrm{ph},\ \mathrm{new}}={W}_{\mathrm{ph}}{e}^{-{\tau }_{\mathrm{abs}}},\end{eqnarray} \tag{ 12 }$$

where W_{ph, new} is the new weight of the packet, W_ph the original weight, and τ_abs the optical depth of the absorption opacity. This smooths the weight reduction of the packets as they travel through the simulation and is mostly useful for when the atmosphere is highly absorbing.

4.2. Forced Interaction

4.3. Emission Biasing

In a nonbiased MCRT scheme, the number of packets simulated in emission for each cell, N_i , of the grid is given proportionally to the total number of packets simulated, N_tot, from the contribution of the cell, L_i [erg s⁻¹], to the total luminosity of all cells together, L_tot [erg s⁻¹],

$$\begin{eqnarray}&&{N}_{\mathrm{ph},{\rm{i}}}={N}_{\mathrm{tot}}\displaystyle \frac{{L}_{i}}{{L}_{\mathrm{tot}}}.\end{eqnarray} \tag{ 13 }$$

However, for cells that contribute negligibly to the total luminosity (L_i ≪ L_tot), very few packets are simulated originating from the cell, leading to increased noise and not properly accounting for their contribution to the end spectrum. To counter this, we apply the hybrid biasing scheme from Baes et al. (2016), which combines a uniform-in-cell biased function with the nonbiased function. We have found that high biasing is generally required to adequately capture the contribution from cooler regions of the planetary atmosphere and is key in producing accurate phase-curve spectra for the large temperature contrasts present in hot Jupiter atmospheres.

4.4. Survival Biasing

4.5. Russian Roulette

4.6. Next Event Estimation

The next event estimation method, also known as the peeloff method, combines the regular MCRT model with geometric ray-tracing applied at every scattering or emitting location of the packet (Yusef-Zadeh et al. 1984; Wood & Reynolds 1999). This greatly reduces sampling noise and is the main tool that provides high-fidelity output for the simulation. The fraction, f, of a packet's luminosity escaping toward the chosen observational direction during a scattering event is given as

$$\begin{eqnarray}&&f={W}_{\mathrm{ph}}\omega {\rm{\Theta }}({\theta }_{\mathrm{obs}})\exp (-{\tau }_{\mathrm{obs}}),\end{eqnarray} \tag{ 14 }$$

where Θ(θ_obs) is the normalized scattering phase function toward the observational direction, ω the single scattering albedo, and τ_obs the total extinction optical depth toward the observational direction.

For emission spectra, should multiple scattering not be required, an initial peeloff from the emission location is enough to produce the spectrum. To increase the efficiency of the calculation, a limiting vertical optical depth can be set to avoid simulating packets from very high optical depths that contribute negligibly to the output spectrum (here typically τ_vert = 30). This focuses the computational effort to produce low-noise results from the emission in important photospheric regions.

4.7.  g-ordinance Biasing

When modeling using correlated-k opacities, the g-ordinance is required to be randomly sampled from the cumulative weight distribution of the g-ordinances (Lee et al. 2019). For an accurate result, the opacity distribution should be adequately sampled across the range of weights. However, for regions of low packet counts, this distribution may not be as evenly sampled as desired.

For unbiased g-ordinance sampling the probability, p(g), is given directly by the weight of the g-ordinance, w_g ,

$$\begin{eqnarray}&&p(g)={w}_{g}.\end{eqnarray} \tag{ 15 }$$

For direct sampling of the k-distribution weights, we can construct a simple biased probability distribution function, q(g), to equally sample all g-ordinances

$$\begin{eqnarray}&&q(g)=\displaystyle \frac{1}{{N}_{g}},\end{eqnarray} \tag{ 16 }$$

where N_g is the total number of g-ordinance points, with the weighting factor

$$\begin{eqnarray}&&{W}_{\mathrm{ph}}(g)=\displaystyle \frac{p(g)}{q(g)}={N}_{g}{w}_{g}.\end{eqnarray} \tag{ 17 }$$

For sampling g-ordinances in emission, the unbiased g-ordinance sampling is given by

$$\begin{eqnarray}&&p(g)=\displaystyle \frac{{w}_{g}{L}_{g}}{{\sum }_{g}{w}_{g}{L}_{g}},\end{eqnarray} \tag{ 18 }$$

where L_g [erg s⁻¹] is the luminosity of that g-ordinance. For the biased sampling, the same uniform probability as Equation (16) can be assumed. The weighting factor is then

$$\begin{eqnarray}&&{W}_{\mathrm{ph}}(g)=\displaystyle \frac{p(g)}{q(g)}=\displaystyle \frac{{N}_{g}{w}_{g}{L}_{g}}{{\sum }_{g}{w}_{g}{L}_{g}}.\end{eqnarray} \tag{ 19 }$$

We have also included a composite biasing method for sampling the g-ordinances. In testing we have found this biasing to be vital when using k-tables that use split quadrature zones, for example, the 8+8 points used in petitRADTRANS (Mollière et al. 2019) and the 4+4 points used in SPARC/MITgcm (Kataria et al. 2013). Due to the low weights for the high split quadrature values, the nonbiased scheme will rarely sample these g-ordinance values, resulting in too low an opacity distribution being sampled. The biasing scheme fixes this issue in a simple manner and ensures a more even sampling of the opacity distribution. We have found that uniform sampled g-ordinance values perform similarly with or without this biasing. However, we recommend this biasing for all cases anyway to ensure that the low weighted g-ordinances are adequately sampled.

In this section, we detail how albedo, transmission, and emission spectra are calculated from the gCMCRT output.

5.1. Transmission Spectra

gCMCRT has the ability to calculate fully 3D transmission spectra including the effects of multiple scattering. The transmission spectra equation is given by (e.g., Dobbs-Dixon & Agol 2013; Robinson 2017)

$$\begin{eqnarray}&&{\left(\displaystyle \frac{{R}_{p,\lambda }}{{R}_{\star }}\right)}^{2}=\displaystyle \frac{1}{{R}_{\star }^{2}}\left({R}_{p,0}^{2}+2{\int }_{{R}_{p,0}}^{\infty }[1-{ \mathcal T }(b)]{bdb}\right),\end{eqnarray} \tag{ 20 }$$

where R_p,λ [m] is the wavelength-dependent radius of the planet, R_⋆ [m] the radius of the host star, R_p,0 [m] the bulk planetary disk radius, ${ \mathcal T }$ the transmission function, and b [m] the impact parameter. Formally the upper limit for the integral in Equation (20) is ∞ . This is replaced by the top-of-atmosphere radius, R_p,TOA [m], as per the simulation boundaries to facilitate numerical calculations.

Following the principles of integration through independent sampling, the result of the integral in Equation (20), I_p , is approximated by simulating a suitably large number of N_ph photon packets that sample the integral function

$$\begin{eqnarray}&&\langle {I}_{p}\rangle =\displaystyle \frac{({R}_{p,\mathrm{TOA}}-{R}_{p,0})}{{N}_{\mathrm{ph}}}\sum _{i}^{{N}_{\mathrm{ph}}}{W}_{\mathrm{ph}}[1-{e}^{-{\tau }_{i}}]{b}_{i},\end{eqnarray} \tag{ 21 }$$

where now each packet contributes its transmission for a randomly sampled impact parameter. This method avoids the slight geometric biasing found in the original method in Lee et al. (2019), where packets were binned in impact parameter.

An initial peeloff (ray-tracing step) is performed at the packet's initial transmission location, giving a baseline signal for the absorption spectra. The packet is then free to multiply scatter through the atmosphere, each time contributing to the transmission through the atmosphere through the next event estimation method. This weights the transmission with the scattering probability toward the observation direction, so packets can contribute to the transmission at different impact parameters from their initial position. Packets are given a random starting location on the transmission annulus of the simulation, ensuring that ray-tracing from the packet's ensuing transit chord will most likely not occult the planetary surface boundary, which would not contribute to the signal.

Should multiple scattering not be required, it can be switched off, which reduces the model to a ray-tracer scheme with random transit chord sampler. This significantly increases the efficiency of the model and is useful for when it is known that multiple scattering will not be an important contribution to the transmission spectrum, which is typically the case for infrared wavelength regimes (e.g., Hubbard et al. 2001).

5.1.1. Limb Darkening

An important consideration for transmission spectra calculations is the change in the stellar disk intensity that passes through the atmosphere as the planet transits across the stellar disk. For our transmission spectra modes we include the capability to take the limb darkening of the star into account.

To calculate this, the sampled transmission chord is mapped using spherical geometry directly to a longitude and latitude, and therefore a zenith angle, μ, onto the spherical stellar surface. This cosine angle is then converted to the limb-darkening fraction, I_ph(μ)/I(0), the ratio of the intensity at μ to the zenith point, using a limb-darkening law, for example, the quadratic law (Kopal 1950)

$$\begin{eqnarray}&&\displaystyle \frac{{I}_{\mathrm{ph}}(\mu )}{I(0)}=1-{c}_{1}(1-\mu )-{c}_{2}{\left(1-\mu \right)}^{2},\end{eqnarray} \tag{ 22 }$$

where c₁ and c₂ are fitting coefficients. We include all limb-darkening laws detailed in J. Southworth's website ⁹ with user-supplied limb-darkening coefficients.

Each packet then has a unique limb-darkening fraction value, which is then used to reduce the weight of the packet through

$$\begin{eqnarray}&&{W}_{\mathrm{ph},\mathrm{new}}={W}_{\mathrm{ph}}\displaystyle \frac{{I}_{\mathrm{ph}}(\mu )}{I(0)},\end{eqnarray} \tag{ 23 }$$

which is then the correctly weighted packet's contribution to the transmission signal.

5.2. Scattered Incident Light

The geometric albedo, A_g , is defined as the ratio of incident to reflected flux at zero phase (e.g., Marley et al. 1999; Heng et al. 2021)

$$\begin{eqnarray}&&{A}_{g}=\displaystyle \frac{j(0)}{\pi {F}_{\mathrm{inc}}},\end{eqnarray} \tag{ 24 }$$

where j(0) is the outgoing flux at 0 phase and F_inc the incident flux. In gCMCRT, the reflected fraction of the energy escaping toward the observable direction is directly calculated, with the incident flux normalized to 1. For a geometric albedo calculation, therefore, the results for a reflection calculation at 0 phase are simply multiplied by π to account for the definition of A_g . The geometric albedo can also be given as a function of wavelength.

For scattered light phase curves, the viewing angle of the detector can be changed for the desired phases. The output is then scaled to the reflected fraction at zero phase, essentially deriving the normalized classical planetary phase function for the simulated planet, given by

$$\begin{eqnarray}&&{\rm{\Phi }}(\phi )=\displaystyle \frac{j(\phi )}{j(0)}.\end{eqnarray} \tag{ 25 }$$

The planet-to-star flux ratio is then given by the well-known formula

$$\begin{eqnarray}&&\frac{{F}_{p}}{{F}_{\star }}(\phi )={\left(\frac{{R}_{p}}{a}\right)}^{2}{A}_{g}{\rm{\Phi }}(\phi ),\end{eqnarray} \tag{ 26 }$$

where R _p [cm] is the planetary radius and a [cm] the semimajor axis.

It is useful to note that in this case Φ(ϕ) can be >1 dependent on the variety of scattering components in the atmosphere, for example, a bright westward cloud, allowing nonsymmetrical phase curves to be modeled.

5.3. Emitted Planetary Light

To calculate emission spectra, each photon packet is emitted from within the atmospheric volume, and its contribution to the outgoing energy toward a certain viewing angle is tracked. For cloud-free simulations, scattering is generally negligible and can be switched off for faster results, especially when focusing on infrared wavelengths.

The luminosity of an individual cell i, L_λ [erg s⁻¹ cm⁻¹], for a given wavelength λ is given by

$$\begin{eqnarray}&&{L}_{\lambda ,i}=4\pi {V}_{i}{\rho }_{i}{\kappa }_{\mathrm{abs},\ i}{B}_{\lambda }({T}_{i}),\end{eqnarray} \tag{ 27 }$$

where V_i [cm³] is the volume of the cell, ρ_i [g cm⁻³] the density of the cell, κ_{abs, i} [cm²g⁻¹] the absorption opacity of the cell (including any cloud absorption opacity), and B_λ (T_i )[erg s⁻¹ cm⁻² cm⁻¹ sr⁻¹] the Planck function at temperature T_i .

The luminosity escaping toward the observable direction, L _p [erg s⁻¹ cm⁻¹] , is then

$$\begin{eqnarray}&&{L}_{p}=\frac{{L}_{\mathrm{tot}}}{{N}_{\mathrm{ph}}}\displaystyle \sum _{i}\displaystyle \sum _{n}^{{N}_{\mathrm{ph}}}f(\phi ,\theta ),\end{eqnarray} \tag{ 28 }$$

where L_tot = ∑L_λ .

The specific intensity emanating from the planet toward a certain viewing angle, I_p [erg s⁻¹ cm⁻² cm⁻¹], is then given by

$$\begin{eqnarray}&&{I}_{p}=\frac{{L}_{p}}{\pi {\left({R}_{p,0}+{z}_{\mathrm{TOA}}\right)}^{2}},\end{eqnarray} \tag{ 29 }$$

where z_TOA [cm] is the top-of-atmosphere altitude. In this case the denominator of this equation is always an area of a circle, as the output of gCMCRT is the luminosity escaping from a spherical hemisphere of the planet. The flux from the planet, F _p [erg s⁻¹ cm⁻² cm⁻¹], is then given by

$$\begin{eqnarray}&&{F}_{p}=\pi {I}_{p},\end{eqnarray} \tag{ 30 }$$

as per the definition of converting specific intensity to flux from a spherical surface.

The planet-to-star ratio F _p /F_⋆ is then given by using a stellar model atmosphere (e.g., Kurucz & Bell 1995), interpolated to the same wavelength grid of the simulation. To produce thermal phase curves, the simulation is run on the GPU several times for different viewing angles and the results reconstructed after the simulation is complete.

By directly accounting for the volume of the emitting regions, gCMCRT avoids the so-called wavelength-dependent "photospheric radius problem," present when modeling emission spectra from nonuniform T–p structures across a globe (e.g., Dobbs-Dixon & Cowan 2017; Fortney et al. 2019). gCMCRT therefore correctly weights the hotter and cooler parts of the atmosphere contribution to the end spectrum, avoiding known biases with blended hotter and cooler T–p profiles (e.g., Feng et al. 2016; Taylor et al. 2021).

We include the ability in gCMCRT to compute 3D contribution functions. The contribution function is generally given by (e.g., Drummond et al. 2018)

$$\begin{eqnarray}&&\mathrm{CF}={B}_{\lambda }(T)\displaystyle \frac{d\exp (-\tau /\mu )}{d{\mathrm{log}}_{10}P},\end{eqnarray} \tag{ 31 }$$

where the ${\mathrm{log}}_{10}P$ term scales the contribution per decade of pressure level. For the height vertical grid, spherical 3D geometry case, and with the random starting positions of MCRT packets, this pressure term is nontrivial to calculate. Instead, in gCMCRT, for albedo spectra, emission spectra, and phase curves, the contribution function calculation performed by tracking the total contribution each cell made to the output spectra and then normalizing by the total escaped energy fraction,

$$\begin{eqnarray}&&\mathrm{CF}(\theta ,\phi ,R)=\displaystyle \frac{f(\theta ,\phi ,R)}{\sum f(\theta ,\phi ,R)}.\end{eqnarray} \tag{ 32 }$$

A similar counter is used for the scattering component, therefore allowing the discretion between the thermal and scattering contributions to the result. This contribution function is different from Equation (31), as it directly gives the fractional contribution of that cell to the end spectra and is not normalized to the column maximum contribution. The same method can also be used to calculate the contribution functions for albedo spectra.

Accompanying gCMCRT, we provide an opacity mixer and interpolator (optools) based on the Lee et al. (2019) model with some updates. This calculator takes in the same flattened 1D array of temperature, pressures, mean molecular weight, and gas species volume mixing ratios as gCMCRT and calculates the gas-phase line opacities from precalculated cross section or k-coefficient tables. Currently this code is parallelized using openMP and requires CPUs to operate; however, plans to overhaul optools for GPU computing are underway.

We create custom k-tables using cross sections calculated using the HELIOS-K (Grimm & Heng 2015; Grimm et al. 2021) opacity calculator, as well as the EXOPLINES opacities (Gharib-Nezhad et al. 2021) spanning YY wavelength H bands over 0.3–30 μm at R = 100 resolution, suitable for producing JWST predictions. These k-tables use an 8+8 g-ordinance scheme, with 8 g-ordinances between g = 0–0.9 and 8 between g = 0.9–1.0, the same as petitRADTRANS (Mollière et al. 2019). The line lists used to generate the cross sections for the k-tables are given in Table 1. The original cross section data were calculated at a resolution of 0.01 cm⁻¹. NEMESIS (Irwin et al. 2008) formatted k-tables available from ExoMol (Chubb et al. 2021) can be also used directly within the opacity mixer, which were benchmarked extensively in Lee et al. (2019). We also include conversion scripts for the ARCiS (Min et al. 2020) and petitRADTRANS (Mollière et al. 2019, 2020) k-tables and the TauREx (Al-Refaie et al. 2021) cross section tables, also available on ExoMol, in the gCMCRT format. Individual species k-tables are mixed using the random overlap method (e.g., Lacis & Oinas 1991; Amundsen et al. 2017). This capability allows the user to take the same opacities commonly used in 1D modeling codes and apply them directly to the 3D model.

We also provide the ability to interpolate premixed k-tables to the atmospheric T–p structure. The GitHub repository contains a premixed 1× solar metallicity opacity k-table, one with TiO and VO and one without, suitable for general post-processing use. Equilibrium condensation is included in the table. Other premixed tables are available from the repository as well.

We take CIA opacities from the HITRAN database (Karman et al. 2019) and options to calculate H⁻ bound-free and free–free opacities following John (1988) and He⁻ free–free opacity following Kurucz (1970). For Rayleigh scattering opacities we include cross section and refractive index data from various sources, namely, Kurucz (1970), Sneep & Ubachs (2005), Irwin (2009), Thalman et al. (2014), and Wagner & Kretzschmar (2008) for H₂O. Thomson cross sections from free electrons can also be calculated as a Rayleigh scattering contribution.

For cloud opacities, single scattering albedos, and asymmetry factors, we use the Mie calculator MieX (Wolf & Voshchinnikov 2004). Various analytical cloud size particle distributions from the literature can also be calculated with optools (e.g., lognormal, Ackerman & Marley 2001; potential exponential, Helling et al. 2008), in addition to numerical size distributions calculated by codes such as DRIFT (e.g., Helling et al. 2008). Bin-based microphysical cloud models such as CARMA (Powell et al. 2018; Gao et al. 2020) can also be processed using optools. CIA, Rayleigh, and cloud opacities are calculated at the bin center wavelengths in correlated-k mode.

For high-resolution spectral modeling we include the option to include the Doppler shift of lines in the line of sight due to wind velocity and planetary rotation. Taking the horizontal, meridional, and vertical velocity (u, v, w [cm s⁻¹], respectively) input from the GCM model, the line-of-sight velocity, v_LOS [cm s⁻¹], is given by (Wardenier et al. 2021, Appendix B)

$$\begin{eqnarray}\begin{array}{rcl}{v}_{\mathrm{LOS}} & = & u\sin ({\theta }_{v})\sin (\phi -{\phi }_{v})\\ & & +\left[v\cos (\theta )-w\sin (\theta )\right]\sin ({\theta }_{v})\cos (\phi -{\phi }_{v})\\ & & -\left[v\sin (\theta )+w\cos (\theta )\right]\cos ({\theta }_{v})\\ & & +\,{\rm{\Omega }}\left({R}_{0}+z\right)\sin (\theta )\sin ({\theta }_{v})\sin (\phi -{\phi }_{v}),\end{array}\end{eqnarray} \tag{ 33 }$$

where Ω [rad s⁻¹] is the rotation rate of the planet, R₀ [cm] the lower boundary planetary radius, z [cm] the cell altitude, θ the latitude of the cell, and ϕ the longitude of the cell. θ_v and ϕ_v are the viewing latitude and longitude, respectively. The opacity in each cell is then interpolated to the Doppler-shifted effective wavelength, λ_eff [μm],

$$\begin{eqnarray}&&{\lambda }_{\mathrm{eff}}={\lambda }_{0}\left(1-\displaystyle \frac{{v}_{\mathrm{LOS}}}{c}\right),\end{eqnarray} \tag{ 34 }$$

where λ₀ [μm] is the rest-frame wavelength and c [cm s⁻¹] the speed of light.

The Doppler shifting of the opacities in each cell is performed at runtime inside gCMCRT. We employ a moving wavelength window method, where, rather than read in all opacity data across the required wavelength range to interpolate to, opacity data across a positive and negative wavelength range from the rest wavelength are retained in memory and Doppler-shifted opacity interpolated to this window. The window is then moved one wavelength forward for the next wavelength iteration, retaining the rest wavelength in the central part of the block. This avoids the large memory requirements when storing large amounts of opacity data across the 3D GCM grid. The size of this moving window can be changed to optimize the calculation, depending on how shifted the opacities are in the line of sight.

gCMCRT contains no native way to produce gas-phase abundances from GCM model output. However, we include 2D T–p grids of gas-phase volume mixing ratios and mean molecular weights in chemical equilibrium (CE) calculated using the ggChem (Woitke et al. 2018) code. These chemical tables can then be interpolated to the GCM T–p profiles using a Python script. This grid is what is used in the current study to find gas-phase abundances when the GCM does not contain chemical information. Alternatively, chemical information can be converted directly from the GCM or other chemical scheme to the gCMCRT profile format.

We first check the algorithm for the albedo spectra mode by comparing the numerical results of gCMCRT to known analytical solutions. In Figure 1 we compare the geometric albedo to the analytical solutions in Heng et al. (2021), who presented scattering functions for certain phase functions with isotropic multiple scattering. We also compare to the Dlugach & Yanovitskij (1974) HG solutions at high g (g = 0.9). Our model is able to reproduce the analytical solutions well, with only larger deviations occurring at very high single scattering albedo (ω > 0.99).

Zoom In Zoom Out Reset image size

For a 3D in situ demonstration we use GCM output from the Roman et al. (2021) grid of cloudy models. We post-process two models: the T_eq = 2250 K, g = 10 m s⁻¹ nucleation-limited model with extended clouds, and the same model but with compact clouds (Roman et al. 2021). We assume the effective particle size given by the vertical size-dependent scheme in Roman et al. (2021). For the cloud opacity calculation we assume a lognormal distribution with variance 0.1 μm as in Roman et al. (2021).

Figure 2 shows the results of post-processing the two GCM simulations. It is clear that the extended clouds and compact clouds produce dramatically different results. The extended clouds show a Mie scattering slope from the optically thick high-altitude clouds, of particle size ∼0.1 μm. This particle size has a maximum scattering efficiency at approximately 0.6 μm, the size parameter of 0.1 μm particles. The compact clouds produce more wavelength-dependent albedo spectra, due to the efficiency of scattering from larger particle sizes deeper in the atmosphere. The spectra show gas extinction features where the gas is more opaque, given by the packets scattered off the compact, deep cloud layer and escaping back through the gas toward the observer.

Zoom In Zoom Out Reset image size

In Figure 2 we also produce the effect of different scattering phase functions for the cloud particles. We perform the MCRT with five phase function tests: one cloud-free (without cloud opacity), one with the TTHG function, the HG function, Rayleigh function, the Draine (2003) function, and an isotropic function for the scattering from cloud particles. From the results, the calculated geometric albedo can vary significantly depending on the chosen cloud particle scattering phase function. The Rayleigh function provides the largest albedo owing to the equal forward-scattering and backscattering fractions. The TTHG and HG function are generally similar, except at optical wavelengths, where the backscattering term in the TTHG becomes more important. The Draine (2003) function performs as expected, with more HG-like behavior at optical wavelengths (and higher g), which returns to the Rayleigh scattering function as g decreases at IR wavelengths.

For emission spectra demonstration we post-process output from an HD 189733b–like THOR model (Deitrick et al. 2020) and an HD 209458b–like simulation from the grid of Baeyens et al. (2021), using the Carone et al. (2020) model (ExoRad). In Figure 3 we show the results of the dayside emission spectra of the two models. Different absorption features can be distinguished in each model. Furthermore, the HD 189733b spectrum is generally lower, as expected given the lower equilibrium temperature of this planet. In Figure 4 we show the contribution functions of the two models at the substellar point. These plots show the fraction that the atmospheric column at the substellar point contributed directly to the end spectrum, which is the contribution from all columns. These show similar contribution functions for each model, showing that the flux emanates from similar pressures at the substellar point.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In this section we perform post-processing of the HD 209458b GCM models presented in Drummond et al. (2020). Drummond et al. (2020) investigated the influence of nonequilibrium chemistry by performing two simulations, one assuming CE and one with a kinetic chemical network.

We post-process both simulations using directly the chemical abundances from the GCM output. In Figure 5 we show the transmission spectrum of the CE and kinetic chemistry models performed in Drummond et al. (2020). We are able to directly reproduce the findings in Drummond et al. (2020), recreating the differences in the CE and kinetic models near 4.5 μm due to the reduction of CO₂ in the kinetics scheme (Drummond et al. 2020) and near 10 μm due to an increase in the NH₃ abundance.

Zoom In Zoom Out Reset image size

As a demonstration of the high-resolution phase curves, we post-process a WASP-77Ab GCM simulation that was performed using the SPARC/MITgcm model (Showman et al. 2009; Parmentier et al. 2021). We chose the IGRINS instrument (Levine et al. 2018) H (1.45–1.8 μm) and K (1.96–2.45 μm) bands, which have a native resolution of R = 45,000. We therefore post-process the simulation at 3× the native resolution at R = 135,000, which was found to be high enough for modeling the cross-correlation signal, with 25,517 wavelength points in the H band and 30,675 wavelength points in the K band. We model the planet following the observational strategy of Line et al. (2021), who observed WASP-77Ab with IRGINS in the pre-eclipse phases between 0.325 and 0.47.

Figure 6 shows the post-processing at high resolution of the full spectrum for the pre-eclipse phases and a zoomed-in portion between 2.3225 and 2.3240 μm. From this figure the gross behavior of the simulated atmosphere can be seen, where both the shifting of the flux of planet with phase increases and decreases in the absolute flux.

Zoom In Zoom Out Reset image size

We perform a synthetic cross-correlation on our modeled emission spectra phase curves, interpolating the modeled phases to a phase grid corresponding to 500 s exposures during the orbit time series with 0.31 < φ < 0.47 of the planet. We produced H₂O and CO spectral templates by computing emission spectra at 0.5 phase, each with only H₂O+ continuum and CO+ continuum opacity sources. We then produced a high-resolution stellar spectrum for WASP-77Ab using SME (Piskunov & Valenti 2017, version 580, private communication), which was then broadened to a projected rotational velocity of $v\sin i=4.0\pm 0.2$ km s⁻¹ (Bonomo et al. 2017), corrected for the Keplerian (K = 0.3234 km s⁻¹; Cortés-Zuleta et al. 2020), as well as the systemic velocity (1.6845 ± 0.0004; Maxted et al. 2013), resulting in spectra in the barycentre rest frame. The modeled emission spectra were then shifted to the same rest frame for each exposure by correcting for the planetary motion and the systemic velocity and added to the high-resolution stellar spectrum. Both spectra were then interpolated onto the wavelength grid of IGRINS, and noise with signal-to-noise ratio (S/N) = 200, assuming minimal S/N as found by Line et al. (2021), and modeled tellurics were added. The spectra were then split into the 54 IGRINS orders and the cross-correlation performed using the techniques following Hoeijmakers et al. (2019), Hoeijmakers et al. (2020), and Prinoth et al. (2022). Figure 7 shows the results of the emission signals, showing a detection of both H₂O and CO in our modeled planetary emission spectra. We report a slight redshift of a few kilometers per second in both molecules from the GCM model results. This is to be expected physically, since in transmission spectra the spectra are blushifted owing to the equatorial jets being weaker on the western limb compared to the eastern limb. When viewed from the dayside hemisphere, a redshift is therefore expected. Line et al. (2021) report a blueshift in their modeling results; however, this may be a result of using 1D models in their analysis (M. Line, personal communication).

Zoom In Zoom Out Reset image size

For this demonstration we use output from the HD 209458b Exo-FMS GCM results (Lee et al. 2021). We take the correlated-k model from Lee et al. (2021). We post-process this simulation to produce high-resolution transmission spectra across the GIANO-B wavelength range (≈0.95–2.5 μm; e.g., Giacobbe et al. 2021). We use a resolution of R = 150,000 (for a total of 142,108 wavelength points) and produce spectra with and without Doppler and orbital shifting assuming that the planet is at an orbital phase of zero. The Doppler-shifted spectrum across the wavelength range is shown in the left panel of Figure 8.

Zoom In Zoom Out Reset image size

We then perform a cross-correlation analysis using the Doppler-shifted spectrum as the mock observable and the non-Doppler-shifted spectrum as the template. The result of this is shown in the right panel of Figure 8. A clear net blueshift of ≈3 km s⁻¹ can be seen in the peak of the cross-correlation function, corresponding to the wind speeds found in the simulation (Lee et al. 2021).

This simple demonstration shows that gCMCRT is suitable for high-resolution transmission spectra modeling studies similar to Miller-Ricci Kempton & Rauscher (2012), Flowers et al. (2019), and Savel et al. (2021). Wardenier et al. (2021) produced a detailed analysis of WASP-76b Fe depletion (Ehrenreich et al. 2020; Kesseli & Snellen 2021) as a function of phase by post-processing a SPARC/MITgcm WASP-76b model using the CPU version of gCMCRT.

With this new version of gCMCRT we have greatly extended the capabilities for RT modeling 3D planets at high spectral resolution. Most 3D high-resolution modeling studies to date (e.g., Miller-Ricci Kempton & Rauscher 2012; Showman et al. 2013; Flowers et al. 2019; Beltz et al. 2021; Harada et al. 2021; Wardenier et al. 2021) have focused on narrow-wavelength H bands, with the exception of Savel et al. (2021). We have shown that gCMCRT can efficiently model across the IGRINS wavelength range (≈1.45–2.45 μm) and GIANO-B (≈0.95–2.5 μm) at moderately high resolutions of R = 135,000 and R = 150,000, respectively, with the number of discrete wavelength points numbering in the hundred thousands, suitable for cross-correlation of model results to elucidate physical processes. This enables a much more detailed physical interpretation of the results of contemporary high-resolution studies, spreading across many more lines and probing the features of more species.

In Section 12 we performed mock cross-correlation detections of H2O and CO. This would be relatively simple to extend to other species (such as HCN) by making template spectra and performing the same cross-correlation analysis. In Section 13 we only cross-correlated the gross model spectrum to find the net wind speed of the model. By generating template spectra of each constituent species, it would be relatively simple to extend our approach here to model the detection of specific molecules and atoms directly from the GCM output, for example, producing predictions of the detected molecules in Giacobbe et al. (2021). Our high-resolution modeling capability provides a useful additional check to GCM models, where theorists and modelers can now directly compare their GCM output to the high-resolution observational data and see whether their predictions about the 3D distribution of chemical species abundance, wind speeds, condensation, and haze patterns are reasonable.

We specifically used outputs from many different modeling groups in this study. This is to demonstrate the applicability of the model to many GCM output formats and to provide the community with baseline conversion scripts they can use to automate the processing of large grids of GCM data.

Our setup and code of gCMCRT enable future implementations of other photon microphysics. For example, heating and cooling rates of the atmosphere can be performed using the Lucy (1999) method. Well-used CMCRT techniques from the astrophysical community such as photoionization and dissociation (Wood et al. 2013) can also be adapted for the planetary regime using this model as a baseline, enabling exploration of complex photon microphysics in 3D. Expanding the model to performing 1D and using different geometries such as Cartesian will further increase the flexibility of the model.

gCMCRT is a GPU-accelerated MCRT code, suitable for global processing of 3D atmospheric data. We develop standard albedo, emission, and transmission spectra modes for typical use in producing observable properties from GCM model output. We also develop a high-resolution spectral mode for emission and transmission spectra, able to capture the Doppler shifting of spectral lines due to rotation and atmospheric winds.

As a fully 3D model, gCMCRT avoids the biases and assumptions present when using 1D models to process 3D structures. Atmospheric layers are weighted appropriately to their contribution to the end spectra.

gCMCRT is well suited to performing the post-processing of large parameter GCM model grids starting to be performed by members of the community. We developed simple pipelines that convert the 3D GCM structures from many well-used GCMs in the community to the gCMCRT format, interpolate chemical abundances (if needed), and perform the required spectra calculation. We include the ability to convert opacity tables from ExoMol to the gCMCRT format, enabling greater intercomparison with well-used 1D models.

The high-resolution spectra modes of gCMCRT provide an additional highly useful capability for 3D modelers to directly compare output to high-resolution spectral data. Simulated wind speed, molecular detections, and other variables can be calculated and compared to the observations through mock cross-correlation analysis of the models.

15.1. Support

We thank M. Line for details and help with the IGRINS WASP-77Ab observations. E.K.H.L. was supported by the SNSF Ambizione Fellowship grant (No. 193448). J.P.W. sincerely acknowledges support from the Wolfson Harrison UK Research Council Physics Scholarship and the Science and Technology Facilities Council (STFC). M.T.R. was supported by a European Research Council Consolidator grant, under the European Unions Horizons 2020 research and innovation program (No. 723890). 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 HPC support staff at AOPP, University of Oxford, and University of Bern are highly acknowledged.