© The Authors 2024

1 Introduction

In the three decades since their first discovery (Wolszczan & Frail 1992; Mayor & Queloz 1995), the inventory of known extra-solar planets has significantly expanded and now counts more than 5700 confirmed objects¹. These vary greatly in their characteristics from small, dense worlds to inflated giant planets, and they encompass a parameter regime beyond the planets found in the Solar System. Techniques used to detect exoplanets vary, but among them the most prolific method to date is the ‘transit method’, where the amount of incident stellar flux received by an observer is reduced by the transit of an exoplanetary companion obscuring part of the stellar disc, which allows the determination of an effective size of the planet.

Studying these transit events provides a unique window into the nature of exoplanetary atmospheres. During the primary transit, the observed stellar flux passes through the atmospheric terminator, where it is attenuated by opacity sources such as molecules or condensates, influencing the observed size of the planet. The wavelength-dependent behaviour of this characteristic ‘depth’ of the transit, which is proportional to the square of the apparent planetary radius in units of the host star’s radius, therefore depends on the nature of the underlying exoplanetary atmosphere. This fact was used in the first successful detections of absorption caused by atomic (Charbonneau et al. 2002; Vidal-Madjar et al. 2003) and molecular (Barman 2007; Tinetti et al. 2007) species in the atmospheres of inflated gas giants, which are called hot Jupiters. Since these initial results, the inventory of characterised exoplanetary atmospheres has increased significantly. The preeminent space-based observatories that have facilitated the characterisation of exoplanet atmospheres are the Hubble Space Telescope (HST) and the Spitzer Space Telescope, which have predominantly provided access to the atmospheres of hot Jupiters. We refer the reader to, for instance, Kreidberg (2018) and Madhusudhan (2019) for comprehensive reviews on the inventory of characterised exo-planet atmospheres.

The launch of the James Webb Space Telescope (JWST; Gardner et al. 2006, 2023) has already started to provide a significant leap forward in our ability to characterise exoplanet atmospheres. One of the first investigations of the atmosphere of a hot Jupiter using the spectroscopic capabilities of JWST has revealed spectral signatures associated with the presence of CO₂ (JWST Transiting Exoplanet Community Early Release Science Team 2023, hereafter referred to as JTERS23). Subsequent observational campaigns have led to the successful atmospheric characterisation of additional hot Jupiters (e.g. August et al. 2023; Bean et al. 2023; Bell et al. 2023; Dyrek et al. 2024). JWST has also improved the possibility of investigating the atmospheres of smaller planets, such as sub-Neptunes (e.g. Kempton et al. 2023; Madhusudhan et al. 2023; Benneke et al. 2024), and detecting atmospheres around terrestrial planets (e.g. Lustig-Yaeger et al. 2023; Moran et al. 2023; Zieba et al. 2023; Kirk et al. 2024). Together with future missions such as Ariel (Tinetti et al. 2018), this leap in accessible exoplanet atmospheres will permit large-scale comparisons of system-level and population-level parameters.

To determine exoplanet atmospheric characteristics, such as the pressure–temperature (p–T) structure and chemical composition based on spectroscopic observations, an inversion process referred to as atmospheric retrieval has found broad application (Madhusudhan & Seager 2009). With this method, the preferred parameter values of an atmospheric forward model are determined through a data-guided statistical sampling process. A variety of atmospheric retrieval codes have been developed that use different parameter exploration methods including NEMESIS (Irwin et al. 2008; Lee et al. 2012), petitRADTRANS (pRT; Mollière et al. 2019; Nasedkin et al. 2024), and TauREx (Waldmann et al. 2015b; Al-Refaie et al. 2021). We refer the reader to MacDonald & Batalha (2023) for a comprehensive overview of atmospheric retrieval framework.

While the statistical sampling of model parameters in atmospheric retrieval is guided by the observational data through Bayesian inference, the reported results of atmospheric characterisation processes also significantly depend on the model parameters themselves. This is exemplified by the varying approaches taken to account for the molecular constituent profiles of atmospheric species in forward models. ‘Free’ chemistry is a heuristic approach, fitting atmospheric volume-mixing ratios (VMRs) individually for each considered species, coarsely determining the presence and abundance values of atmospheric absorbers through a large number of free parameters (e.g. Changeat et al. 2019; Al-Refaie et al. 2022). Self-consistent chemical models rely on an underlying chemical network, and are commonly derived under the assumption of thermochemi-cal equilibrium to infer atmospheric molecular profiles through the minimisation of Gibbs-free energy (e.g. Stock et al. 2018; Woitke et al. 2018; Agúndez et al. 2020). While these models depend on few free parameters, they are also restricted to their model assumptions.

The p–T profile is another important characteristic in the context of exoplanet atmospheres. Measurements within the Solar System have illustrated the complex pressure-temperature structure of atmospheres accessible to in situ investigations (e.g. Seiff et al. 1998; Fulchignoni et al. 2005; Koskinen et al. 2015; Limaye et al. 2017), and they underline the fact that characterisation results for exoplanets strongly depend on the assumptions of simplified models. In atmospheric retrieval, p–T profiles represent another aspect of forward models with a range of underlying assumptions. Analogous to the chemical profiles, a ‘free’ p-T profile retrieves temperature (and potentially associated pressure) values within the atmospheric domain on a purely heuristic basis. The simplest form of this is an isothermal (or one-point) prescription, which has a pressure-independent temperature value as its only free parameter. Multi-point p–T profiles allow for a large degree of flexibility in the determination of the atmospheric temperature structure at the cost of an increasingly large number of free parameters, especially if both the pressure node location and the associated temperature are chosen as free parameters. The other commonly used method of implementing the p–T structure into atmospheric forward models includes parametric p–T profiles (e.g. Madhusudhan & Seager 2009; Guillot 2010). Retrieving information about the p–T structure of an exoplanetary atmosphere is more commonly associated with the analysis of emission spectra, which probe large pressure regimes of the day side of exoplanetary atmospheres. Line et al. (2012, 2013) have investigated the impact of using different parametric and heuristic p–T prescriptions on the retrieval of synthetic emission spectra of hot Jupiters, showing a signal-to-noise-ratio- (S/N) and spectral-resolution-dependent preference for the use of a parametric profile. Blecic et al. (2017) compare the performance of different parametric p–T profiles in retrievals of synthetic hot-Jupiter emission spectra generated from 3D radiative-hydrodynamic simulations.

Compared to this, the retrieval of atmospheric spectra observed during a primary transit are more commonly associated with being sensitive to the characterisation of molecular constituents and the signatures of clouds and hazes. An isothermal prescription for the atmospheric terminator region probed in primary transit observations has found ample application in the analysis of HST data (e.g. Tsiaras et al. 2018; Roudier et al. 2021). However, Rocchetto et al. (2016) show that for a synthetic, cloud-free hot-Jupiter transmission spectrum, under the assumption of JWST-equivalent quality with larger wavelength coverage and increased spectral resolution, the implementation of an over-simplified p–T prescription can significantly bias the retrieved atmospheric VMRs, when compared to a retrieval using a parametric p–T profile. A comprehensive review of outstanding challenges in atmospheric retrievals can be found in, for instance, Barstow & Heng (2020).

In this work, we systematically investigated the influence of p–T prescriptions of increasing complexity on atmospheric retrievals of synthetically generated transmission spectra. We quantify the drawbacks of using an over-simplified, isothermal assumption in connection with the transmission spectrum of a hot Jupiter of the quality produced by JWST. We modelled these spectra after the recent observations of the planet WASP-39 b, accounting for different assumptions about the atmospheric p–T profile, cloud coverage, and transit-depth uncertainty. Furthermore, we determined a feasible level of complexity for the choice of a heuristic p–T profile in the retrieval process, investigating different assumptions about the underlying atmospheric characteristics.

2 Methods

To generate atmospheric forward models and perform atmospheric retrievals, we made use of the open-source atmospheric retrieval framework TauREx (Waldmann et al. 2015b; Al-Refaie et al. 2021), specifically the TauREx3.1 version² (Al-Refaie et al. 2022). TauREx performs radiative transfer calculations for transit or eclipse geometries, as well as atmospheric parameter retrieval as a fully Bayesian parameter inference network. TauREx has found extensive applications in the retrieval of exoplanet atmospheric properties – including the characterisation of transmission (e.g. Tsiaras et al. 2018; Edwards et al. 2021, 2023; Gressier et al. 2022; Saba et al. 2022) and emission (e.g. Changeat et al. 2022; Edwards & Changeat 2024) spectra – encompassing a variety of planetary types ranging from Earth-sized planets to hot Jupiters. To generate exoplanet transmission spectra from atmospheric forward models, we made use of absorption cross-sections from the ExoMol project (Tennyson et al. 2020; Chubb et al. 2021), the HITRAN (Gordon et al. 2022) and HITEMP (Rothman et al. 2010) archives, and the Rayleigh scattering incorporated in TauREx (Cox 2015). We list individually used opacity contributions in Table A.1. To perform parameter estimations, we made use of nested sampling through MultiNest (Feroz et al. 2009; Buchner et al. 2014), with a standardised value for the evidence tolerance of 0.5 and 700 live points (Feroz et al. 2009) in all cases.

To investigate the possibility of recovering more complex pressure-temperature structures from primary transit observations taken by JWST, we performed atmospheric retrievals on a set of synthetic transmission spectra. We modelled our synthetic dataset after the NIRSpec PRISM observation of the hot Jupiter WASP-39 b (PID: 1366, PI: N. Batalha, Co-PIs: J. Bean and K. Stevenson) to mimic the transmission spectrum of an inflated hot Jupiter. In the construction of atmospheric forward models, we used different assumptions about the cloud deck, underlying p–T profile, and associated precision of the transmission spectrum. Using these forward models as observational data, we performed atmospheric retrievals using a range of multi-point p–T profiles with fixed pressure nodes, ranging from an isothermal (or one-point) prescription to an eleven-point profile. The latter represents a self-retrieval validation process, as a fixed-pressure eleven-point profile was used to generate the synthetic input spectra (we refer the reader to Sect. 2.1.1 for more details). In all cases (including the setup of the synthetic transmission spectra), the p–T profile is defined through fixed pressure-temperature nodes. The atmospheric p–T profile is then created using a smoothed linear interpolation on the atmospheric pressure grid, using a standard smoothing window of ten atmospheric layers. We evaluated the preference of different model setups based on the reported Bayesian evidence values and their performance when retrieving the correct underlying atmospheric parameter values.

2.1 Synthetic input data

We generated synthetic atmospheric transmission spectra using the planetary and stellar parameters given in Table 1 to mimic a hot-Jupiter system akin to WASP-39 b. We define the extent of the atmospheric pressure domain used in these forward models as log₁₀(p [bar]) ∈ [1; −9], using 110 layers. As a standardised setup for all cases, we constructed a primary H₂ −He-dominated atmosphere and additionally include H₂O, CO, CO₂, H₂S, and CH₄ as atmospheric absorbers with constant VMRs throughout the atmosphere. Their associated values are given in Table 2. We also include contributions from collision-induced absorption (CIA) of H₂–H₂ and H₂–He (see Table A.1 for references of individual opacity data sources) and Rayleigh scattering (Cox 2015).

We binned the high-resolution forward models generated with TauREx to the resolution of NIRSpec PRISM instrument configuration. In this process, we take into consideration that the full exposure of the NIRSpec PRISM observation of WASP-39 b (after 5 groups) is saturated at wavelengths below approximately 2.2 µm. Therefore, we restricted the synthetic transmission spectra generated in this work to the same wavelength regime to mimic the information content present in the actual observational data (we refer the reader to JTERS23 and Rustamkulov et al. 2023, as well as Schleich et al., in prep. henceforth referred to as Paper II, for more details about the saturated region of this spectrum). The sample of synthetic transmission spectra used as the input dataset is differentiated by three aspects – the underlying pressure-temperature profile (two cases), the cloud-top pressure of a flat-opacity cloud deck (three cases), and the precision of the transit depth values (two cases). In total, this results in twelve distinct atmospheric forward models.

We point out that in the construction of the synthetic transmission spectra, we did not scatter our spectral data points. As we are interested in investigating the potential biases of using simplified pressure-temperature characterisations in the retrieval of atmospheric properties, added randomisation through a self-scattering step of the spectrum could introduce additional, non-controllable biases into the retrieval results. We therefore used the non-scattered instance of the transmission spectrum as an idealised case representing the averaging of results from a larger set of randomised instances (e.g. Feng et al. 2018; Changeat et al. 2019).

2.1.1 Pressure–temperature profiles

We investigated two cases of the underlying pressure-temperature structure as reference cases. These are shown in Fig. 1 and categorised as a less complex profile, represented by a positive temperature-lapse rate (i.e. a negative temperature gradient with increasing altitude, referred to as ‘monotonie’ from here on), and a more complex case characterised by a temperature inversion between a positive temperature-lapse rate in the lower and negative temperature-lapse rates in the upper atmosphere (referred to as ‘inverse’ from here on). Both cases were constructed as eleven-point p–T profiles, where the pressure nodes were log-uniformly distributed in the pressure domain of our atmospheric model. The shapes of the monotonie and inverse profiles are inspired by the cases of the hot Jupiters WASP-76 b and WASP-77 A b, respectively, from Changeat et al. (2022). We used these two profile shapes to investigate the performance of an isothermal prescription, as well as multi-point prescriptions of the p–T profile for two cases of varying complexity. The synthetic transmission spectra generated with the monotonie and inverse input pressure-temperature profile are shown in the top and middle panels of Fig. 2, respectively.

Fig. 1 Constructed pressure-temperature profiles, which are used to create synthetic transmission spectra for our test cases. Temperature (in K) is shown on the x-axis and pressure (in bar) on the y-axis. The (p. T)-pairs (filled markers) are log-uniformly distributed in the pressure regime of the atmosphere. The shapes of the two profiles characterise a simpler, purely positive temperate lapse rate (solid line) and a more complex structure including a temperature inversion point (dashed line). We note that the (p, T)-pair at p = 10−5 bar is the same for both cases.

2.1.2 Atmospheric cloud deck

We constructed the synthetic transmission spectra under the assumption of a flat-opacity cloud deck, defined as an optical depth step-function: $$\tau (\lambda ,z)=\{\begin{array}{ll}\infty \hfill & \text{ if }p(z)\ge p_{\text{cloud }}\hfill \\ 0\hfill & \text{ if }p(z)<p_{\text{cloud }}\hfill \end{array},$$(1)

where τ(λ, z) denotes the optical depth at a certain wavelength, A, and altitude, z; p(z) is the pressure at altitude z, and р_сloud is the designated cloud-top pressure. Sing et al. (2016) show a distribution from clear to cloudy atmospheres in a comparative analysis of the transmission spectra of ten hot Jupiters, where clouds and hazes mute spectral features in the respective atmospheres.

To cover a wide range of possible modulations to the transmission spectrum, we included three different cases of the cloud-top pressure, more specifically log₁₀(p_cloud [bar]) ∈ {0, −3, −5}. The influence of these cloud decks on the final transmission spectrum is illustrated in the first two panels of Fig. 2. In the case of a low-pressure (or high-altitude) cloud deck, the molecular absorption features show significantly lower amplitudes, which are much more pronounced in the case of a high-pressure (or low-altitude) cloud deck.

Fig. 2 Synthetic forward model parameters that influence generated transmission spectra. (Top) Transit depth (in %) against wavelength (in µm) for the monotonie pressure–temperature profile. (Middle) Same, but for the inverse pressure–temperature profile. (Bottom) Associated transit depth uncertainty (in ppm) from a PandExo simulation (solid line) and for the DRS case (dashed line).

2.1.3 Transit depth precision

The third varying parameter we included in our model spectra is the precision of the transmission spectrum, σ_td, which we generated in two different ways. The first one uses error bars generated with the PandExo (Batalha et al. 2017), a noise-simulation tool for observations with JWST. We used standard assumptions for our simulations with PandExo (without an associated noise floor) with the planetary and stellar parameters given in Table 1 and assuming observational parameters corresponding to the NIRSpec PRISM observation of WASP-39 b. As we excluded a noise floor for this case, we see this as a more optimistic assumption on the associated transit depth error bars. We therefore refer to this as the ‘optimistic’ noise case henceforth.

We also included precision values associated with a Eureka!-based reduction of the NIRSpec PRISM dataset of WASP-39 b. A more detailed description of this data reduction is given in Paper II, but the error-bar estimation on the transit depth from this includes an additional point-scatter fitting parameter for the transit depth (e.g. JTERS23; Rustamkulov et al. 2023). We determined the transit depth precision values from this by fitting a fourth-order polynomial to the data reduction results to smooth over outliers and derive a more general error bar trend with respect to wavelength, which is referred to as ‘data-reduction scaled’ (DRS) henceforth. This represents a more conservative approach to determining the error bars on the transit depth, resulting in values larger by a factor of approximately 2.3 when compared to the precision values from PandExo. We included both error-bar cases to investigate and compare the influence of transit-depth S/Ns on the retrieval accuracy of the synthetic data. Both precision estimates considered in the construction of the synthetic transmission spectra are shown in the bottom panel of Fig. 2.

2.2 Retrieval performance metrics

We evaluated the performance of a retrieval with the chosen pressure-temperature characterisation based on two metrics: the Bayesian evidence, E_M, associated with the model solution; and the accuracy of the retrieval with respect to the parameter ground truths.

2.2.1 Bayes factor

We use E_M to calculate the Bayes factor, $$B_{10}=\frac{E_1}{E_0},$$(2)

where E₀ and E₁ refer to the evidence values associated with the underlying models M₀ and M₁, respectively. The Bayes factor serves as a metric to quantify the evidence, based on the underlying data, in favour of a chosen model, M₁, with respect to a reference model, M₀ (or null-hypothesis). We evaluate the natural logarithm of the Bayes factor, ln B₁₀, which is equivalent to Δ ln E, following the prescription from Kass & Raftery (1995) based on individual threshold values for ln B₁₀ given in Table 3.

We note that the Bayes factor is a relative metric comparing the fitting performance between different model assumptions, but this does not provide any conclusion regarding the absolute fitting performance of any single model. We omit an absolute metric from this work. The construction of the synthetic retrievals performed here tests for potentially inherent biases stemming from assumptions about the atmospheric p–T profile, but it does not include any noise apart from assumed observational uncertainties. This implies that the fits performed here can reproduce the input data (meaning the transmission spectrum) to an indistinguishable degree, irrespective of the assumed multi-point p–T profile. We investigated how, for oversimplified p–T characterisations, this leads to compensatory deviations in other model parameters. We illustrate an overview of absolute fit performances in Appendix B, but we summarise this as follows: the residuals of retrievals performed with any of the multi-point profiles, apart from the isothermal profile, are indistinguishable from each other. While we see a minor excess in residuals larger than the 1σ-equivalent error-bar size in retrievals performed using the isothermal profile, 99% of them fall below the associated 2σ-value, and we see no individual residual with an absolute value larger than 3σ.

2.2.2 Retrieval accuracy

Additionally, we evaluated the performance of each retrieval against the known input parameter values of the underlying synthetic spectra. To evaluate the accuracy of the retrieved parameters, we calculate the distance between the estimated parameter (θ_p,est) and the known true value (θ_p,GT) with $$\text{Δ}{\theta }_{\text{p}}={\theta }_{\text{p},\text{est}}-{\theta }_{\text{p},\text{GT}}$$(3)

for the planetary reference radius, R_p, the VMRs given in Table 2, and the cloud-top pressure, p_cloud, corresponding to the individual synthetic forward model cases described in Sect. 2.1. This value will be positive if the retrieved parameter value overestimates the underlying true value and negative if the underlying true values are underestimated.

We scale Δθ_p using a factor representing the attributed uncertainty in the estimated parameter value, employing two different prescriptions. Firstly, we calculated a centred credible interval (CCI) containing a 3σ equivalent of the sample parameter population. For this, we determine the centre of the credible interval to be the median of the posterior distribution (or the 50% quantile) and the edges to be defined by the quantile values encompassing 99.7% of the posterior distribution. Additionally, we use the 1σ value given by TauREx as a scale factor, which is computed as the CCI encompassing 68% of the posterior distribution. An example of the parameter estimation evaluation is shown in Fig. D.1 for the case of a synthetic transmission spectrum generated with an underlying inverse p–T profile, a cloud-top pressure of 1 bar, and DRS transit-depth error bars. The asymmetry of the upper and lower boundaries around the median value of the retrieved parameter derives from the underlying non-Gaussianity of the marginalised posterior distribution. We provide a more detailed view of all accuracy values associated with individual parameters in the supplementary repository for this work³. In the case of the 3σ-equivalent CCI, we evaluate the retrieval accuracy on a binary scale, with the true parameter value being within the credible interval or not. For the 1σ intervals provided by TauREx, we compute the scaled parameter distance outright.

Fig. 3 Transmittance for synthetic spectra generated with monotonic (left panel) and inverse (right panel) input pressure-temperature profiles. Both panels show wavelength (in µm on the x-axis) and pressure (in bar) on the y-axis. For illustrative purposes, the transmittance in these plots is calculated without the cloud-deck cases. The cloud-top pressure levels, pcloud, considered in this work (as described in Sect. 2.1.2) are marked with a solid (pcloud = 1 bar), dashed (pcloud = 10−3 bar), and dotted (pcloud = 10−5 bar) black line, respectively.

3 Results and discussion

To quantify the performance of different pressure-temperature prescriptions on the recovery of atmospheric characteristics, we performed a homogenised set of retrievals on the sample of synthetic transmission spectra using a variety of multipoint pressure–temperature profiles. As a reference case, we performed retrievals using a free p–T profile with eleven pressure nodes, which represents a self-retrieval on the constructed input cases described in Sect. 2.1.1. Firstly, we compared the performance of this to retrievals that utilise an isothermal prescription, equivalent to fitting one temperature parameter. We then iteratively increased the number of pressure nodes, performing retrievals using a two-point profile (which represents the assumption that a temperature gradient in the atmosphere can be retrieved), a four-point profile (allowing for the retrieval of a potential temperature-inversion point), as well as a six-point and an eight-point profile (to investigate the potential of recovering an even more detailed p–T structure than with the two- or four-point profile).

3.1 Retrieving synthetic spectra

In all cases, we used TauREx with the general setup described in Sect. 2, and we used the same forward model setup as in the construction of the transmission spectra; that is, the pressure domain, assumption of homogeneous VMRs, and atmospheric opacity sources. For all cases, we retrieved the planetary radius at the bottom of the atmospheric pressure regime, R_p; molecular VMRs, log₁₀(X_VMR); and cloud-top pressure, log₁₀(p_cloud); representing seven free parameters shared by all models. In addition, for each model we fitted the free parameters associated with the p–T prescription. In the case of an isothermal profile, this means retrieving the temperature parameter T_iso, while in the case of the multi-point profiles, we retrieve the temperature values, T_n, at fixed pressure nodes, p_n, within the atmosphere. These pressure nodes are set to the top of the atmosphere (TOA) and bottom of the atmosphere (BOA) for the first and last node of each case, respectively. The remaining nodes are distributed close to log-uniform within the pressure domain – we separated them as equally as possible between the top and bottom of the atmospheric pressure domain, but make an informed choice to slightly cluster them around the expected probed atmospheric pressure region. This is illustrated in Fig. 3, which shows the transmittance of the atmospheric forward model for both the spectra generated with the inverse, as well as with the monotonic pressure-temperature profile, calculated without the respective cloud-deck cases. The probed atmospheric pressure regime, excluding the cloud deck, extends approximately between 1 mbar and 10 µbar. The explicit values of log₁₀(p_n [bar]) for these multi-point profiles are as follows:

• {1, −3, −7, −9} for the four-point profile

• {1, −1, −3, −5, −7, −9} for the six-point profile

• {1, −1, −2, −3, −4, −5, −7, −9} for the eight-profile

We note that the pressure nodes for the self-retrieval performed with an eleven-point p–T profile are the same ones shown in Fig. 1. In total, this results in a number of free parameters between eight (isothermal case) and eighteen (eleven-point pressure-temperature profile). We list these retrieved parameters, and their associated prior distributions, in Table 4.

We illustrate the evaluation of models based on the Bayes factor, ln B₁₀, in Fig. 4, showing the Bayes factor in reference to the self-retrieval using the eleven-point profile, as a function of the number of the temperature points in each model (exemplified on the DRS noise case). In the majority of cases, the isothermal model is strongly disfavoured against all other multi-point models, illustrating the fact that the isothermal prescription of the p–T profile is not capable of accurately capturing the complexity of either of the underlying pressure–temperature profiles shown in Fig. 1. The exception to this is the synthetic transmission spectra generated with a low-pressure cloud top (p_cloud = 10⁻⁵bar). In these cases, the Bayes factor does not sufficiently distinguish between any of the models, owing to the fact that a flat-opacity cloud deck at pressure levels this low obscures significant fractions of the molecular signatures, masking the impact of the p–T profile on the molecular contributions. This can also be seen in Fig. 3, where the low-pressure cloud deck as a flat-opacity step function obscures a significant fraction of the wavelength-dependent transmittance at higher pressures, effectively masking the probed pressure regime.

Apart from the disfavoured performance of the isothermal p–T profile, the behaviour of the Bayes factor is clearly separated between the cases of an underlying monotonic p–T profile (left panel of Fig. 4) and an underlying inverse p–T profile (right panel of Fig. 4). In the former case, the two- and four-point profiles are strongly to moderately preferred, while the Bayes factors for profiles with a higher degree of complexity converge towards zero. However, in the inverse case, the two-point profile is no longer preferred against the self-retrieval, while there is still a moderate preference for the four-point profile in the high-pressure cloud case. The model performance evaluation for the optimistic noise case, shown in Fig. C.1, provides conclusions in line with the DRS case, with the distinction that the smaller transit-depth error bars associated with the optimistic case lead to more distinctly differentiated Bayes factors for each model. This is in line with the expectation that the increased S/N of the transmission spectrum in these cases provides higher information content. The optimistic noise case provides a clearer preference for the two-point and four-point profiles in the monotonic input p–T case and a clearer Bayes factor peak for the four-point profile used in the retrievals on spectra generated with the inverse p–T profile.

From the evaluation of ln B₁₀, we can conclude that the isothermal profile should be disregarded compared to more complex p–T characterisation when retrieving transmission spectra of a quality in line with what is expected from JWST observations. This is in agreement with the conclusions presented in Rocchetto et al. (2016), which showed that in retrievals of transmission spectra modelled for expected JWST observations of hot Jupiters, assuming an isothermal p–T profile leads to a significant bias in the retrieved atmospheric VMRs.

The accuracy of the retrieved parameters corroborates the preference for a more complex prescription of the pressure-temperature profile. Retrievals performed using an isothermal profile consistently lead to biased parameters, in several cases with misleadingly high precision. We illustrate this in Fig. 5, which shows the retrieval accuracy evaluated using the 3σ-equivalent CCI for synthetic transmission spectra containing a high-pressure cloud-top (for both assumed noise cases). In all cases, performing retrievals with the isothermal prescription fails to retrieve the VMRs of the main atmospheric molecular species. The exception to this (in the monotonic case) is CH₄, which, due to the low underlying abundance value and lack of distinct absorption signatures, has significantly broader posterior distributions. Following from this, the associated parameter estimation uncertainty values are significantly larger, and they more clearly encompass the underlying true value. In the case of the other retrieved molecular VMRs, the value of Δθ_p for the isothermal case is centred at around twice the size of the 3σ CCI.

When comparing the accuracy of the retrieved parameters between retrievals performed on synthetic spectra generated under the DRS noise assumption (as exemplified by the top two panels in Fig. 5) and on spectra generated using the optimistic noise assumption (as exemplified by the bottom two panels in Fig. 5), it is also apparent that in the latter case, the values of Δθ_p are, on the scale of the 3σ CCI, larger than for the equivalent DRS case. This is influenced by the higher reported retrieval precision, reflected in a smaller centred credible interval. We provide a complete overview of the accuracy evaluation for all parameters and cases, including a comparison between the credible intervals of the retrieved parameters, in the supplementary material to this work⁴.

As suggested by the Bayes factor analysis for model preference, distinguishing between favoured multi-point characterisation for the pressure–temperature profile is less clear. In the case of a monotonic underlying p–T structure, we find equal performance for all multi-point profiles tested, regardless of the accuracy scale factor used in evaluating them. The pressure-temperature structure of this profile, shown in Fig. 1, suggests that this is to be expected as a simple atmospheric temperature gradient, produced by the two-point profile, can be judged as sufficient to reproduce the overall structure of the atmospheric p–T profile. Considering either approach in evaluating the retrieval accuracy shows the equivalently successful performance of all multi-point profiles in retrieving the correct true-parameter values. Based on the Bayes factor evaluations in Figs. 4 and C.1, the two-point and four-point profiles are the most favoured ones; the two-point profile could be preferred overall due to it having fewer free parameters.

This picture changes in the case of the inverse input p–T profile. When comparing the structure of this underlying profile (see also Fig. 1) with the expected p–T structure in retrieving a two-point profile, it is apparent that the two-point profile does not provide enough complexity to capture the temperature inversion point. This leads to wrongly retrieving the planetary reference radius and cloud-top pressure in several cases here, one of which is illustrated in the left panels of Fig. 5. Distinguishing the performance of retrievals using a pressure–temperature profile with more than two points based on the accuracy alone proves to be more inconclusive. As would be expected for models with an increasing number of free parameters, we see no more distinguishable difference in the retrieved parameter accuracy for any of these cases.

The connection between the preference of different multipoint profiles to the underlying pressure-temperature structure of the input spectrum is further illustrated in Fig. 6, which shows the retrieved p–T profiles of the isothermal, two-point, and four-point assumption in the context of the DRS noise case. With an underlying inverse p–T profile, a multi-point profile with more than two pressure nodes is necessary to capture the input p–T structure. A two-point profile, which can only produce a temperature gradient under the assumptions taken here, cannot sufficiently reproduce the inversion point of this profile. On the other hand, in the case of the monotonic input case, it becomes quite apparent that a two-point profile is already sufficient in capturing the (at first order) gradient-like structure of the input p–T profile. In either case, the isothermal profile is, by its construction, unable to reproduce any vertical temperature structure in the atmosphere, and it fails to accurately retrieve the majority of the main atmospheric characteristics. This can also be seen in the bottom panels of Fig. 6; they show that, apart from the isothermal profile, each of the retrieved multi-point profiles is in good agreement with the underlying p–T structure used as the input for the corresponding synthetic transmission spectrum. We show the retrieved pressure-temperature profiles for the synthetic cases associated with the optimistic noise case in Appendix E.

We note that in all retrievals run under an isothermal assumption, a distinct difference between the accuracy of the retrieved atmospheric parameters is visible when comparing the input cases using the inverse p–T profile (illustrated in the two left panels of Fig. 5) and the monotonic p–T profile (illustrated in the two right panels of Fig. 5). With an underlying inverse profile, the value of Δθ_p for all molecular VMRs is consistently negative, implying that the parameters are underestimated during the retrieval process. In the cases with an underlying monotonic p–T profile, these values are consistently positive, implying that the parameters are overestimated. As shown in Fig. 3, the probed pressure regime in both cases extends approximately between 1 mbar and 10 µbar. Comparing this to the retrieved p–T profile for the isothermal case (Figs. 6 and E.1), it is apparent that in the monotonic case the temperature profile is lower than the underlying ground truth, while in the inverse case, the opposite is true. This anti-correlation between the biases in temperature and atmospheric VMRs is attributable to the influence of the atmospheric scale-height in the transit geometry. If the temperature is overestimated, the scale height of the atmosphere increases. Subsequently, this also increases the path length within the atmospheric layer. In the calculation of the optical depth, the increased path-length is compensated by lowering the chemical abundances (see Eq. (5) in Al-Refaie et al. 2021, represented by the column density). If the temperature is underestimated, the opposite effect is apparent.

Fig. 4 Bayesian evidence comparison for retrievals in all cases of the synthetic dataset associated with the DRS noise case, split into the synthetic spectra generated using the monotonic input p–T profile (left) and the inverse input p–T profile (right). Both panels show the number of retrieved temperature points (corresponding the implemented multi-point profile) on the x-axis and the model-associated Bayes factor on the y-axis. Data points are separated into spectra generated with high (pcloud = 1 bar, circles), medium (pcloud = 10−3 bar, squares), and low (pcloud = 10−5 bar, diamonds) cloud-top pressure input values. The grey shaded area denotes the Bayes factor threshold given in Table 3. The black cross denotes the reference model in all cases (eleven-point p–T profile), and Bayes factor values fulfilling |ln B10| > 6.0 are denoted by a black arrow.

Fig. 5 Accuracy of retrieved parameters, showing individual parameters on the x-axis and accuracy calculated based on the 3σ-equivalent centred credible interval on the y-axis. The top two panels represent retrievals performed on a synthetic spectrum generated under the DRS noise assumption, while the bottom two panels show the same for the optimistic noise case. (Left) Results for retrievals performed on a synthetic spectrum generated with the inverse input p–T profile. (Right) Same, but for a synthetic spectrum generated with the monotonic p–T profile. In all panels, the reported accuracy values are colour-coded by the p-T profile used in the retrieval. The black dashed lines denote the boundaries of the 3σ CCI, and associated accuracy values of |Δθp| > 5.0 are cut off for visualisation purposes.

Fig. 6 Comparison between retrieved pressure–temperature profiles and corresponding ground truths, showing temperature (in K) on the x-axis and pressure (in bar) on the y-axis, for synthetic cases using the error bar size scaled from data-reduction results. The first three panels show the inverse p–T profile as the ground truth, and are ordered by increasing cloud-top pressure: log10(pcloud [bar]) = {0, −3, −5). The second three panels show the same, but for the monotonic p–T profile as the ground truth.

3.2 Outlook

We note that the synthetic retrieval investigation described above only represents a preliminary analysis the performance of an isothermal prescription against that of more complex multipoint profiles. In this work, we chose the underlying pressure– temperature profiles with arbitrary complexity to cover a simpler, positive temperature-lapse rate and a more complex structure including a temperature-inversion point, respectively. Similarly, we modelled the chemical structure of the mock atmosphere after previously described results for the hot Jupiter WASP-39 b (JTERS23) and did not incorporate physics-informed choices of the atmospheric chemical structure based on the chosen p–T profile. Additionally, in this work we only looked at the case of a primordial, H₂ / He-dominated atmosphere, which, due to their low mean molecular weight, are prime targets for detailed atmospheric characterisation efforts. A more comprehensive study aiming to address this question for a broader variety of p–T structure and planetary characteristics would complement the conclusions drawn here.

In this work, we constructed an input transmission spectrum using VMRs that are constant with altitude. This guarantees that, when performing atmospheric retrievals under the same assumption of constant profiles, we do not introduce additional biases in oversimplifying the vertical VMR structure during the retrieval. However, the assumption about constant VMRs might be considered an oversimplification, along the same lines as an isothermal p–T profile. As shown in Fig. 5, there is a degeneracy between the retrieved temperature and molecular VMR values. Through their influence on the atmospheric scale height, the biases in temperature and VMR values are anti-correlated. This implies that an oversimplified chemical profile used for atmospheric retrievals could induce a similar bias to the retrieved pressure– temperature structure. Changeat et al. (2019) showed that, in the context of simulated observations of JWST and Ariel, using constant chemical profiles when the underlying true chemical profile is more complex introduces this degenerate bias in retrieving an isothermal p–T profile. While a test of this manner is beyond the scope of this work, additional investigations to catalogue the effects of these sources of biases in atmospheric retrievals will help inform best practices for the characterisation of exoplanet atmospheres in a consistent manner.

We also point out that, in the case of multi-point pressure– temperature profiles, the location of the pressure nodes influences the performance of the associated retrieval, especially with respect to reproducing ground-truth p–T profiles, as is the case here. Fixing the pressure nodes inherently includes the possibility that a significant “feature” of the underlying p–T structure (such as the inversion point of the inverse profile tested here) is missed. This could be counteracted by including the pressure-nodes of the profile as free parameters as well. However, this effectively doubles the number of free parameters associated with any multi-point profile, significantly increasing the computational requirements of running such retrievals (Changeat et al. 2021). Another approach to mitigate this problem could be provided by determining a pressure-region of maximum contribution function to the optical depth (or transmittance) of the atmosphere and clustering any pressure-nodes around this relevant region. Waldmann et al. (2015a) described an approach of this manner for the retrieval of exoplanetary emission spectra, performing a multi-stage approach to dynamically adjust the complexity of the retrieved p–T profile. A detailed study of the application of this to transmission spectra is outside the scope of this work, but it would provide additional insight into best practices for determining a free p–T structure of exoplanet atmospheres from transmission spectroscopic observations.

Applying atmospheric retrievals to real datasets adds observational noise characteristics as a potential additional source of bias in the atmospheric characterisation process. A recent, in-depth analysis of the information content of JWST transit observations of WASP-39 b by Lueber et al. (2024) has explored the question of p–T complexity in the case of real data. By performing a Bayes factor analysis on atmospheric retrievals of the individual near-infrared transmission spectra taken by JWST, the authors found a preference for isothermal prescriptions in several cases. Based on the results shown here, we advocate for implementing p–T prescriptions of at least the complexity of a two-point profile. If atmospheric retrievals using more complex p–T prescriptions show disagreement in the retrieved chemical profiles when compared to retrievals performed using an isothermal prescription, it would indicate an unresolved underlying bias. We present one possible source of this, although other sources of model biases could also impact transmission results. For instance, oversimplified chemical profiles (Changeat et al. 2019), 3D effects (e.g. Espinoza & Jones 2021; Pluriel et al. 2022), and poorly constrained planetary masses (e.g. Changeat et al. 2020; Di Maio et al. 2023) have already been investigated in the past. To quantify the importance of such biases, retrieval results should systematically be evaluated against complex prescriptions in controlled scenarios.

4 Conclusions

We investigated the recoverable complexity of pressure-temperature profiles associated with transmission spectra observed with JWST based on synthetic transmission spectra modelled after the hot Jupiter WASP-39 b. We investigated several scenarios, including two different underlying pressure– temperature structures, several different cloud-top pressure values, and two different assumptions about the associated transit depth uncertainty to cover a large range of possible input transmission spectra. We tested the performance of atmospheric retrievals using isothermal and multi-point pressure–temperature profiles with two, four, six, and eight fixed-pressure nodes. We evaluated the performance of these individual retrievals based on two metrics: the Bayes factor and the accuracy of the retrieved parameters in comparison with the known true values. The Bayes factor was evaluated in reference to retrievals performed with an eleven-point pressure–temperature prescription, which represents a self-retrieval of the input p–T profile.

Based on these two metrics, we find that the isothermal profile is an insufficient approximation for atmospheric retrievals of transmission spectra with JWST-level expected precision. We find that the Bayes factor analysis consistently rejects retrievals run with the isothermal prescription, and we see that these retrievals consistently lead to wrongly recovered atmospheric parameters. In comparing the performance of the individual multi-point profiles that are more complex than the isothermal case, the Bayes factor analysis is able to distinguish preferred profiles based on the underlying input p–T case. Both the two- and four-point profiles show a moderate-to-strong preference when the underlying pressure–temperature profile shows a purely positive lapse rate. In the case of an underlying p–T profile showing an inversion point, we find that the four-point prescription is the only preferred profile. Compared to this, more complex multi-point profiles are not preferred in any cases, most likely owing to the fact that their additional free parameters penalise the Bayesian evidence evaluation. In comparing the accuracy of the retrieved parameters, we sporadically find incorrect values for the two-point profile, but see no clear distinguishing factor for the more complex prescriptions.

This suggests that, among the p–T characterisations investigated in this work, the four-point profile is the least complex model that correctly retrieves the known input parameters in all cases and leaves flexibility in the retrieval to capture potentially more complex pressure-temperature structures. While a finegrained selection process of an optimal multi-point p–T profile based on the expected nature of the transmission spectrum could still be improved in future work, we recommend not using an isothermal prescription as it has been shown to wrongly bias the retrieved atmospheric VMRs. This bias also depends on the underlying p–T structure being positive in the case of a positive temperature lapse rate, but negative in the presence of a temperature inversion point around the probed atmospheric pressure regime.

Data availability

We provide an online repository with all data products used in this work, as well as supplementary data products and plots, at https://doi.org/10.5281/zenodo.10497342.

Acknowledgements

We thank the anonymous referee for their helpful comments and suggestions. S. Schleich extends his gratitude to G. Van Looveren, N. Pawellek, and M. Voyer for insightful discussions on the representation of the results. This project was funded by the FGGA Emerging Field Grant 2021. We acknowledge financial support by the University of Vienna and Österreichis-che Forschungsgemeinschaft (ÖFG). Parts of this work were achieved using the Vienna Scientific Cluster (VSC). We gratefully acknowledge the use of open-source libraries for the Python programming language: matplotlib (Hunter 2007), numpy (Harris et al. 2020), and pandas (McKinney 2010).

Appendix A Atmospheric opacity data

Table A.1 shows an overview of the opacity data and associated references used within the retrievals conducted in this work.

Appendix B Absolute retrieval fit performance

In evaluating the atmospheric retrievals of the synthetic transmission spectra described in Sect. 2.1, we use the Bayes factor, ln B₁₀, as a relative metric to compare the performance of individual model. As the Bayes factor is a relative metric, comparing the evidence value, ln E, associated with each model, it makes no comment about the absolute fit performance for each run. To review to total performance of each model, differentiated only by the characterisation of the p–T profile, we investigate the total distribution of individual fit residuals, δ_i, scaled to the associated data point error bars, $${\delta }_i=\frac{X_i-O_i}{{\sigma }_i},$$(B.1)

where for each data point of the synthetic transmission spectrum, O_i, we calculate the residual associated with the best-fit model point, X_i, and scale it by the error-bar size, σ_i, associated with the data point. Figure B.1 (top panel) illustrates this distribution for each pressure–temperature profile used in our retrieval runs, where the residuals for each of the twelve synthetic input spectra are combined. The distributions of residuals for the two-, four-, six-, eight-, and eleven-point profiles are indistinguishable from each other. In comparison to that, the residuals for all retrievals using the isothermal profile show a broader wings towards larger values. However, we still find that, in a cumulative count of all residual values from isothermal retrievals, 99.2% fall within |δ| < 2.0, and 95% fall within |δ| < 1.0 (bottom panel of Fig. B.1). This is attributable to the way in which the synthetic spectra in our work are set up – Our investigation aims to find inherent biases stemming from the utilisation of different p–T profile complexities, which is why we omit this scatter process to not introduce uncontrollable biases (Feng et al. 2018; Changeat et al. 2019). As we do not add any scatter do the data points in the construction of the transmission spectra, we enable the fitting process to very closely reproduce the shape of the underlying cross-sectional data used to generate the spectra. Based on this, we deem it not necessary to discuss the absolute performance of each fitting process in the retrievals we conduct, and restrict our analysis to the relative model comparison through the Bayes factor, as well as the accuracy evaluation in the known underlying input parameters of each synthetic transmission spectrum.

Fig. B.1 Distributions of retrieval residuals (scaled to the observational error bar, σobs), counted in bins with a width of 0.25 σobs. (Top) Total number of counts per bin. (Bottom) Cumulative fraction of residuals. Both panels show the collected distribution of residuals from retrievals applied to each of the twelve synthetic transmission spectra, but counted separately for each pressure–temperature prescription.

Appendix C Additional Bayes factor analysis

In addition to the results of synthetic retrievals on the mock atmospheric transmission spectra created with error bars scaled by data reduction results, Fig. C.1 shows the same analysis on spectra created with error bars generated with PandExo.

Fig. C.1 Same as Fig. 4, but for transmission spectra generated using the transit-depth error bars resulting from the optimistic noise case.

Appendix D Retrieval parameter accuracy

Supplementary to the parameter accuracy evaluation discussed in Sect. 2.2.2, we show an example of the 3-sigma-equivalent centred credible interval and the 1σ parameter uncertainty reported by TauREx in Fig. D.1. We refer to the supplementary repository of this work for a full set of all individual accuracy evaluations⁵.

Fig. D.1 Accuracy of retrieved parameters, in comparison to known input value, for the synthetic transmission spectrum generated with the inverse p–T profile, a cloud-top pressure of 1 bar, and the DRS noise case. In all panels, the parameter value is shown in the x-axis, and varying retrieval cases are listed on the y-axis. For each retrieval case, the median retrieved value (black dot), as well as the 1σ uncertainty provided by TauREx (black error bar), and the boundaries of the 3σ-equivalent CCI (grey error bars) are shown. The known true parameter value is indicated with a red dashed line.

Appendix E Additional retrieved p–T profiles

Complementary to Fig. 6, we show the retrieved isothermal, two-point, and four-point p–T profiles for all input transmission spectra generated with the ‘Pandexo’ noise case in Fig. E.1, together with the underlying input pressure–temperature profiles, and the pressure-layer of the associated flat-opacity cloud deck.

Fig. E.1 Same as Fig. 6, but for the synthetic spectra generated using the ‘Pandexo’ noise case.

References

All Tables

All Figures

Fig. 1 Constructed pressure-temperature profiles, which are used to create synthetic transmission spectra for our test cases. Temperature (in K) is shown on the x-axis and pressure (in bar) on the y-axis. The (p. T)-pairs (filled markers) are log-uniformly distributed in the pressure regime of the atmosphere. The shapes of the two profiles characterise a simpler, purely positive temperate lapse rate (solid line) and a more complex structure including a temperature inversion point (dashed line). We note that the (p, T)-pair at p = 10−5 bar is the same for both cases.

In the text

	Fig. 2 Synthetic forward model parameters that influence generated transmission spectra. (Top) Transit depth (in %) against wavelength (in µm) for the monotonie pressure–temperature profile. (Middle) Same, but for the inverse pressure–temperature profile. (Bottom) Associated transit depth uncertainty (in ppm) from a PandExo simulation (solid line) and for the DRS case (dashed line).
In the text

Fig. 3 Transmittance for synthetic spectra generated with monotonic (left panel) and inverse (right panel) input pressure-temperature profiles. Both panels show wavelength (in µm on the x-axis) and pressure (in bar) on the y-axis. For illustrative purposes, the transmittance in these plots is calculated without the cloud-deck cases. The cloud-top pressure levels, pcloud, considered in this work (as described in Sect. 2.1.2) are marked with a solid (pcloud = 1 bar), dashed (pcloud = 10−3 bar), and dotted (pcloud = 10−5 bar) black line, respectively.

In the text

Fig. 4 Bayesian evidence comparison for retrievals in all cases of the synthetic dataset associated with the DRS noise case, split into the synthetic spectra generated using the monotonic input p–T profile (left) and the inverse input p–T profile (right). Both panels show the number of retrieved temperature points (corresponding the implemented multi-point profile) on the x-axis and the model-associated Bayes factor on the y-axis. Data points are separated into spectra generated with high (pcloud = 1 bar, circles), medium (pcloud = 10−3 bar, squares), and low (pcloud = 10−5 bar, diamonds) cloud-top pressure input values. The grey shaded area denotes the Bayes factor threshold given in Table 3. The black cross denotes the reference model in all cases (eleven-point p–T profile), and Bayes factor values fulfilling |ln B10| > 6.0 are denoted by a black arrow.

In the text

Fig. 5 Accuracy of retrieved parameters, showing individual parameters on the x-axis and accuracy calculated based on the 3σ-equivalent centred credible interval on the y-axis. The top two panels represent retrievals performed on a synthetic spectrum generated under the DRS noise assumption, while the bottom two panels show the same for the optimistic noise case. (Left) Results for retrievals performed on a synthetic spectrum generated with the inverse input p–T profile. (Right) Same, but for a synthetic spectrum generated with the monotonic p–T profile. In all panels, the reported accuracy values are colour-coded by the p-T profile used in the retrieval. The black dashed lines denote the boundaries of the 3σ CCI, and associated accuracy values of |Δθp| > 5.0 are cut off for visualisation purposes.

In the text

Fig. 6 Comparison between retrieved pressure–temperature profiles and corresponding ground truths, showing temperature (in K) on the x-axis and pressure (in bar) on the y-axis, for synthetic cases using the error bar size scaled from data-reduction results. The first three panels show the inverse p–T profile as the ground truth, and are ordered by increasing cloud-top pressure: log10(pcloud [bar]) = {0, −3, −5). The second three panels show the same, but for the monotonic p–T profile as the ground truth.

In the text

Fig. B.1 Distributions of retrieval residuals (scaled to the observational error bar, σobs), counted in bins with a width of 0.25 σobs. (Top) Total number of counts per bin. (Bottom) Cumulative fraction of residuals. Both panels show the collected distribution of residuals from retrievals applied to each of the twelve synthetic transmission spectra, but counted separately for each pressure–temperature prescription.

In the text

	Fig. C.1 Same as Fig. 4, but for transmission spectra generated using the transit-depth error bars resulting from the optimistic noise case.
In the text

Fig. D.1 Accuracy of retrieved parameters, in comparison to known input value, for the synthetic transmission spectrum generated with the inverse p–T profile, a cloud-top pressure of 1 bar, and the DRS noise case. In all panels, the parameter value is shown in the x-axis, and varying retrieval cases are listed on the y-axis. For each retrieval case, the median retrieved value (black dot), as well as the 1σ uncertainty provided by TauREx (black error bar), and the boundaries of the 3σ-equivalent CCI (grey error bars) are shown. The known true parameter value is indicated with a red dashed line.

In the text

	Fig. E.1 Same as Fig. 6, but for the synthetic spectra generated using the ‘Pandexo’ noise case.
In the text