{'fraction_non_alphanumeric': 0.02547604556125554, 'fraction_numerical': 0.017650639074863056, 'mean_word_length': 5.064460848932244, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 0, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Hierarchical Bayesian Atmospheric Retrieval Modeling for Population Studies of Exoplanet Atmospheres: A Case Study on the Habitable Zone

With the growing number of spectroscopic observations and observational platforms capable of exoplanet atmospheric characterization, there is a growing need for analysis techniques that can distill information about a large population of exoplanets into a coherent picture of atmospheric trends expressed within the statistical sample. In this work, we develop a Hierarchical Bayesian Atmospheric Retrieval (HBAR) model to infer population-level trends in exoplanet atmospheric characteristics. We demonstrate HBAR on the case of inferring a trend in atmospheric CO2 with incident stellar flux, predicted by the presence of a functioning carbonate-silicate weathering negative feedback cycle, an assumption upon which all calculations of the habitable zone (HZ) rest. Using simulated transmission spectra and JWST-quality observations of rocky planets with H2O, CO2, and N2 bearing atmospheres, we find that the predicted trend in CO2 causes subtle differences in the spectra of order 10 ppm in the 1-5 um range, underscoring the challenge inherent to testing this hypothesis. In the limit of highly precise data (100 stacked transits per planet), we show that our HBAR model is capable of inferring the population-level parameters that characterize the trend in CO2, and we demonstrate that the null hypothesis and other simpler trends can be rejected at high confidence. Although we find that this specific empirical test of the HZ may be prohibitively challenging in the JWST era, the HBAR framework developed in this work may find a more immediate usage for the analysis of gas giant spectra observed with JWST, Ariel, and other upcoming missions.


INTRODUCTION
The habitable zone (HZ) provides a tangible starting point in the search for habitable exoplanet surface environments and life beyond the solar system (Kasting et al. 1993;Kopparapu et al. 2013;Kaltenegger 2017;Meadows & Barnes 2018). Despite the HZ's scientific lineage rooted in Earth system science (e.g., Hart 1978Hart , 1979, understanding the persistent habitability of Earth over geologic time remains an ongoing interdisciplinary investigation (e.g., Goldblatt & Zahnle 2011;Charnay et al. via atmospheric CO 2 buffering (Berner 2003), but which is also an assumed ingredient in habitable zone calculations (Kasting et al. 1993;Williams & Kasting 1997;Kopparapu et al. 2013). In the carbonate-silicate weathering cycle, atmospheric CO 2 warms the climate via the greenhouse effect. An increase in volcanic outgassing produces more CO 2 , which further warms the climate. However, the weathering rate of continents also increases with temperature, thereby increasing the rate at which CO 2 is removed from the atmosphere and ultimately subducted back into the mantel. Thus, the temperature dependence of the weathering rate provides the climate-stabilizing negative feedback that helps to maintain habitable surface temperatures against increases in CO 2 and the solar luminosity over geologic timescales (Glaser et al. 2020).
The link between Earth's long-term climate evolution and our perspective on exoplanet habitability provides a compelling opportunity to observationally test such hypotheses on the nature of planetary habitability using the population of exoplanets. Bean et al. (2017) outlined a statistical comparative planetology approach to empirically test the habitable zone hypothesis by recognizing that the climate model calculations for the HZ form a set of predictions that can be tested in the future using the growing sample of known likely-rocky exoplanets. Specifically, if the carbonate-silicate weathering feedback mechanism operates roughly as expected by climate theory, rocky exoplanets should exhibit a trend of increasing atmospheric CO 2 from the inner edge of the HZ to the outer edge such that temperate surface temperatures are maintained. Lehmer et al. (2020) used a coupled climate and weathering model to investigate the dependence of this trend on practical geophysical and physiochemical differences that are likely to exist between exoplanets. They reaffirmed the existence of such a trend in their models, but found that significant scatter may make it difficult to distinguish observationally, instead suggesting that the 2D distribution of planets in the flux-CO 2 phase space may yield a more reliable test of the HZ hypothesis. Thus, although no single exoplanet can offer a definitive test of the habitable zone, the entire exoplanet ensemble provides a population that, in theory, can.
However, testing for a statistical comparative planetology trend in atmospheric composition will be challenging because the trend itself is not directly observable, but must be properly identified by synthesizing the results of many individual inferences. Since exoplanet atmospheric compositions must be inferred from spectroscopic observations using retrieval models, any trends in atmospheric composition fall into the category of multilevel or hierarchical inference problems. While numerous trends in exoplanet atmospheric composition have been suggested as a means to understand exoplanet habitability (e.g., Turbet et al. 2019;Checlair et al. 2019;Bixel & Apai 2020;Checlair et al. 2021;Bixel & Apai 2021), a consistent framework to tackle the hierarchical atmospheric retrieval problem-going from spectroscopic observations to population-level atmospheric trends-has not been presented.
In this paper, we present a novel retrieval methodology that enables inferences from observations of multiple planets to be combined and synthesized to constrain population-level atmospheric characteristics, and we apply the model to an idealized population of potentially habitable exoplanets to test for the predicted carbonatesilicate weathering CO 2 trend. This is achieved using a new hierarchical Bayesian atmospheric retrieval (HBAR) modeling approach using the importance sampling formalism from Hogg et al. (2010). This use case is a natural extension of a classical hierarchical Bayesian parameter estimation problem (Gelman et al. 2013), which have been highly successful in numerous exoplanet population studies (e.g. Hogg et al. 2010;Rogers 2015;Wolfgang et al. 2016), and recently applied to the atmospheric characterization of hot Jupiters using Spitzer eclipse measurements (Keating & Cowan 2021). However, such methods have yet to be implemented for exoplanet atmospheric retrievals. Typically, hierarchical models are used to properly account for and determine the underlying population-level prior distributions from which an entire population of astrophysical objects are sampled. For instance, the mass-radius relationship has been constrained for sub-Neptune sized planets using mass and radius inferences across an ensemble of known exoplanets (Wolfgang et al. 2016). While completely novel to exoplanet atmospheric retrievals, a similar approach has been applied to Earth remote sensing aerosol retrievals by leveraging a hierarchical model with a built-in spatial dependence to capture spatial smoothness (Wang et al. 2011).
In Section 2, we describe our standard and hierarchical retrieval methods. In Section 3, we present the results of our atmospheric retrieval modeling, and then use them to infer CO 2 trends with our HBAR model and perform a population-level model comparison. In Section 4 and Section 5 we provide a discussion and conclusion of our findings, respectively.

METHODS
We begin by describing our nominal exoplanet atmospheric retrieval model in Section 2.1, which shares many common traits with other retrieval codes in the literature. We then detail the hierarchical modeling approach and how it interfaces with the standard retrieval framework in Section 2.2.

Nominal Retrieval Model
We use the Spectral Mapping Atmospheric Radiative Transfer for Exoplanet Retrieval model (smarter) to solve the Bayesian inverse problem on simulated terrestrial exoplanet transmission spectra (Lustig-Yaeger 2020;Lustig-Yaeger et al. 2021, in prep). We provide a brief description of the smarter model below, which is limited to the essential components for this work, but refer the reader to Lustig-Yaeger (2020) and Lustig-Yaeger et al. (2021, in prep) for a complete description of the smarter retrieval model and its rigorous validation using exoplanet-analog observations of Earth's infrared transmission spectrum.

Forward Model
smarter relies on the Spectral Mapping Atmospheric Radiative Transfer model (smart) as the core of the forward model used to simulate line-by-line transmission spectra for transiting exoplanets (Meadows & Crisp 1996;Crisp 1997;Misra et al. 2014) using the ray tracing formalism described in Robinson (2017). In turn, smart leverages the DIScrete Ordinate Radiative Transfer (DISORT; Stamnes et al. 2017) model to solve the radiative transfer equation. The Line-By-Line ABsorption Coefficient code (lblabc; developed by D. Crisp;Meadows & Crisp 1996) is used to calculate molecular vibrational-rotational absorption coefficients for input into smart radiative transfer calculations. lblabc combines information about the atmospheric state with HITRAN line-parameter and isotope information from the HITRAN2016 line list (Gordon et al. 2017) to calculate gas absorption coefficients as a function of pressure, temperature, and wavenumber. Collisionally-induced absorption (CIA) data are used for CO 2 -CO 2 (Moore 1972;Kasting et al. 1984;Gruszka & Borysow 1997;Baranov et al. 2004;Wordsworth et al. 2010;Lee et al. 2016) and N 2 -N 2 (Lafferty et al. 1996;Schwieterman et al. 2015b).
For simplicity in this study we consider one-dimensional atmospheres composed of N 2 , CO 2 , and H 2 O with isothermal temperature-pressure (TP) profiles and evenly-mixed gas abundances. While plainly limited, this combination of gases is consistent with the climate modeling work that underpins the HZ (e.g. Kasting et al. 1993;Kopparapu et al. 2013). We allow the (log 10 ) volume mixing ratios (VMRs) of CO 2 and H 2 O to freely vary within the forward model, but set the N 2 VMR to the residual VMR such that the VMRs of all three gases sum to unity, as in numerous retrieval studies (e.g. Feng et al. 2018;Krissansen-Totton et al. 2018;Barstow et al. 2020). In addition to the (log 10 ) VMRs of CO 2 and H 2 O, we fit for the isothermal temperature T 0 (in Kelvin), the solid-body surface reference planet radius R 0 (in R ⊕ ), and the surface reference pressure P 0 (in log 10 Pascals). Although we do not formally include clouds in this study, the reference radius and pressure can be used to account for an opaque gray cloud top. In total, we use a five parameter state vector, ω, for our forward model, subject to the following uninformative priors on ω: where U(lower, upper) denotes a uniform distribution with finite probability between the lower and upper bounds. We use the function g(ω) to denote the forward model transformation of parameters ω into wavelength dependent spectroscopic units that can be directly compared to the data.

Inverse Model
We use the dynesty nested sampling code (Speagle 2020;Skilling 2004) to solve the Bayesian inverse problem for the posterior probability distribution function (PDF) of our forward model parameters given transmission spectrum observations. A standard χ 2 loglikelihood function is used to calculate the probability of the transmission spectrum data (∆F = (R p /R s ) 2 ) given the model parameters, where σ j is uncertainty on the spectrum for the jth observed wavelength. Adding the log-likelihood to the logarithm of the aforementioned uninformative priors yields the unnormalized log-posterior that can be sampled with dynesty. We run dynesty with 1000 live points and take model convergence to be achieved when the estimated contribution of the remaining prior volume to the total evidence (Ẑ) falls below ∆ lnẐ < 0.5 between consecutive iterations. This procedure yields K equally weighted samples from the posterior distribution, where K ≈ 16,000 for our experimental setup, as we will see in Section 3.2.

Hierarchical Modeling Approach
We employ the importance sampling hierarchical model described in Hogg et al. (2010) originally presented for inferring the eccentricity distribution of exoplanets. This model has numerous advantages over a traditional fully coupled hierarchical approach. First, all atmospheric retrievals can be pre-computed, independently, and potentially in parallel, using traditional methods and posterior sampling approaches. This assumes that there are no likelihood covariances between parameters from different planets, but benefits computationally from not being required to sample the corresponding high dimensional spaces.
Following closely with the derivations presented in Hogg et al. (2010), for any exoplanetary spectrum n, there are ω n parameters that we attempt to infer as defined in Section 2.1. Consider that we have N exoplanets n (1 ≤ n ≤ N ), each of which has M n transmission spectrum measurements ∆F nj (or wavelength resolution elements). For every exoplanet n, the set of spectroscopic measurements is modeled using a radiative transfer forward model. This is given by where the function g n (ω n ) is the spectroscopic forward model described previously in Section 2.1.1, which is parameterized in terms of the five aforementioned dimensions of ω n , and has an additive noise component drawn from a normal distribution (N ) with variance σ 2 nj (the uncertainty of the jth observed wavelength of the nth star-planet system). The model of all N planets has 5 × N continuous parameters in the larger list {ω n } N n=1 . We note that it is the high dimensional 5 × N parameter space that makes a fully coupled hierarchical model (e.g. using PyMC3) computationally inefficient, as it must simultaneously infer all 5 × N parameters.
The likelihood L n for the five parameters ω n for the nth exoplanet spectrum is the probability of the data D n for planet n given the parameters ω n L n ≡ P(D n |ω n ).
This is the standard Bayesian likelihood function previously defined in Section 2.1.2. Now for each planet n, suppose that we have inferred (using the inverse model described in Section 2.1.2) or been provided with a Kelement sample from a posterior PDF created from the likelihood and an uninformative prior PDF P 0 (ω n ): where Z n is the marginal likelihood or evidence, and is simply a normalization constant. We expect that the prior PDF P 0 (ω n ) is uninformative. For every exoplanet n this posterior sampling is a collection of K equally weighted samples k, each itself a set of five parameters ω nk . The joint likelihood L of all parameters for all exoplanets n in the dataset is simply the product of the individual likelihoods: As stated in Hogg et al. (2010), this makes the assumption that there are no likelihood covariances between the parameters of different exoplanets n. Now we want to reframe the problem slightly and instead consider the likelihood L α for the set of (hyper)parameters α that are used to define an updated prior probability on the CO 2 abundance P α (f CO2 ). It is this updated prior that will be used to encode and characterize population-level trends in atmospheric parameters, which, once known, is a better choice of prior on CO 2 than our original uninformative choice. The probability of the entire data ensemble given the populationlevel parameters α is then given by: This joint likelihood can be expressed as the product of N marginalization integrals over parameters ω n , which can be factored further by assuming that the data D n depend on α only through ω n 1 (i.e., P(D n |ω n , α) = P(D n |ω n )), such that We recognize the first term in the integrand as the individual likelihood for the nth planet spectrum, but the second term-the probability of ω n given α-is critical and allows us to re-weight the integral using our new hierarchical prior, Equation (14) divides out the contribution from CO 2 to the original uninformative prior and multiplies through by the new CO 2 prior that depends on the hyperparameters α. Substituting Equation (14) into Equation (13) and recognizing that the product of the individual likelihood times the original uninformative prior is the posterior via Equation (7), we arrive at the following N multidimensional integrals: Note that we have dropped the dependence of Equation (15) on the marginal likelihood Z n because typical posterior inference applications only evaluate Equation (7) up to the unknown normalization constant. While Equation (15) looks computationally exhausting, the fact that we have already obtained posterior samples simplifies the integral substantially. As articulated in Hogg et al. (2010), since all probability integrals can be approximated as sums over samples, we can employ the posterior sampling approximation to obtain where k runs over all posterior samples K, and the sum simply contains the ratio of the new prior PDF that we want to infer P α (f CO2nk ) to the uninformative prior PDF that was used in the original retrieval inference. Although the individual posteriors P(ω n |D n ) do not explicitly appear in Equation (16), they are implicitly contained within the distribution of K samples. With the likelihood L α as defined in Equation (16), it is straightforward to infer posteriors on the population-level parameters α using Bayes' Theorem, where P(α) is the (hyper)prior PDF for the hyperparameters α.
Within this section we have kept the importance sampling derivations as agnostic as possible to the specifics of the population-level trend(s) under consideration to ensure that the methods can be readily adapted to other problems. Critically, we have not yet specified the form of the CO 2 trend, the hyperparameters that define α, or their respective hyperpriors. We refer the reader to Section 3 (particularly Equation (18), Equation (19), and Equation (20)) for our specific implementation and the subsequent results.
In the original Hogg et al. (2010) formulation of importance sampling, it was assumed that the original exoplanet data was not in-hand, but that the posteriors had been obtained from a colleague or another research group for further analysis. While this is certainly a circumstance that may motivate the use of importance sampling for atmospheric retrievals, we found that it was crucial, perhaps necessary, to perform the hierarchical analysis in a subsequent step following the completion of a uniform set of retrievals, due to the excessive computational expense of simultaneously inferring the population parameters α along with all of the 5 × N individual system atmospheric parameters ω n . This intractability stems from the inherent computational expense of retrieval codes, which solve the radiative transfer equation at each step in the spectral inference. However, retrieval codes with exceptionally fast forward models may prove important for exploring the cost-benefit analysis of HBAR methods that simultaneously infer individual and population parameters. This may be an opportunity for machine learning augmented retrievals (e.g., Nixon & Madhusudhan 2020;Himes et al. 2020;Hayes et al. 2020).

RESULTS
First, in Section 3.1, we present a uniform set of spectral models for an idealized population of terrestrial exoplanets that by design exhibit the silicate weathering CO 2 trend of interest. Second, in Section 3.2, we use our new HBAR model to infer the atmospheric CO 2 trend across the population of synthetic exoplanets. Third, in Section 3.3, we conduct a population-level model comparison to determine the robustness of the inferred trend relative to the null hypothesis and other functional forms for the silicate weathering relation.

Spectral Models
We generate a set of transmission spectra that will allow us to empirically test for the existence of the habitable zone as described in Bean et al. (2017). To limit the number of confounding factors in this study, we assume that the set of N exoplanets with observed spectra are identical to one another except for their stellar irradiation and the quantity of CO 2 and N 2 in their atmospheres, which follows the silicate weathering feedback trend that we impose. We assume the planets possess 1 bar Earth-like atmospheres composed only of N 2 , CO 2 , and H 2 O. We use globally averaged Earth vertical thermal and H 2 O profiles (from Robinson et al. 2010Robinson et al. , 2011Schwieterman et al. 2015a) to satisfy the assumption that each planet is habitable, but neglect changes in these profiles that would be expected from self-consistent climate modeling across the HZ to limit our focus to observables due to CO 2 . Furthermore, we assume that all simulated planets are Earth-sized (1 M ⊕ , 1 R ⊕ ) and orbit TRAPPIST-1 with the same transit duration as TRAPPIST-1e (Gillon et al. 2017;Agol et al. 2021), which simply provides a tangible point of comparison to judge the plausibility of such an analysis in the context of current exoplanet targets and observing capabilities.
To simplify the underlying model for our population trend, we fit an analytic function to the predicted CO 2 volume mixing ratios calculated by Bean et al. (2017) using a 1D radiative-convective climate model. We used the following "Gaussian-like" functional form where S ⊕ is the stellar irradiation incident on the planet relative to Earth and µ, σ 1 , and δ are free parameters which we determine to be 0.04727, 0.5372, and 0.4376 respectively by minimizing the squared residuals. We assume that the remainder of the atmospheric volume is filled with N 2 , and then calculate the mean molecular weight of the atmosphere self-consistently. Figure 1 shows our resulting transmission spectrum models at 1 cm −1 wavenumber resolution (left panel) that correspond to the assumed trend in CO 2 with stellar irradiation (right panel) due to the carbonate-silicate weathering feedback mechanism for N = 20 theoretical exoplanets. By design, the spectra exhibit differences due solely to the volume mixing ratio of CO 2 (and implicitly N 2 ), which are small relative to the total transit depth. Two competing effects shape the observable characteristics of the spectra shown in Figure 1: the CO 2 optical depth and the atmospheric mean molecular weight. As the CO 2 abundance increases the CO 2 optical depth increases, and the weak CO 2 bands, primarily seen between 1 − 2 µm, increase in absorption strength. The opposite is seen for the saturated CO 2 bands at 2.7 and 4.3 µm. The increase in CO 2 causes the saturated bands to decrease in absorption strength as the mean molecular weight of the atmosphere increases from N 2dominated (28 g/mol) to CO 2 -dominated (44 g/mol), and the atmospheric scale height decreases correspondingly. In general, these subtle spectral differences must be sufficiently resolved in each observed spectrum for the population-level model to infer a meaningful trend.
We used the PandExo JWST noise model (Batalha et al. 2017; to simulate synthetic transmission spectrum observations using the Near-Infrared Spectrograph (NIRSpec) Prism instrument (Bagnasco et al. 2007;Ferruit et al. 2014). We used the same PandExo simulation setup as Lustig-Yaeger et al. (2019) assuming the partial saturation strategy for the NIRSpec Prism  and no assumed noise floor. Figure 2 shows the precision of the Prism spectra for TRAPPIST-1e using 10 and 100 stacked transits, compared against the standard deviation among the spectra shown in Figure 1. Based on the fact that the 1σ spectral uncertainties for 100 transits is of similar magnitude to the deviations caused by CO 2 , we conclude that approximately 100 observed transits may be required, for each planet, to obtain spectra with high enough precision to clearly resolve the CO 2 bands in sufficient detail to distinguish between the atmospheres and resolve the trend. While this is an objectively large number even for a single target, we adopt the spectral uncertainties corresponding to 100 stacked transits for each target to ensure that the next stage in the analysis will contain enough information to properly test this population trend with our HBAR model.

An Empirical Test of the Habitable Zone
We performed a uniform set of retrievals on the 20 spectra shown in Figure 1 that evenly span the habitable zone range of stellar irradiation with spectral resolution and noise appropriate for 100 transits with JWST/NIRSpec Prism. Figure 3 shows a representative corner plot of the 1D and 2D marginalized posterior distributions for the inferred physical planetary parameters for the case with the ∼83% atmospheric CO 2 -the fourth largest CO 2 abundance in the sample. The covariance between the isothermal temperature and CO 2 abundance shows a degeneracy due to the dependence of the atmospheric scale height on the temperature and mean molecular weight. Similarly the reference radius and pressure show an expected degeneracy that represents the set of radii and pressures that maintain the transmission spectrum continuum near the observed value. The retrieved H 2 O abundance is consistent with stratospheric values and, notably for the purpose of this investigation, the CO 2 abundance is constrained to within ±0.5 dex. The upper right of Figure 3 shows the median model transmission spectrum obtained from fit- Figure 1. Transmission spectrum models (left) for a sample of 20 hypothetical rocky exoplanets generated following trend in CO2 with stellar irradiation predicated on the assumption of a functioning carbonate-silicate weathering feedback mechanism (right). Planets near the outer edge of the HZ at lower stellar irradiation are predicted to possess atmospheres with higher CO2 abundances than those near the inner edge of the HZ, and such variations in CO2 manifest in observable features in the transmission spectrum.  Figure 1 with varying CO2 (black) compared with the 1σ observational uncertainties for 10 (magenta) and 100 (orange) stacked transits of TRAPPIST-1e with JWST's NIRSpec Prism. The spectroscopic differences caused by the silicate weathering CO2 trend are less than 10 ppm for Earth-sized exoplanets transiting TRAPPIST-1, indicating that approximately 100 stacked transits may be necessary to clearly resolve the deviations in CO2 mixing ratios across HZ instellations.
ting the synthetic JWST data with bounding envelopes to represent the upper and lower 1σ and 3σ credible intervals. The spectral models shown were derived from 500 random samples from the posterior distribution.
Although the retrieval results shown in Figure 3 are only for one representative case from our sample, the other 19 retrievals show similar results with expected differences caused by the different underlying CO 2 abundance and the propagation of random Gaussian scatter in the spectrum through the inference procedure. On average each retrieval with dynesty yielded K ≈ 16,000 equally weighted samples from the respective posteriors. Next, the ensemble of posteriors obtained from the individual planets will be used to infer population-level parameters, in an attempt to retrieve the silicate weathering feedback trend that we injected into the sample.
With the posteriors from our uniform retrieval analysis in hand, we now wish to infer population-level parameters for the CO 2 versus stellar irradiation trend by running MCMC on the HBAR importance sampling model. The HBAR likelihood function is given by Equation (16), where the original prior P 0 (f CO2nk ) on the logCO 2 abundance is uninformative U(−12, 0) and the updated prior P α (f CO2nk ) is calculated from the analytic relationship provided in Equation (18). Specifically, the updated prior is a function of hyperparameters α and is taken to be normally distributed, where σ 2 is the standard deviation of the Gaussian distribution that lends high probability to values of the CO 2 population trend that lie closest to the original CO 2 posterior samples. Thus, for our HBAR importance sampling model, we have the free hyperparameters α ≡ [µ, σ 1 , δ, σ 2 ] that we seek to infer, subject to parameters retrieved from fitting the transmission spectrum using the smarter model. The subplots along the diagonal show the input atmospheric profiles as a function of pressure that were used to generate the synthetic spectrum for comparison against the vertically homogeneous (isothermal and evenly mixed) atmospheric models that were retrieved. The upper right inset shows the median retrieved spectrum bounded by 1σ and 3σ credible intervals derived from posterior samples. This particular retrieval result was for the case with ∼83% CO2 and is representative of similar results obtained for the other planets in our synthetic sample. the following uninformative hyperpriors on α: where N 1/2 (0, 1) refers to a half-normal distribution with a mean of 0 and a standard deviation of 1. At this point, we are ready to evaluate our HBAR model and infer the hyperparameters α. To recap, we now have a specific population-level model for the CO 2 trend (Equation (18)) that is used to define a new populationlevel prior on the CO 2 abundance (Equation (19)). This allows the likelihood L α to be obtained by evaluating Equation (16). Finally, the uninformative hyperprior (Equation (20)) can be multiplied by the likelihood, as in Equation (17), to infer the desired α posterior PDF.
We used MCMC with emcee (Foreman-Mackey et al. 2013) to infer posterior samples of the hyperparameters α using 20 walkers. We ran the chain until it reached a length of approximately 50× the integrated autocorrelation time, as suggested by the emcee code documentation 2 . Figure 4 shows the MCMC results from our HBAR importance sampling model. The lower left set of panels in Figure 4 show the 1D and 2D marginalized posteriors for the population-level hyperparameters in our HBAR model. Relative to their uninformative hyperpriors, the hyperparameters are well constrained by the inference. The median retrieved trend in CO 2 with stellar irradiation is calculated from the posterior samples and is shown in the upper right panel of Figure 4, bounded by the 1σ and 3σ credible intervals. The upper right panel also shows the retrieved 1σ CO 2 constraints for all 20 independent spectrum retrievals plotted as a function of stellar irradiation. It may be conceptually useful to imagine that we have directly fit the purple population model to the black error bars in CO 2 abundance, while in fact we have actually taken the full set of multidimensional posteriors into consideration in our numerical evaluation of Equation (16) and Equation (17). The true underlying CO 2 trend is also shown for reference, where the characteristic decline in CO 2 abundance with stellar irradiation is well resolved by the population-level inference.

Model Comparison
2 https://emcee.readthedocs.io/en/stable/tutorials/autocorr/ One of the benefits of performing an importance sampling HBAR meta-analysis is that multiple different hypothesized population-level atmospheric trends can be investigated and compared without the need to re-run the computationally expensive retrieval models. We now compare three simpler population models to our previously obtained result to demonstrate how such models can be discriminated. These population-level models for the CO 2 versus stellar irradiation relation include a linear trend, a log-linear trend, and a flat non-trend representing the null hypothesis. Figure 5 compares four different best-fitting models to the population-level CO 2 trend. The Bayesian Information Criterion (BIC) is calculated for each model using the following relation where κ is the number of free population-level parameters estimated by the model (i.e., the number of dimensions in α), N is taken to be the number of individual planet spectra in the population, and lnL α is maximum of the log-likelihood function obtained through optimization. When selecting between multiple models, the model with the lowest BIC is taken to be preferred. We subtract the preferred model BIC from every other model to obtain ∆BIC values for comparison. The ∆BIC values in Figure 5 indicate that the "Gaussian-like" model is preferred and there is strong evidence against all of the models with higher BICs. This is the expected result because we generated the synthetic planetary models using the Gaussian-like population trend, but this serves as a useful demonstration of how population-level trends can be compared using an importance sampling HBAR framework.

DISCUSSION
We conducted an idealized simulated search for a trend in CO 2 abundance with stellar irradiation predicted by the carbonate-silicate weathering feedback for exoplanets in the HZ. Our approach used a novel hierarchical Bayesian model-a first of its kind for exoplanet atmospheric retrievals-to infer population trends in atmospheric characteristics that may prove useful well beyond the scope of this work. We elaborate on our scientific and methodological findings in the following subsections.

Practical Challenges to an Empirical Test of Carbonate-Silicate Weathering
We found that the CO 2 trend predicted by carbonatesilicate weathering will be challenging to infer, potentially limiting the observational feasibility of this sta- Figure 4. Corner plot showing 1D and 2D marginalized posterior distributions for the four hyperparameters used to describe the carbonate-silicate weathering feedback population-level trend for rocky HZ exoplanets. The median inferred population trend is displayed in the upper right inset (solid purple line) with shading to denote 1σ and 3σ credible intervals. Using the importance sampling HBAR method we are able to retrieve statistically robust constraints on the parameters characterizing the silicate weathering population trend in CO2 abundance with stellar irradiation. tistical comparative planetology test of the HZ. As presented in Bean et al. (2017), the use of a large statistical sample of exoplanets is enabled by making relatively low precision spectroscopic observations of each planet. However, as we have shown, the changes to the transmission spectrum of an Earth-like planet (which this investigation must detect) caused by changes in CO 2 abundance (which are consistent with the predicted trend) are quite small at < 10 ppm for TRAPPIST-1elike planets. As a result, to infer the silicate weathering trend at high confidence, we found that precise transmission spectra were required that corresponded to 100 stacked transit observations with JWST of each planet in our idealized 20 planet sample, despite the optimistic assumption that each planet was TRAPPIST-1e-like (while only one such planet is known to exist). Acquiring such a sample of high-precision terrestrial exoplanet transmission spectra is not feasible in JWST's nominal mission lifetime, nor would it be made appreciably easier using a successor to JWST, such as the Origins Space Telescope (OST) concept (Origins Space Telescope Study Team 2019), due to the high probability that yet-to-be-discovered transiting HZ rocky exoplanets will be found in systems less amenable for atmospheric characterization than those transiting TRAPPIST-1 (Gillon et al. 2020).
However, our results do not rule out the possibility of using JWST quality transmission spectra to detect an increase in CO 2 abundance with decreasing insolation. Using our novel HBAR framework, we demonstrated that an optimistic simulated survey would be capable of ruling out the null hypothesis for the silicate weathering feedback CO 2 trend at high confidence. Thus, it stands to reason that fewer transits per planet would be able to resolve the CO 2 population trend with less confidence, while remaining statistically robust. While our results do not reveal the spectral precision of such a transition in statistical confidence, hierarchical Bayesian methods are well-suited to resolve population trends that are not as vividly resolved at the individual level, in particular, for noisy datasets where priors dominate the inference.
The difficulty of precisely measuring CO 2 abundances starkly contrasts against the relative ease with which CO 2 detections are predicted for terrestrial exoplanet transmission spectra with JWST. Numerous reports suggest that CO 2 may be an optimal molecule to target to detect the presence of rocky exoplanet atmospheres Lustig-Yaeger et al. 2019;Fauchez et al. 2019;Pidhorodetska et al. 2020). These proposed atmospheric detections hinge upon the strong and saturated CO 2 bands at 4.3 µm and 15 µm, which are largely insensitive to significant changes in CO 2 abundance (Barstow et al. 2016;Wunderlich et al. 2020). This places a limit on the CO 2 abundance precision that can be retrieved from the spectrum. One method for overcoming this limitation and inferring more precise CO 2 abundances is to resolve and detect (or confidently non-detect) the weaker CO 2 bands in the NIR that are not saturated using precise transmission spectra, as we have shown here.
While we have focused exclusively on transmission spectroscopy, other spectroscopic methods for exoplanet atmospheric characterization may prove more successful at detecting population trends in CO 2 abundance. For example, a next-generation direct-imaging mission that can obtain spectra of Earth-like exoplanets around Sunlike stars (as recommended by the National Academies of Sciences, Engineering, and Medicine 2021), such as the Large UV/Optical/IR Surveyor (LUVOIR; LU-VOIR Mission Concept Study Team 2019) or the Habitable Exoplanet Observatory (HabEx; HabEx Study Team 2019) with access to the 1.6 µm and 2 µm CO 2 bands, may offer more leverage for precise CO 2 abundance retrievals. However, this will need to be demonstrated in a future study since these weak CO 2 bands were omitted from the seminal retrieval work of Feng et al. (2018) due to the relative insignificance of CO 2 in the Earth's visible and NIR spectrum.
Nature is likely to produce more complicated atmospheric trends than what we have investigated here. This may hold particularly true for habitable exoplanets where the presence of stable surface liquid water is appreciated to be dependent on a complex interplay of stellar, planetary, and planetary system-wide factors that may indeed produce an elusive population trend that spans many dimensions (Meadows & Barnes 2018). To this end, recent work by Lehmer et al. (2020) demonstrated that the carbonate-silicate weathering feedback trend in CO 2 with incident flux may be log-linear in form with significant scatter due to individual planet considerations, such as land area for weathering and CO 2 outgassing fluxes. The log-linear trend differs from the nonlinear trend from Bean et al. (2017) considered in this work because the Lehmer et al. (2020) model included temperature and CO 2 feedbacks that cause the surface temperature to decline with semimajor axis (Kadoya & Tajika 2014), rather than remain fixed at 289 K throughout the HZ (Bean et al. 2017). The exact nature of the predicted CO 2 -flux trend in the HZ is unlikely to change our results due to the relative consistency of the two similar hypotheses compared to our posterior constraints on CO 2 . Moreover, these predictions are all likely to be incorrect at some level due to their exclusive reliance on geophysical evidence from Earth. This only further motivates the need for methods that allow us to update our understanding of comparative planetology trends using exoplanet data. Future work could leverage the HBAR framework presented here to investigate the 2D population density trend suggested by Lehmer et al. (2020) and potentially incorporate a third dimension for the surface temperature to better capture the climatic feedbacks expected throughout the HZ (Kadoya & Tajika 2014;Lehmer et al. 2020).
Similarly, Seales & Lenardic (2021) used coupled geophysical models to study the temporal onset of habitability and found that variations in tectonic efficiency from one planet to another may produce a predictable distribution in the CO 2 abundance for an ensemble of planets with the same absolute age. Thus it may be possible to expand the population trends studied here into the system age dimension to relate an inferred distribution in CO 2 back to predictions from geophysical models.

The Future of HBAR
We have taken a first step towards developing a hierarchical Bayesian atmospheric retrieval model for tracking population trends in exoplanet atmospheres through the complicated, non-linear, and degenerate problem of fitting exoplanet spectra. Standard atmospheric retrievals are well known to be computationally expensive due to the requirement of a radiative transfer forward model to fit spectroscopic observations. By using the relatively simple importance sampling method for hierarchical Bayesian modeling from Hogg et al. (2010), we effectively avoid the computationally taxing need to perform retrievals on each planet's spectrum simultaneously, as would be required for a standard hierarchical model. Instead, importance sampling allows for population-level inferences in a straightforward meta-analysis of the posterior samples obtained from a uniform set of standard atmospheric retrieval results. This effectively decouples the computationally expensive retrieval modeling from the hierarchical modeling, such that the hierarchical problem can be readily solved once a uniform set of retrieval results (posteriors) are in-hand. Thus, intrigued readers may find that they already have all of the necessary ingredients to characterize population-level trends within their existing retrieval results.
Characterizing population-level trends in exoplanet atmospheres has use-cases well beyond the habitable zone and offers a critical capability for advancing comparative planetology with current, upcoming, and future telescopes. For example, population studies of extrasolar gas giants are already underway. Seminal work by Sing et al. (2016) analyzed the spectra of 10 hot Jupiters observed with HST and found that the planets exhibited a continuum from clear to cloudy atmospheres which may suggest that clouds and hazes, rather than water depletion during formation, are the cause of weaker-than-expected H 2 O absorption features. However, subsequent uniform retrieval analyses by Barstow et al. (2017) and Pinhas et al. (2019) complicate this picture as their retrieved H 2 O abundances suggest subsolar oxygen and/or supersolar C/O ratios with no clear correlations identified. Additionally, Tsiaras et al. (2018) conducted a population study of 30 gaseous exoplanets and found that about half of the sample had detectable atmospheres via H 2 O absorption features. Future work on this hot Jupiter sample may benefit from the HBAR model described here to infer and compare population trends for different proposed formation and evolutionary pathways. Moving towards smaller planets, Changeat et al. (2020) conducted a uniform retrieval analysis to demonstrate that the ESA-Ariel mission (Tinetti et al. 2016) will be sensitive to trends between the atmospheric chemistry and planetary parameters for a population consisting of mostly sub-Neptune and Neptune size planets (Edwards et al. 2019). As more telescopes dedicated to exoplanet atmospheric characterization come online, and as the number of observed exoplanet spectra grows, the use of HBAR modeling may become crucial for comparative planetology.

CONCLUSION
We implemented a first-of-a-kind hierarchical Bayesian atmospheric retrieval model to characterize populationlevel trends in exoplanet atmospheres. We argue that hierarchical Bayesian models are well suited for this task due to the sophisticated inference methods (retrievals) required to transform the observed spectra of exoplanets into meaningful atmospheric characteristics. In particular, the HBAR model that we implemented using importance sampling offers a computationally tractable approach for performing such multi-level inferences because it requires only the posteriors from a uniform set of traditional retrievals, which can all be performed independently. Marginalizing over the full multidimensional posteriors allows the importance sampling HBAR method to propagate complicated parameter covariances through to the population-level hyperparameters; conserving information that may be lost when analyzing retrieved atmospheric trends using only 1D marginalized posteriors or traditional statistical moments. While this study by no mean represents the end-all-be-all of HBAR modeling, we have taken the first few steps towards a computationally tractable HBAR model that will benefit over time from further application and refinement by the exoplanet community.
We tested the importance sampling HBAR framework on an empirical probe of the HZ using simulated transmission spectra of rocky planets with an injected trend in CO 2 abundance with stellar irradiation that is consistent with predictions for a functioning carbonate-silicate weathering negative feedback cycle. We demonstrated that the HBAR method can be used to (1) accurately constrain population-level parameters that characterize the silicate weathering trend and (2) discriminate between multiple different hypothetical population trends. However, we found that such precise spectroscopic measurements would be required to sense the CO 2 trend in terrestrial exoplanet atmospheres that inferring this particular statistical comparative planetology trend may be infeasible using upcoming missions with transmission spectroscopy capabilities. Nonetheless, the use of the HBAR methods presented here may prove to be an important ingredient for future comparative planetology studies as new theories of planetary atmospheric formation, evolution, and habitability are forged in the crucible of exoplanet demographics.