Modeling Atmospheric Lines By the Exoplanet Community (MALBEC) version 1.0: A CUISINES radiative transfer intercomparison project

MALBEC is part of the CUISINES model intercomparison work group. CUISINES stands for Climates Using Interactive Suites of Intercomparisons Nested for Exoplanet Studies and provides a framework for different model intercomparison projects, as well as opportunities to find synergies (pairings) between different work groups. It currently consists of nine different projects, including chemical, circulation, energy balance, and planet specific models. For more information, see https://nexss.info/cuisines. The framework started with the TRAPPIST-1 Habitable Atmospheres Intercomparison (THAI) project, which presented a trilogy of GCM studies investigating dry and wet climate for TRAPPIST-1e (Turbet et al., 2022; Sergeev et al., 2022; Fauchez et al., 2022). CUISINES was inspired by the THAI success, and led to two international workshops and eight more exoMIPs at the time of writing. Three other exoMIPs — SAMOSA, CAMEMBERT, and FILLET — have currently published protocol papers describing their efforts:

• •

  The Sparse Atmospheric Modeling Sample Analysis, or SAMOSA, aims to provide guidance for what type of planets need to be modeled using full 3D general circulation models and where more simple 1D equivalent models are sufficient (Haqq-Misra et al., 2022). The parameters investigated are insolation and surface pressure and it compares the outcome scenarios between the ExoCAM 3D GCM and ExoPlaSim 1D models (Wolf et al., 2022; Paradise et al., 2022).

• •

  Comparing Atmospheric Models of Extrasolar Mini-Neptunes Building and Envisioning Retrievals and Transits, CAMEMBERT, focuses on inter-comparing the results from eight different GCMs simulating the atmospheres of GJ1214b and K2-18b (Christie et al., 2022).

• •

  Functionality of Ice Line Latitudinal EBM Tenacity, FILLET, studies energy balance models, and in particular looks at climate ‘end’ states using different values for isolation, obliquity, starting temperatures and atmospheric composition (Deitrick et al., 2023).

MALBEC serves several roles within the CUISINES framework: as a down-stream framework to produce spectra and simulate observables from the 3D GCM exoMIPs (THAI, CAMEMBERT, MOCHA), constrain the detectability of outputs from photochemical codes (PIE), and verify the radiative fluxes from energy balance models (FILLET). Central in all of these is the MALBEC configuration file (see below), which provides a standardized and accessible format for both input and output. Namely, “up-stream” models that generate profiles of the atmospheric variables (e.g., GCMs, RCE models) can produce output files in the MALBEC format that can then be read by the “down-stream” radiative transfer codes. This standardized file provides a bridge for testing the outputs of different models, and is how we ensure the test cases in this protocol paper are interpreted similarly by all codes. Additionally, MALBEC can provide two metrics to other exoMIPs: the detectability of recorded differences between up-stream codes, and the spread or uncertainty in the down-stream codes.

The objective of this paper is to present the MALBEC experimental protocol that will be used to compare 10 already participating models. A clear and descriptive protocol, separated from the result papers, as well as broad participation is essential for a meaningful model intercomparison project. This paper is structured as follows. In Section 2, we outline the motivation for the MALBEC exoMIP, while Section 3 describes the physical and computational considerations that came across during the development of the study, as well as the input/output files. In Section 4 we describe the models currently onboard MALBEC, and the experimental protocol is addressed in Section 5. The MALBEC repository, the output formats, standards and archiving are discussed in Section 6. The summary is given in Section 7.

Throughout the wavelength range under investigation in MALBEC, many processes affect the radiation propagating through an atmosphere. Table 1 highlights these and in which spectral regions they are most relevant. Spectroscopic linelists are at the core of radiative transfer calculations, and in many cases, they are the dominant source of discrepancy between models. Considerations regarding the inclusion of isotopes, of high-energy lines for warm atmospheres, of lineshape models and wing-cutoffs can have substantial effects on the computed spectra. In addition, how these spectroscopic linelists are integrated into the RT model, such as line-by-line, cross-sections and correlated-k methods, can lead to different levels of fidelity.

Based on the spectroscopic parameters (including strength and energies) of the different ro-vibrational molecular transitions, the absorption of a particular gas can be determined for a temperature-pressure combination. The absorption lineshapes also depend on the properties of the dominant collisional partner in the atmosphere (e.g., air, CO${}_{2}$, H${}_{2}$, He). Different databases are available, with HITRAN (Gordon et al., 2022) being one of the most comprehensive databases, containing 55 molecules and 148 isotopologues at different collisional regimes, and about 300 in experimental cross-sections. The focus of HITRAN is Earth’s atmospheric conditions, and for that reason molecular states and transitions that are relevant at different atmospheric temperatures and pressures are not always included. For higher temperatures, the HITRAN team created the HITEMP database (Rothman et al., 2010), which contains spectroscopic information of amongst other \ceH2O, \ceCO2, CO (Li et al., 2015), \ceN2O (Hargreaves et al., 2019) and \ceCH4 (Hargreaves et al., 2020), molecules particularly relevant for exoplanet atmospheres. Additionally, the ExoMol repository (Tennyson et al., 2016; Chubb et al., 2021) is specifically focused on exoplanets, brown dwarfs, and cool stars, and hosts a large spectroscopic databases of molecules, which are particularly suited for high temperature and non-LTE regimes. Specialized databases that combine several repositories also exist, such as HELIOS-K (Grimm et al., 2021) [https://dace.unige.ch/opacity], PSG databases Villanueva et al. (2022) [https://psg.gsfc.nasa.gov/helpatm.php], or MAESTRO (Batalha & Rooney, 2020) [https://science.data.nasa.gov/opacities/app].

RT codes have different ways of handling linelists, and line shape models - including different assumptions regarding line wing cutoffs - can have substantial impact on the computed spectra. The most direct method is to calculate the absorption by line-by-line shape modeling, which is very accurate but computationally expensive, and it is not always possible in the UV/optical since linelists typically do not accurately cover this region. Some RT codes rely on UV/optical cross section tables to complement the line-by-line/correlated-$k$ tables. The MPI-Mainz provides extensive cross section data for many species in this spectral range (Keller-Rudek et al., 2013). For low resolution studies (i.e., resolving powers less than 5000), the absorbance is often saved in correlated-$k$ tables, which capture the distribution of absorbance in a spectroscopic bin in a selected number of points. Correlated-$k$ tables are pre-calculated for a grid of temperature and pressure combinations and interpolated upon running RT calculations. Alternatively, cross-section data can be used as well, which are pre-calculated absorption spectra for particular temperature-pressure conditions, but the RT calculation need to be performed at very high spectral resolution to ensure high fidelity (Garland & Irwin, 2019). Isotopes are available both in HITRAN and ExoMol. Even when using the same databases/methods, the inclusion and treatment of isotopologues among the RT codes varies greatly, mainly due to the assumed linelist, and the considered isotopic abundances. For most isotopes, spectral shifts are minor and are typically lost under the main isotope’s absorbance features (e.g., ${}^{13}$CH${}_{4}$), but in some cases, absorbance of the heavier isotopes is in absorbance minima of the main isotope and can thus be important to consider (e.g., HDO).

At short wavelengths, Rayleigh Scattering has a strong wavelength dependency, resulting in a smooth increase from the near-IR, visible to UV, following $\lambda^{-4}$. The magnitude depends on the polarizability of the atoms/molecules, but in principle all (dominant) gases in an atmosphere scatter light in this way. The formulation used in RT codes is often the one described in Sneep & Ubachs (2005). The tabulated values are available from NIST¹¹1https://cccbdb.nist.gov/pollistx.asp, yet different codes implement this process in slightly different ways. In general we have encountered relatively small differences on this process across the “casks” compared here.

Primarily in the infrared, but also in the UV and in the visible, Collision Induced Absorption (CIA) arises from temporary bi-molecular complexes, resulting in relatively broad features whose intensities scale quadratically with pressure. Description of the phenomena and how to calculate the corresponding absorption cross-sections can be found in Karman et al. (2019). The HITRAN CIA repository ²²2https://hitran.org/cia/ lists the most generally used CIA lists. Depending on the main emphasis of each RT model, the number of CIA databases and considered repositories varies greatly. RT models with an emphasis on Earth-like planets implement CIA lists well beyond HITRAN, including water continuum and capturing several \ceN2-x and \ceO2-X CIA bands, while gas-giant RT models mainly consider several theoretical and lab-based CIA databases for \ceH2-X and He-X pairs.

During the development of this protocol, we investigated a number of different treatments of the computational aspects of the radiative transfer calculations, to understand the fundamental differences between codes. To ensure all models interpreted the different tests (“casks”) similarly, a common input file was agreed upon. This “configuration” file, called the MALBEC file, contains all the information required to calculate the different test cases, in a consistent format across the different “casks”, and can be ingested by each RT code within MALBEC.

Not all radiative transfer codes discretize the atmosphere in the same way, so a method was developed to homogenize the layering schemes (pressure and temperature versus height) for all RT codes. In addition to providing anchor points at specific altitudes to describe the local atmospheric pressure, temperature, and molecular/aerosols abundances, we also provide polynomial fits to these quantities, so each code can discretize the atmosphere in any preferred manner. This is relevant both for the radiative transfer calculations in each layer, and for the path lengths of light through each layer (i.e., ray tracing). Particularly for ray-tracing calculations in spherical refractive atmospheres, it is important to sample at sufficiently high resolutions, either by having enough layers or by creating sub-layers. If insufficient layers are used this may lead to large differences. In the configuration file we provide “anchor” values, defining the specific pressure, temperature and abundances at specific points in the atmosphere. These points quantify the properties of the atmosphere at that exact pressure point, yet when performing ray-tracing calculations to determine the path length of individual layers, each model may sub-divide the light path between anchor points into much finer points by using a sub-layering scheme in order to ensure an accurate sampling of the large variability of the parameters in-between “anchor” points. Additionally, we also provide analytical fits (3rd order vs. log${}_{10}$ pressure [bar]) to those profiles, so the user can create their own vertical discretization, ensuring accurate results with their respective model.

The configuration file defines the following:

• •

  Planetary and Stellar parameters: Planetary radius [m], surface gravity [m/s${}^{2}$], phase at which the planet is observed, planet temperature [K], albedo, emissivity, star-planet distance [m], stellar type, stellar temperature [K] and stellar radius [m]. Albedo and emissivity are defined to be constant across the simulation wavelength range, and the surface scattering model is assumed to follow the Lambert diffuse principle.

• •

  Opacity sources opacity sources used: linelists, CIAs, Rayleigh, etc.

• •

  Atmosphere-columns: this entry specifies the variables included in the simulation and the anchor points. Pressures are in [bar], temperatures in [K], mean molecular weight in [g/mol], height in [km], and molecular abundances in [mol/mol].

• •

  Layers: the specifics of each layer, containing the pressure, temperature, mean molecular weight, altitude, and molecular abundances following the definitions described above.

• •

  Analytical functions layering: For a number of relevant species, 3rd order polynomial parameters describing the abundance with the log${}_{10}$ of the pressure are given.

• •

  Bottom of the atmosphere (BOA): the atmosphere is defined to start at the first anchor point, with surface properties as listed in the configuration file (surface temperature, albedo, emissivity). The planetary gravity and radius refer to this point.

• •

  Top of the atmosphere (TOA): the atmosphere is defined to end at the last anchor point, and therefore the atmospheric radiative transfer is only computed between BOA and TOA.

• •

  Spectral range and sampling: four parameters define the range and resolution of each simulation. The minimum wavelength defines the center wavelength of the first “pixel”, while maximum wavelength defines that of the last “pixel”. The separation between pixels can be a fixed wavelength if the unit is “µm”, while if the unit is “RP”, the separation assumes a constant resolving power ($\frac{\lambda}{\delta\lambda}$).

• •

  Output units: for transit simulations, the output model units are in ($\frac{R_{planet}}{R_{star}})^{2}$, while for direct imaging simulations, the model outputs are in [W/sr/m${}^{2}$/µm].

• •

  Opacities: the models should employ all opacity terms suitable for each test case, for instance high-energy linelists would be beneficial when simulating hot atmospheres as in case T3B (see below). This includes molecular opacity (for all isotopes), Rayleigh scattering/extinction, water continuum (e.g., MT$\_$CKD) and Collision-Induced-Absorptions (CIA). Molecular opacities should be generated employing the most appropriate linelist available to each model for that collisional regime (e.g., \ceH2, air broadening). All simulations as part of this intercomparison are at moderate resolutions (with resolving powers equal or lower than 500), and therefore correlated-$k$ or other pre-computing methods are encouraged.

APOLLO is a 1D radiative transfer and MCMC spectroscopic retrieval code designed for modularity and internal model intercomparisons. The main description of APOLLO can be found in Howe et al. (2022), and the model can be used to synthesize both transit and emission spectra. It is built on a philosophy of free retrieval that is agnostic about specific microphysics and as such allows the user to select from several options for prescriptions for temperature-pressure profiles, cloud structures, and molecular inventories. While multiple radiative transfer prescriptions are also available, by default, APOLLO simulates radiative transfer (reflected and emission) through an atmosphere using the two-stream approximation (Toon et al., 1989a) with hemispheric closure (Meador & Weaver, 1980). APOLLO uses the Freedman et al. (2014) cross section tables on a species-by-species basis. A weighted sum of the cross sections for each species is then computed based on the mixing ratios to compute optical depths layer by layer, including CIAs and aerosols. For transit spectroscopy, the starlight is assumed to be parallel. For emission spectra, the emitted light is integrated over an eight-point Gaussian quadrature. APOLLO can perform retrievals by an MCMC method using the emcee Python package (Foreman-Mackey et al., 2013). APOLLO also includes an “ensemble” method to rapidly generate and evaluate regularly-spaced forward model grids.

BART (Harrington et al., 2022; Cubillos et al., 2022; Blecic et al., 2022) is an open-source exoplanet retrieval code which pairs the transit RT code (Rojo, 2006) with the Multi-Core Markov-Chain Monte Carlo (MC3; Cubillos et al., 2017) Bayesian framework. BART performs 1D radiative transfer under the assumption of hydrostatic balance and LTE. For transmission spectra it performs ray tracing, and for emission spectra it assumes a plane-parallel atmosphere. Its programmatic design focuses on flexibility, enabling BART to replicate the setups of other codes and achieve agreement in results (e.g., Himes & Harrington, 2022). The RT routines are implemented in C for speed, while BART is implemented in Python for ease of use and portability. Both transit and BART are generally executed via configuration files, where the user may specify the various parameters of the run, including the spectral range and resolution, atmospheric species and profiles, Rayleigh extinction and/or cloud models, and the observing geometry (transit, eclipse, or direct observations). At present, BART does not consider multiple scattering when computing emission or transmission spectra. Transit’s package to ingest linelists allows for HITRAN/HITEMP (Rothman et al., 2010; Hargreaves et al., 2020; Gordon et al., 2022), ExoMol (Tennyson et al., 2016; Chubb et al., 2021; via Repack, Cubillos, 2017), TiO (Schwenke, 1998), and the \ceH2O list of Partridge & Schwenke (1997). Transit utilizes a line-by-line opacity sampling approach to reduce the sizes of extensive linelists while still retaining an accurate RT calculation (see Rojo, 2006, Harrington et al., 2022, and Cubillos et al., 2022 for more details on other computational cost-saving measures). In Figure 2 of test cask T3B, we utilize the HITEMP linelists for H${}_{2}$O and CH${}_{4}$; HITRAN CIA databases for H${}_{2}$-He, H${}_{2}$O-H${}_{2}$O, CH${}_{4}$-CH${}_{4}$, CH${}_{4}$-He, and H${}_{2}$-CH${}_{4}$; and the H${}_{2}$-H${}_{2}$ CIA database of Borysow et al. (2001); Borysow (2002).During a retrieval, BART calculates high-resolution spectra and bins them according to the instruments used to observe the exoplanet. MC3 calculates the effective sample size of the posterior resulting from the inference as well as its induced uncertainty in credible regions determined from that posterior (see Harrington et al., 2022 for more details).

gCMCRT is an open source 3D spherical geometry hybrid ray tracing and Monte Carlo radiative-transfer code which uses GPUs to accelerate calculations. The main descriptions of gCMCRT can be found in Lee et al. (2019) and Lee et al. (2022). gCMCRT is primarily aimed at accurately post-processing 3D GCM data, especially when multiple scattering from cloud particles is important. gCMCRT directly simulates the path of photon ‘packets’ (usually $\sim$ 10${}^{6}$) through the prescribed 3D spherical geometry structure. This is a sub-grid model, with packets traversing through the 3D volume of each cell. For transmission spectroscopy, packets are fired at the dayside of the planet and traced through the transmission limbs, which are then used to estimate the transmission spectra integral through random sampling. For emission spectra, packets originate in the 3D spherical volume itself, and are traced towards a certain viewing angle (planetary phase). The fraction of energy escaping towards this direction is then tracked and used to calculate the emission spectra. For albedo calculations, the scattered fraction of the incoming starlight is calculated towards a certain viewing direction. gCMCRT can also perform multiple-scattering in 3D geometry by sampling scattering phase functions (for example, the Rayleigh or Henyey-Greenstein phase functions) when a packet undergoes a scattering event. gCMCRT contains a well tested opacity mixer (optools), which can mix and prepare gas phase abundances, correlated-$k$ tables, line-by-line, CIA, Rayleigh and cloud opacities and scattering properties. gCMCRT is able to perform high-resolution spectral modelling including Doppler shifting of lines, allowing simulated GCM data to be used as part of high-resolution cross-correlation studies.

petitRADTRANS (pRT for short) is a publicly available package to calculate transmission and emission spectra of clear and cloudy exoplanets (see Mollière, P. et al., 2019; Mollière et al., 2020, for more details). The core routines are written in FORTRAN and wrapped in a Python package. The model can run at high-resolution ($\lambda/\Delta\lambda\leq 10^{6}$) using a line-by-line treatment, as well as at low resolution ($\lambda/\Delta\lambda\leq 10^{3}$) using a correlated-$k$ approach. Opacities can be rebinned at any resolution needed by the user, either directly in petitRADTRANS (for high resolution), or through the exo-k package (Leconte, 2021) which is called internally (for low resolution). The default set of opacities provided by the maintainers includes opacities from HITRAN/HITEMP and ExoMol, generally spanning a temperature range between 80 and 3000 K, and a pressure range between $10^{-6}$ and $10^{3}$ bar. Further opacity tables can be downloaded from ExoMol directly in a pRT-friendly format, or calculated by the user through the ancillary scripts provided.

petitRADTRANS allows the user to consider different cloud treatments, such as gray clouds, wavelength-dependent power laws, or using opacity (in $cm^{2}/g$) of real condensates of various compositions and shapes. Arbitrary opacity functions and particle setups can also be employed. Users can also provide low-res (pseudo-continuum) opacities to petitRADTRANS to be included in the calculation, thus customizing the cloud treatment. petitRADTRANS includes the calculation of Rayleigh scattering and collision-induced absorption. It was recently updated to include the surface direct and thermal scattering for rocky exoplanets (see Appendix A in Alei et al., 2022, for details). Specifically for the computation of multiple scattering of reflected starlight and thermal emission, pRT uses the Feautrier method, which is a third-order method that allows the treatment of the radiative transfer equation in the diffusive regime, and for which pRT iterates to solve for the multiple scattering source by employing the ALI and Ng acceleration. Users can specify wavelength-dependent surface reflectivity and emissivity, allowing the scattering treatment of various surface compositions of rocky exoplanets. In its latest version, it also allows the user to run retrievals using the nested sampling approach. Thanks to the MALBEC collaboration, opacities of additional isotopes for many atmospheric absorbers (e.g., \ceCO2, CO, \ceH2O) were computed, as well as the collision-induced absorption tables for water (MT_CKD, Kofman & Villanueva 2021; Mlawer et al. 2012). These will be soon included in the main opacity repository of petitRADTRANS. For the purpose of this intercomparison, the version of petitRADTRANS was frozen to 2.6, though the package is in the process of being upgraded to version 3, with additional features and quality-of-life improvements.

The Planetary Spectrum Generator, Villanueva et al. (2018, 2022), is a publicly available radiative-transfer suite, generalized to any type of planetary atmosphere/surface across a broad range of wavelengths (50 nm to 100 mm, UV/Vis/near-IR/IR/far-IR/radio) and observational geometry (e.g., transit, direct-imaging, in-situ, orbiter, limb). The suite includes a 3D orbital calculator for most bodies in the Solar system and all confirmed exoplanets, and also includes a broad range of atmospheric and chemical models. The ray-tracing integration across the atmosphere is performed in PSG via a sub-layering spherical/refractive iterative algorithm, enabling high-fidelity results even when operating with a limited number of layers. Several radiative-transfer modules operate within PSG, including CEM (Cometary Emission Model), a non-LTE exospheres solver, and PSGDORT, the multiple-scattering radiative-transfer module of PSG that employs the discrete ordinate method to model atmosphere/surface scattering processes. PSGDORT is based on DISORT v2.1 (Stamnes et al., 2000) which were adapted for non-LTE and optimized to operate with a variety of spectral grids (e.g., line-by-line, correlated-k, surface scattering grids) as employed by the PSG radiative transfer algorithm. PSGDORT also includes correction for pseudo-spherical geometry as described by (Dahlback & Stamnes, 1991). The PSG radiative transfer modules take into account Rayleigh and Raman scattering, refraction, non-LTE processes, molecular/atomic absorptions, and collision induced absorption processes. Ro-vibrational absorption are calculated line-by-line for resolving powers higher than 5000 and correlated-$k$ tables are used for lower resolutions.

PSG operates with several databases, including the latest HITRAN/HITEMP distribution (linelist, CIAs, cross-sections), GEISA, and ExoMol. For the PSG simulation results presented in Figure 2 of cask T3B, we utilized correlated-k tables based on the HITEMP linelists of H${}_{2}$O and CH${}_{4}$. PSG complements the UV/optical range, by integrating UV/optical cross sections for dozens of species collected from a range of spectral databases. Most of the UV cross-sections originate from the MPI-Mainz Spectral Atlas (Keller-Rudek et al., 2013), which have been parsed, combined and formatted to provide a comprehensive coverage of the 0.01 to 1 µm wavelength range. Additional cross-sections include those of \ceO3 by Serdyuchenko et al. (2014), \ceCO2 by Venot et al. (2018), the Herzberg \ceO2 continuum bands as well as the \ceO2–\ceO2 absorption bands (Wulf bands) by Fally et al. (2000), the Herzberg \ceO2 band system (Mérienne et al., 2001), and several \ceO2-X CIAs (Fauchez et al., 2020). The MT_CKD water continuum is parameterized in PSG by computing CIAs tables for \ceH2O-\ceH2O and \ceH2O-\ceN2 from that water model (Kofman & Villanueva, 2021; Mlawer et al., 2012). PSG also includes a retrieval framework, optimal estimation and Multinest methods, and a noise simulator. The tool can be operated in several ways: via a public web interface³³3https://psg.gsfc.nasa.gov, via a flexible application program interface (API) accessible from any scripting language (e.g., python, IDL), and it can be also operated by an installable suite via the Docker virtualization system. Additionally all opacity tables, spectroscopic databases, and core fundamental source codes are publicly availble at the site and the PSG github respository.

Originally developed to easily handle, convert, and interpolate radiative opacities from various sources, the Exo_k python library (Leconte, 2021) now includes a 1D radiative transfer module able to model the transmission and reflection/emission spectra of a wide range of atmospheres. Thanks to its original philosophy, Exo_k can include any molecule or opacity source for which opacities can be provided in the form of either cross sections or correlated-$k$ coefficients tables in one of the many formats handled by the library. In particular, data can be provided at any resolution as the library implements a novel algorithm able to bin-down even $k$-coefficient tables (Leconte, 2021). A high numerical efficiency is ensured by the use of the the “just in time” compilation capability of the numba library.

The reflection/emission algorithm is an improved version of the two-stream framework from Toon et al. (1989b) that accounts for multiple scattering and scattering asymmetry (see the appendix of Chaverot et al. (2022) for details). Mie scattering can be modeled for any type of aerosols by providing the aerosol optical properties. The change in atmospheric scale height due to both gravity and mean molecular weight variations are taken into account, both in emission and transmission.

Pytmosph3r is the 3D implementation of Exo_k that can be used to model observables from both parameterized 3D atmospheric structures and outputs from a 3D global climate model (Caldas et al., 2019; Falco et al., 2022). Using Exo_k for all the opacity calculations, it inherits all the aforementioned features. In transmission, rays of light are shot through the limbs of the planet, possibly crossing many atmospheric columns of the atmospheric model and thus accounting for its heterogeneity, both along and across the limb. In emission/reflection, the two-stream algorithm from the 1D model is used in each atmospheric column to compute the radiative intensity in the direction of the observer – thus accounting for limb darkening. Intensities are then weighted by solid angle and integrated over the visible sphere to yield the observed flux at the desired observer’s location.

The Spectral Mapping and Radiative Transfer (SMART) code is a highly flexible one-dimensional, plane-parallel, line-by-line, multi-stream, and multiple-scattering radiative transfer model originally presented by Meadows & Crisp (1996) and Crisp (1997). SMART solves the radiative transfer equation using the discrete ordinates method through the DISORT FORTRAN code (Stamnes et al., 1988). SMART ingests gaseous absorption coefficients generated by its Line-By-Line Absorption Coefficients (LBLABC) companion model, as well as optional supplemental cross sections typically used to compensate for gaps in line-list opacity data, and separately sourced collisionally-induced absorption (CIA) coefficients. LBLABC uses the HITRAN linelists to calculate absorption coefficients (Gordon et al., 2022). Usually cross-section data derived from the MPI-Mainz database (Keller-Rudek et al., 2013) or other sources is used for UV-Vis wavelengths where linelist data is often absent; however, cross-section data can also be used at any wavelength where appropriate. The same code is used for both shortwave and longwave radiation. The version of SMART employed here uses the ray tracing model of Robinson (2017) to simulate transit transmission observations, which includes the effects of refraction. SMART has been repeatedly verified against terrestrial solar system objects (e.g. Tinetti et al., 2005; Robinson et al., 2011; Arney et al., 2014) and has a long history of being used to generate exoplanet spectra (e.g. Charnay et al., 2015; Lincowski et al., 2018; Lustig-Yaeger et al., 2019). SMART is highly customizable with user selected input gases and high internal resolution. The model can also be configured for clouds and aerosols.

“Suite Of Community RAdiative Transfer codes based on Edwards & Slingo (1996)” (SOCRATES) is an open-source radiative transfer code with a two-stream multiple scattering solver for reflected starlight and thermal emission. It uses the correlated-$k$ method to solve for gaseous absorption, with further details covered in Edwards & Slingo (1996) and Manners et al. (2022). Opacity data are sourced from a range of databases such as HITRAN, HITEMP and ExoMol (see Amundsen et al., 2014, 2016; Goyal et al., 2018, for full details). For cask T3B, correlated tables were constructed from ExoMol linelists. At the moment, the SOCRATES code was written to calculate transmission spectra for the global grid from the Unified Model (UM) — a 3D GCM developed and used by the Met Office, and at present SOCRATES cannot calculate transmission spectra outside a GCM. However, a standalone version of this calculation is planned for a near-future release of SOCRATES. Within each GCM column the transmitted solar flux is calculated using a spherical shell geometry, replacing the standard plane-parallel assumption (detailed description can be found in Lines et al., 2018). The total optical depth through the atmosphere is found by summing the mass along the slant path through each shell multiplied by the corresponding mass extinction coefficient for the layer. Each model column is treated independently and each layer is assumed to be homogeneous with identical optical properties in that column. The fluxes are calculated in the column facing away from the star in the limb. The resulting fluxes in each column are added up to represent the full transmission spectrum. With this approach, the variations of fluxes in the limb perpendicular to the observer can be fully captured but those along the line-of-sight can only be approximated. The assumption of homogeneous spherical shell might also lead to a bias if there are variations across the terminator.

TauREx (Tau Retrieval for Exoplanets) is an open-source program initially developed by Waldmann et al. (2015a, b) for both transmission and thermal emission spectroscopy. The current version in use – now completely Python-stacked – corresponds to TauREx 3 and is substantially described in Al-Refaie et al. (2021). TauREx is a fully Bayesian, line-by-line radiative transfer and retrieval framework. It can be used with any molecular line-lists in hdf5 and in pickle format, including the latest data from the Exomol project (Tennyson et al., 2016; Chubb et al., 2021), HITEMP (Rothman et al., 2014) and HITRAN (Gordon et al., 2017). Its modular, object-oriented architecture, combined with an extensive TauREx 3 Python library, makes it particularly flexible and customisable to the user. TauREx models a 1D-atmosphere and computes the path integral from a single cross-section where each (interpolated and weighted) opacity has been fused. The opacity integral is computed based on a finite difference scheme. In addition, TauREx provides a full implementation of Rayleigh and Mie scattering via extinction, as well as simple (Grey) cloud models. It has been particularly tailored for the modeling and interpretation of transit spectra, so it only considers extinction processes and it does not solve for scattering emission when computing thermal emission. The user builds both transmission and emission models by defining a temperature profile, pressure parameters, a chemistry model, and contributions to the optical depth (molecular absorption, CIA, Rayleigh, Mie scattering, simple clouds). The wavenumber grid of the forward model is selected at runtime from whichever loaded opacity with the highest resolution. Every other opacity is then resampled to the chosen opacity’s grid before any computation begins; this grid becomes the native grid of the forward model.

In this study, the T-P profile, as well as chemical abundances at each layer, were provided as inputs (thus not requiring the use of any TauREx chemical model). For the TauREx simulation results presented in Figure 2 of case T3B, we used ExoMol absorption cross-sections for H${}_{2}$O (Polyansky et al., 2018) and CH${}_{4}$ (Yurchenko et al., 2017), CIA opacities for H${}_{2}$-H${}_{2}$ (Fletcher et al., 2018), H${}_{2}$-He (Abel et al., 2011), CH${}_{4}$-X (from HITRAN), as well as Rayleigh scattering contributions (Cox, 2015).

PICASO, the Planetary Intensity Code for Atmospheric Scattering Observations, is a radiative transfer code originally designed for reflected light planetary spectra (Batalha et al., 2019), and has since been expanded to include calculations of transmission and emission planetary spectra (Batalha & Rooney, 2020; Batalha et al., 2021), 3D spectral modeling (Adams et al., 2022), phase curves (Robbins-Blanch et al., 2022), 1D climate modeling (Mukherjee et al., 2023), as well as fitting spectroscopic data to models (The JWST Transiting Exoplanet Community Early Release Science Team et al., 2022). The code is written in Python and publicly available on GitHub⁴⁴4https://github.com/natashabatalha/picaso, and has been used in a wide range of studies for exoplanets and brown dwarfs (e.g., Mukherjee et al. 2021; Foote et al. 2022; Lew et al. 2022; Kozakis et al. 2022). PICASO is based off codes originally described in McKay et al. (1989), Marley et al. (1999), and Cahoy et al. (2010), being built off of methods from Toon et al. (1977, 1989a) to compute reflected starlight and thermal emission as observed at full phase, with updates from Cahoy et al. (2010) for observations at any phase angle and including multiple-scattering via the two-streams approximation.

The code is accompanied by a database of opacities resampled down to R = 20,000 and 60,000 from line-by-line calculations on the order of R$\sim$10${}^{6}$ (see methodology in Gharib-Nezhad et al. 2021 for details). This database spans wavelengths from 0.3 to 14 µm and can either be downloaded in entirety or queried using sqlite3. The resampled opacity database draws from multiple databases such as HITRAN and ExoMol, along with many other sources for linelists, CIAs, and cross sections, with the code providing helpful tools to extract the relevant references for individual model runs. If desired, PICASO additionally provides functions for creation of correlated-$k$ tables for individual gases. Within the code there are options to add clouds in the form of grey cloud coverage or cloud model output (Rooney et al., 2022), as well as user supplied parameteriz]d arbitrary Rayleigh scattering.

The core experiments (T0A, T3B) are defined specifically to test basic assumptions within each model and to further assist with the homogenization of the results. The purpose of the MALBEC exoMIP is not for all models to agree, but to understand how certain aspects in the parameterization of the models impact the simulations. As such, issues with layering, opacities, refraction, line-lists, scattering considerations, impact the results. As long as each model is parameterized as close as possible to the MALBEC configuration file, then each model result will reveal how parameterizations and modeling issues impact the spectra.

Some of the most common differences in RT models is how each code discretizes and defines the atmospheric layering. For instance, some models expect values describing the average properties across the layer, while others expect the properties at the border regions. If the number of layers is finite and reduced, then differences in these assumptions can lead to substantial errors. Such type of error primarily originates from the path length the light traces through the atmosphere. To investigate this, the T0A experiment explores how a transit is modelled for the same atmosphere but when considering different numbers of layers (20, 50, 100, 1000). It is demonstrated that if the RT model implements some sort of sub-layering algorithm (e.g., internally dividing the input layers to better capture the change of densities and opacities), these effects can be largely mitigated and practically removed as we can see in Fig. 1 right panels where the transit depth difference between the 1000 layer cases and the lower layer cases tends towards 0 ppm when sub-layering ray-tracing is performed. As noted by Malik et al. (2017), sub-layering can also have a strong effect on the numerical convergence of Radiative Equilibrium (RE) and Radiative Convective Equilibrium (RCE) models, which are used to establish and model the thermal structure of planetary atmospheres.

For the MALBEC intercomparison, the configuration files provide anchor values, values at the border of the layers, and analytical fits to each parameter across the atmosphere, so the modelers can explore and optimize the calculations of the layering as needed. For instance, if the models require layer quantities (not border values), providing straight average values of the two border values would be then more accurate than providing border values.

The T0A is a highly simplified atmosphere designed to test the treatments of layering and sub-layering of the atmosphere, and investigate the effects this has on the final absorbance spectrum. It is based on the atmosphere of GJ1214b, where the temperature decreases linearly from 1000K at 100 bar to 400K at the top of the atmosphere at 1 nano-bar. The experiment explores two scenarios: a case with a simple parameterized water cross-section, increasing with wavelength, and a case employing realistic water opacities. The simple cross-section increases exponentially from $10^{-26}$ cm${}^{2}$/molecule at 0.3 $\mu$m to 8 $10^{-19}$ cm${}^{2}$/molecule at 20 $\mu$m, and is designed to probe different altitudes with wavelength and to highlight the layering effects in the final transmission spectrum. More details about the cross-sections employed in this test can be found at the MALBEC repository. In Figure 1, we show the comparison between the parameterized water absorption (top panels) and using the water cross-section in the bottom panels (\ceH2O; No CIA, no continuum, no Rayleigh for simplicity). This serves to show that our simple opacity covers a dynamic range similar to water.

The effects of not including sub-layering in the RT models can be quite noticeable, as shown in the middle panels of Figure 1, reaching absolute errors beyond 500 ppm at 2 $\mu m$ on the transit radius when no sub-layering is included in the model. Interestingly, the error is not an offset, but more of a trend: the higher in the atmosphere the transit probes, the larger the error. This manifests in smaller errors at weaker water absorption (windows) versus where absorption is stronger (bands). This discretization and ray-tracing error can be largely mitigated when employing a sub-layering algorithm, as shown in the right panels of Figure 1 and computed with PSG, which subdivides each layer in 10 sub-layers when computing the ray-traced column densities.

T3B is a more advanced core experiment, although still with a limited amount of absorbers, but integrating a complex temperature profile, varying molecular abundances and all possible opacity effects (e.g., Rayleigh scattering, CIA, molecular absorptions). In Figure 2, we show preliminary results for this case computed using the models reported in Table 2 and in section 4, in which we consider a base pressure of 1 bar, a planetary diameter of 36318 km, and a stellar diameter of 306108 km (more simulation details available at the MALBEC CKAN repository). One can see a general agreement yet notable differences at certain spectral regions remain. Many of the differences originate from the considered linelists - Assumptions regarding linelists, CIAs and isotopic species can lead to notable differences (Fortney et al., 2020; Niraula et al., 2022). A detailed report exploring this in detail and for the other “casks” will be reported soon following this initial protocol report.

The transit experiments (T0A, T1A, T1B, T3A, T3B, T3C, T4A, T4B) are meant to test the functionalities of the models for the interpretation of transit spectra, which is one of the prime techniques to characterize exoplanets. The resulting spectra in this mode will be primarily defined by the considered opacity processes and the ray-tracing algorithm employed by each model. As such, we have defined a set of experiments to test these models at extreme regimes, highly-compact and dense atmospheres (e.g., TRAPPIST-1 e, Gillon et al. (2017)), extended and light atmospheres (e.g., GJ1214b, Charbonneau et al. (2009)), and hot-Jupiters (e.g., Barstow et al. 2020). By testing transit simulations at different temperatures we expect to better understand the differences due to various line-list databases (e.g., HITRAN vs. HITEMP vs. ExoMol), while the different ambient compositions (e.g., high molecular weight atmospheres vs. \ceH2-rich) would allow us to explore the effects on the assumptions regarding collisional broadening and CIAs. Results from these experiments are primarily tailored to understand the impact of the RT models when interpreting HST, JWST and ARIEL transit data.

As we explore more in detail the origin of the discrepancies of the spectra presented in Figure2, linelists and the methods considered to model opacities by each model are the dominating factors. The case of T3B is interesting in that it simulates a hydrogen-rich atmosphere at relatively high temperatures (350-900K), requiring specializing linelists. On the other hand, models primarily developed to model terrestrial planets and temperate planets are set to operate with HITRAN (mostly air collisional parameters and tailored for $<$300K), which is limited in operability for this specific case. When including the millions of transitions active at these high temperatures, correlated-k and pre-computed high-resolution opacity tables are needed, yet how these databases are computed, the source linelists and the way they are handled by each RT model can lead to substantial differences. As users use these models in their default operational regimes, they should be mindful of the limitations of the linelists configured for the simulation. Detailed comparisons between linelists and opacity modeling methods (e.g., line-by-line, correlated-k, and cross-sections) have been discussed elsewhere as previously discussed, and as part of this model intercomparison, we primarily want to identify spectral regions in which models tend to disagree and why.

In Figure 3, we show a zoom of the near infrared region, in which one can see notable differences between the models, and different levels of transit intensity particularly near the “valley” regions. These valley regions are primarily affected by water and methane opacity effects. In order to demonstrate the effects of linelists in this spectral region of this specific test case, we ran one specific RT model (PSG) but employing different linelists and methods. If measured at extremely high-resolution, the transit spectrum would look substantially spiky as shown with the thin gray trace in panels (b) and (c), yet as we convolve the spectrum to a lower resolution (RP=200), the corresponding spectrum (LBL HITRAN) looks substantially smooth, with small wings and baseline considerations in the line-by-line (LBL) analysis notably affecting the resulting convolved spectrum. The high-energy linelist (HITEMP), which is more suitable for this test case than HITRAN, fills substantially more the valleys though, while the cross sections by Freedman et al. (2014), which rely on a different set of linelists, show a general better agreement to HITEMP than to HITRAN. The line-by-line (LBL) and correlated-k (CK) methods when considering the same source linelist and resolution should in principle agree perfectly, yet small spectral grid interpolation effects do lead to subtle differences as shown for LBL-HITRAN and CK-HITRAN.

A particularly challenging spectral region is the optical and UV, in which linelists normally provide incomplete coverage, and where different RT models treat molecular opacities and Rayleigh scattering in different ways. As such, simulation results differ substantially in this spectral region as the panel (a) in figure 4 shows. Some RT models (e.g., PSG) compute Rayleigh cross-sections as a summation of the individual molecular cross-sections, and based on the individual polarizability of the encompassing molecules. Subtle differences in how each model handles Rayleigh scattering or the assumed molecular polarizability values would lead to differences in the generated spectra. In panel (b), we perturb the PSG model by scaling the polarizability values by 0.9 and 1.1, and one obtains a level of difference comparable to the observed model differences. More complex is the treatment of molecular opacities at these wavelengths, since only sparse and limited molecular spectroscopy is available. At these high frequencies, the energy structure of molecules become increasingly non-quantized into a broad swath of energy levels. As such, classical linelists do not accurately characterize these high-frequency spectral regions, and laboratory measurements of molecular cross-sections are typically needed to fully capture the opacities at these wavelengths. For instance, the line-by-line and correlated-k simulations in PSG are complemented in this spectral range with laboratory cross-sections collected from a range of spectral databases, with most of these cross-sections originating from the MPI-Mainz Spectral Atlas (Keller-Rudek et al., 2013). The shaded regions in panel (c) shows how this complements, which is particularly notable for methane near 0.7 $\mu$m. As better and more complete spectroscopy becomes available over time for a wide range of collisional and excitation regimes, the RT models will need to be updated and adapted to ingest these. More details regarding the databases, the methods, and the line-shape considerations employed by each model, and the associated implications and recommendations, will be presented in a subsequent paper that discusses these findings in detail.

For the direct imaging experiments (D1A, D3B, D5A, D5B), the simulations are to be carried out when the planets are in quadrature (phase 90), the prime phase employed by coronagraphy, and sampling the night-hemisphere (phase 180) to test thermal emissions, and as seen from the star (phase 0, geometric albedo). These modes will allow us to better separate issues when modeling thermal and reflected light, each being susceptible differently to the assumptions regarding multiple-scattering, observational angle, disk discretization and spherical considerations. The simulations are to be carried out employing the highest fidelity available for each model, including multiple-scattering, the highest number of scattering/Legendre terms possible and employing realistic stellar and surface flux models. The tests designed with different levels of complexity. D1A is the simplest test: it has no aerosols, a generally thin atmosphere (with N${}_{2}$, CO${}_{2}$ and H${}_{2}$O starting at 1 atm), and the reflected and thermal emissions are mostly dominated by the modeling of the surface, which is assumed to be simply a Lambert scattering solid surface with a constant albedo and emissivity. Test D3B explores how models handle optically thick atmospheres, in which the deep bottom is defined by opacity functions and not by a solid surface, and it will test how each model handles deep gas atmospheres as those of gas giants. The D3B atmosphere is free of aerosols, so the opacity terms are relatively simple to model and defined by molecular extinction and Rayleigh scattering, yet the treatment of reflected radiation for the quadrature phase will explore how the models discretize and sample the disk across a wide range of illumination regimes in an optically thick multiple scattering non-Lambertian deep atmosphere. The D5A test explores a much more spectrally complex atmosphere of Earth, with strong signatures of H${}_{2}$O, O${}_{3}$, O${}_{2}$, CH${}_{4}$, CO${}_{2}$ and CO, in an ambient air mixture (O${}_{2}$/N${}_{2}$). In particular, this test will be used to analyze the treatment of specialized UV cross-sections and modeling of Rayleigh multiple-scattering across models and RT algorithms. Case D5B is based on D5A but adds the complexity of clouds, and the test is designed to produce a realistic Earth spectrum as targeted by future observatories. One type of cloud is considered in this test, yet the size distribution is assumed to vary across altitude together with its abundance. A specific target of this test is to explore how each model handle handles size distributions of aerosols, optical constants, and multiple-scattering for a wide range of incidence and emission angles.

We expect more differences in the reflected component of the spectrum, in which the effects of multiple-scattering are more notable, and where differences in the treatment of UV opacities and stellar models/templates could become significant. These experiments are specifically tailored to test the differences of RT models in the characterization of large exoplanets with JWST via secondary transit/phase analysis, of gas giants with the Roman Space Telescope, and of Earth-like planets via direct imaging. In particular, future flagship observatories (e.g., Habitable Worlds Observatory, LIFE mission) have as a prime target the characterization of Earth-like planets via direct imaging, so better understanding the current state of RT modeling of these planets is of great importance when designing these future telescopes.