Importance of Sample Selection in Exoplanet-atmosphere Population Studies

To select a subset of the MTS targets, we explore different ranking metrics using planet insolation flux F_insol, planet radius R_p , host star effective temperature T_eff, and the JWST integration time ${t}_{\exp }$, with the goal of addressing the population-level question: is there an observable transition in atmospheric composition that accompanies varying planet radii and flux? We choose a ranking metric that results in a sample with less variation in T_eff (a rough proxy for the high-energy radiation environment of the planet; Linsky 2014) and F_insol (see Figure 1). Our final, purely empirical merit function applied to down-select the MTS sample is:

$$\begin{eqnarray}\begin{array}{rcl}{\mathrm{merit}}_{\mathrm{mock} \mbox{-} \mathrm{atm}} & = & {t}_{\exp }^{-2}\\ & & \times \,{ \mathcal N }({R}_{p};1.7,1)\\ & & \times \,{ \mathcal N }({T}_{\mathrm{eff}};4000,200)\\ & & \times \,{ \mathcal N }(\mathrm{log}({F}_{\mathrm{insol}});1.5,0.3)\end{array}\end{eqnarray} \tag{ 1 }$$

where ${ \mathcal N }$ (variable; μ, σ) represents a normal distribution with mean μ, standard deviation σ, and lower limit of 0, in the same units as the corresponding variable (R_⊕, K, and log(F_⊕), respectively). For each target, the JWST ${t}_{\exp }$ is calculated as the time required to achieve a 30 ppm spectral precision sampled at R = 100 with NIRSpec G395H at 4 μm. We use PandExo (Batalha et al. 2017b) to compute an optimal duty cycle for each observation, where each transit event is assumed to have a total time of twice the transit duration. The strength of the dependency on ${t}_{\exp }$ of −2 was chosen on a trial-and-error basis to ensure a balanced JWST large program of a reasonable size (the two approved JWST large GO programs for transiting exoplanet science are 142 and 75 hr). If the ${t}_{\exp }$ exponent was too steep, the sample would be biased toward bright targets, with short transit durations (which translates to short orbital periods). If the exponent was too shallow, the program size would not be feasible within a single large program. A value of −2 offered a balance between these two end-cases. Though it is beyond the scope of the analysis, exploring how this exponent affects the chosen sample, and the results of the population analysis, is an important question that we leave to future work.

Zoom In Zoom Out Reset image size

Applying this ranking metric to the MTS list as of 2020 April results in the following TOIs (in rank order): 260.01, 776.02, 836.01, 562.01, 134.01, 175.01, 836.02, 687.01, 776.01, 455.01, 175.02, 186.01, 174.01, 174.02, 402.01, 402.02. We removed TOI 687.01 due to the uncertain period of the planet, and TOI 186.01 because it would saturate the NIRSpec detector. We remove TOI 174.01 and .02 because they have not been validated as a confirmed planets. This results in a sample of 12 targets as shown in Figure 1 (left panel). We fix the number of targets to twelve such that it could feasibly fit into a single large JWST program. Though quantifying how well population parameters can be retrieved as a function of sample size would be an important exploration, it is beyond the scope of this analysis. The total exposure time for this sample (${\sum }_{n=1}^{12}{t}_{\exp }$) is 94.2 hr.

Our second sample represents a binned approach to sample selection, wherein planets of different sizes are equally represented on the target list and are otherwise prioritized by how easy they are to observe. For each planet we calculate a JWST integration time ${t}_{\exp }$ using the same formulation as the QS sample. Then we separate these planets into three R_p bins: 1.0–1.67, 1.67–2.34, and 2.34–3.0 R_⊕. Within each bin we rank the targets by their ${t}_{\exp }$ (discounting targets that would saturate the detector by exceeding 100% full well within two frame times, known as a "hard saturation"), and fill in a list of 12 targets by selecting the top four targets in each bin. The radius bins are roughly chosen to represent planets that are rocky, planets that are likely sub-Neptunes, and those that are in between.

We also tried implementing an insolation flux dimension (4 flux × 3 radius bins), such that there was one target per bin. We first attempted a log-normal spacing between 4 and 4000 F_⊙. However, the "best" target in the highest insolation flux & radius bin (K2-66 b) has ${t}_{\exp }\,=$ 101 hr resulting in an unfeasible total sample time of 201 hr. We also tried extending the bounds of the insolation flux parameter space to cover both lower and higher fluxes. Each of these attempts resulted in samples that had at least two empty radius-flux bins (out of 12 total). In those cases, additional questions arose regarding how to reassign targets. For example, if one radius-flux bin is empty is it more optimal to assign a target to an adjoining bin? Or, is it better to choose the next highest ranked target (e.g., by t_exp) regardless of what bin it may fall in? Though these are interesting questions, addressing them was ultimately beyond the scope of the analysis in this paper.

Therefore, we ultimately opted for a binning method with a single dimension (radius). We emphasize that this method has the potential to create a highly skewed sample (e.g., if all the targets in a single radius bin fell into a narrow flux range). Our 3Bin-sample selection method naturally resulted in a sample that was fairly evenly spaced in flux. Overall, we note that this "missing-bin" problem is common when the number of targets with mass constraints (colored points in Figure 1, middle and right) is small, and therefore we view our choice as a viable test for this analysis.

We choose the same number of planets as the QS sample to ensure that our population relation retrieval results is not driven by differences in sample size. The resulting 3Bin targets, shown in Figure 1 (middle panel) include GJ 357 b (TOI-562.01), LTT 1445 A b (TOI-455.01), L 98-59 d (TOI-175.02), Kepler-21 b, HD 86226c (TOI-652.01), HD 213885 b (TOI-141.01), GJ 9827 d, TOI-776 b, HD 15337c (TOI-402.02), HIP 116454 b, HD 106315 b, and TOI-824 b. We note that five of the targets in this sample overlap with the QS sample. The total exposure time for the 3Bin sample is 106.9 hr.

Our third sample represents an approach to sample selection wherein targets are hand-picked to evenly cover the R_pl–F_insol plane. Hand-picking has the potential to prioritize planets that are novel or unusual in terms of their size or insolation flux. Since we have no direct ranking on observing feasibility, we filter the sample to include only targets with ${t}_{\exp }\leqslant 15$ hr (discounting targets that would hard saturate within two groups), and plot the planet R_p and F_insol values. A 15 hour ${t}_{\exp }$ limit ensures no one planet dominates the total time allocated to the program. Then we choose a sample of 12 targets (again, the same size as our QS sample) across this parameter space. The ESBE targets, shown in Figure 1 (right panel), are TOI-421 b, HD 97658 b, LTT 3780c (TOI-732.02), HD 213885 b (TOI-141.01), HD 86226c (TOI-652.01), GJ 9827 d, TOI-776 b, HD 15337 b (TOI-402.01), GJ 9827c, GJ 357 b (TOI-562.01), LTT 1455 A b (TOI-455.01), and L 98-59 d (TOI-175.02). We note that five (seven) of the targets in this sample overlap with the QS (3Bin) samples, respectively, and that four planets are common across all three samples (see Figure 1). The total exposure time for the ESBE sample is 89.7 hr.

To begin our simulation study, we first define a simple C/O ratio based only on the abundances of three molecules:

$$\begin{eqnarray}&&{\rm{C}}/{\rm{O}}=\displaystyle \frac{{\chi }_{{\mathrm{CO}}_{2}}+{\chi }_{{\mathrm{CH}}_{4}}}{2\ast {\chi }_{{\mathrm{CO}}_{2}}+{\chi }_{{{\rm{H}}}_{2}{\rm{O}}}}\end{eqnarray} \tag{ 2 }$$

where χ represents the abundance of each molecule. We focus specifically on these three molecules as they are the expected dominant sources of opacity in NIRSpec G395H's 3–5 μm region for the planets explored here (see Figure 1). Of course, there are other C and O bearing species that have the potential to affect the C/O ratio, and that would give valuable context clues to the nature of these planet atmospheres. For example, the vertically distributed abundance of CO (along with CH₄) is sensitive to many parameters such as vertical mixing, gravity, temperature, and metallicity (Zahnle & Marley 2014). Additionally, CO along with HCN, C₂H₂, and C₂H₆ have been identified as molecules that could help distinguish the existence of a shallow surface (<10 bar) typical of rocky planets, versus a deep surface (>100 bar) typical of gaseous planets (Yu et al. 2021). These molecules also have absorption bands in the 3–5 μm region. However, as we motivate in Section 3.2, including these physical processes would require complex atmospheric modeling that is beyond the scope of this analysis. Instead, we bypass complex modeling and directly choose how to vary C/O with respect to R_p and F_insol. Our results are therefore based on the specific choices for this injected population relation, which we motivate below.

Next, we must choose how to vary C/O with respect to R_p and F_insol. From in-depth retrieval and information content studies that utilize simulated JWST data (Greene et al. 2016; Batalha & Line 2017), the expected precision obtained on log(C/O) will likely be on the order of 0.5–1.5 dex, if not upper/lower limits. Simply put, for small planets we may only obtain order of magnitude constraints on C/O ratio. Therefore, only the most basic of population relations can be determined. Given our primary goal of determining the feasibility of unearthing a possible trend with JWST-quality spectra and how that depends on the selected sample, we need a relation that can be informative even in the presence of large error bars and that allows for significant variability among the individual planets' C/O ratios. One such population model is a logistic function that classifies a planet as having a carbon- or oxygen-dominated atmosphere ($\mathrm{log}\,{\rm{C}}/{\rm{O}}\geqslant 0$ and $\mathrm{log}\,{\rm{C}}/{\rm{O}}\lt 0$, respectively), where the atmospheric state depends on the R_p and F_insol as

$$\begin{eqnarray}&&\mathrm{ln}\displaystyle \frac{P(\mathrm{log}\,{\rm{C}}/{\rm{O}}\lt 0)}{P(\mathrm{log}\,{\rm{C}}/{\rm{O}}\geqslant 0)}={m}_{r}{R}_{p}+{m}_{f}{F}_{\mathrm{insol}}+b.\end{eqnarray} \tag{ 3 }$$

In this logistic relation, $P(\mathrm{log}\,{\rm{C}}/{\rm{O}}\lt 0)\equiv {P}_{{\rm{O}}}$ is the probability that the atmosphere is oxygen-dominated, $P(\mathrm{log}\,{\rm{C}}/{\rm{O}}\geqslant 0)\equiv {P}_{{\rm{C}}}$ is the probability that the atmosphere is carbon-dominated, and m_r , m_f , and b are the population parameters that define how the relative likelihood of being carbon- or oxygen-dominated varies with respect to the planet's radius and insolation flux. Note that the left-hand side of the relation is not C/O itself, but the probability that a planet has a C/O ratio below 1 relative to the probability that a planet has a C/O ratio above 1. By focusing on the probability that the target variable (C/O in our case) has a value in a certain range instead of focusing on the value of the variable itself, logistic relations define boundaries in the parameter space of the regressor variables (R_pl and F_insol) that simultaneously enable classifications (carbon-dominated if P_O,i < 0.5 or oxygen-dominated if P_O,i ≥ 0.5) and allow the target variable to take on a wide range of values across the population. This parameterization therefore enables us to assess the presence of transitions in planetary atmospheres as a function of radius and insolation flux, even with uncertain atmospheric C/O measurements.

To use Equation (3) in our simulation study, we must first specify values for the population parameters m_r , m_f , and b. To do this, note that Equation (3) describes a plane in the (R_pl, F_insol, ln (P_O/P_C)) space; to define a specific plane, we choose three points that are reasonable based on our current understanding of exoplanet atmospheres and solve the resulting system of three equations for m_r , m_f , and b. Specifically, these points are: ⁸ P_O = 0.5 at (R_pl = 1R_⊕, F_insol = 1000F_⊕), P_O = 0.01 at (R_pl = 3R_⊕, F_insol = 1000F_⊕), and P_O = 0.99 at (R_pl = 1R_⊕, F_insol = 1F_⊕), which gives m_r = −2.30, m_f = −1.53, and b = 6.89. For potentially rocky planets, this P_O represents the hypothesis that those planet's atmospheres are more likely to be oxidized, similar to present day atmospheres of the terrestrial planets in our own solar system. For example, CO₂, not CH₄, is the major carbon-bearing species in Venus, Earth, and Mars (Wayne 1991). For potentially gaseous planets, this P_O corresponds to the hypothesis that these planet's atmospheres are more likely to have near solar C/O ratios. However, we emphasize that these probabilities do not correspond to verified hypotheses in the study of exoplanet atmospheres. Additionally, we note that these probabilities test the optimistic case that P_O actually has a transition with radius from low (near zero) to high (near 1) probability. We acknowledge that many other scenarios are possible, including the scenario where there is no such transition or a weak transition. However, choosing such cases would prevent us from addressing our goal of determining if population-level inferences can be made by studying a sample of planets with JWST. We view our choice in P_O as a starting point.

From this fully specified population relation we then plug in each planet's R_pl and F_insol to obtain P_O,i and P_C,i = 1 − P_O,i . Next, we draw the planet's $\mathrm{log}\,{\rm{C}}/{\rm{O}}$ from a uniform distribution with bounds of (−2, 0) or (0, 2), where the <0 or the >0 range is chosen in proportion to the drawn value of P_O,i . With a planet's specific C/O value in hand, next we must use it to determine the individual molecular ratios CO₂/H₂O, CH₄/H₂O, and CO₂/CH₄. To do this uniquely, we must independently specify at least one of these three molecular ratios; we choose CO₂/CH₄, which we draw from a uniform distribution. We specifically do not choose to draw atmospheric ratios using the assumption of chemical equilibrium since this physical condition will not apply to the full sample of planets here (namely those that are potentially rocky, and/or cool; T_eff < 1000 K). To determine the bounds of this distribution, we first consider that the ability to detect a species in transmission can be roughly determined by assessing the ratio of the cross-sections, σ, weighted by the molecular abundance, ξ. At a given wavelength, λ, if ${\sigma }_{\lambda ,{\mathrm{CO}}_{2}}{\xi }_{{\mathrm{CO}}_{2}}\gg {\sigma }_{{\lambda }_{,}{\mathrm{CH}}_{4}}{\xi }_{{\mathrm{CH}}_{4}}$, CO₂ will dominate the spectrum. Across 3–5 μm, the median ratio of the cross-sections is σ_λ CO₂/σ_λ CH₄ ∼ 10 (similarly, σ_λ CH₄/σ_λ H₂O ∼ 15). Therefore, we choose −2 < log(CO₂/CH₄) < 0, which creates a diverse set of resultant spectra spanning cases with detections of both carbon-bearing species to detections of only a single carbon-bearing species.

With this population relation defined, and C/O and CO₂/CH₄ values drawn for each planet, we have begun to specify the atmospheric state of each planet, which we continue in Section 3.2.

We choose a simple methodology for simulating observed planetary spectra, rather than creating "self-consistent" models (relying on converging temperature, chemical models, and cloud profiles based on initial boundary conditions). This is in line with many theoretical studies that have sought to determine the detectability of super-Earth and sub-Neptune atmospheres in the JWST era (Batalha & Line 2017; Morley et al. 2017; Batalha et al. 2018; Chouqar et al. 2020; Guzmán-Mesa et al. 2020). For example, Morley et al. (2017) created a grid of models with Earth-, Titan-, and Venus-like atmospheres, and Batalha et al. (2018) chose a fixed background gas scenario (e.g., H₂O-rich versus H₂-rich) and varied the remaining trace species in fixed increments. These modeling choices are especially necessary for planets that straddle the super-Earth/sub-Neptune regime, because mass/radius cannot serve as a reliable proxy for H/He envelope mass below about 2.2 R_⊕, at which point internal structure models with significant fractions of heavier gases like H₂O are also able to fit the observed exoplanet masses and radii (Valencia et al. 2006; Rogers 2015).

Here we choose a modeling framework that affords us the opportunity to address the main goal of this work—the importance of sample selection. We use a double-gray analytical parameterization for the temperature–pressure profile (Guillot 2010). We only consider four-molecule atmospheres (H₂, H₂O, CO₂, and CH₄). Similar simplicity has been given to other studies investigating hypothetical atmospheres (e.g., 100% H₂O in Greene et al. 2016; H₂–H₂O in Batalha et al. 2017a). We assume each atmospheric composition to be well-mixed (i.e., uniform with altitude). These well-mixed values are chosen based on the injected population relation.

We note that opting for this simplicity ignores two important factors in determining population-level trends that would exist in nature. First, our simplicity may exclude physical processes that could potentially provide vital context clues regarding the nature of super-Earth/sub-Neptunes. For example, Yu et al. (2021) identified seven chemical species that could help distinguish the existence of shallow versus deep surfaces. Second, our simplicity may exclude physical processes that would make it more difficult to establish trends in atmospheric parameters. For example, the stellar UV flux from each parent star determines to what degree photochemical processes will drive chemical abundances of key molecules such as CH₄ (e.g., Hu 2021). Both processes analyzed in Yu et al. (2021) and Hu (2021) require robust photochemical modeling with complete chemical networks. Therefore, despite these likely effects, incorporating them in a uniform manner is not trivial and beyond the scope of this analysis.

Given the chosen model, and with the chosen C/O, CO₂/CH₄, and H₂O/CH₄ uniquely assigned (see Section 3.1), the final parameters to choose are bulk H₂ fraction and cloud parameters. Higher bulk H₂ fraction increases the scale height of the atmosphere because of a decreased mean molecular weight, and thus increases the magnitude of spectral features. Increased cloud optical depth/decreased cloud pressure has the effect of muting spectral features. For JWST-quality data, the effect of these two parameters on the spectra are degenerate in the infrared (⪆2 μm; Benneke & Seager 2012; Batalha et al. 2017a). Additionally, even for hot Jupiters it is not yet clear how to predict the degree of cloud coverage as a function of planet parameters (e.g., equilibrium temperature and gravity; Wakeford et al. 2019; Alam et al. 2020; Gao et al. 2020). The same is true for super-Earths/sup-Neptunes, which are undoubtedly more difficult targets for high SNR spectra (Crossfield & Kreidberg 2017; Dymont et al. 2022).

Therefore, we start by exploring the cloud-free case with a bulk H₂ fraction fixed to 99%. This H₂ value is chosen to correspond to a Neptune-like H₂/He-fraction (100 × Solar metallicity, Karkoschka 1998). For completeness, we also compute spectra with a bulk H₂ fraction of 90% as well as a "cloudy" case. In our "cloudy" scenario we insert a gray opacity source at 0.1 bars—a case where the observation is limited by the tropopause of the planet, which is defined at 0.1 bars in all solar system planets (Robinson & Catling 2014). We find that if we prescribe a lower H₂ bulk fraction and included the effect of clouds as gray opacity source at a fixed pressure, there are 4–5 planets in each sample with no detectable features. Therefore, all three samples will be similarly encumbered by clouds and increased mean molecular weight via an effective decrease in sample size. We discuss this limitation in our concluding remarks.

With the given prescription for atmospheric abundances and temperature–pressure profiles, we use the PICASO radiative transfer tool (Batalha et al. 2019; Batalha & Rooney 2020) to compute the transmission spectra. Of importance for this analysis are the opacities of H₂O, CH₄ and CO₂, for which we use Polyansky et al. (2018), Yurchenko & Tennyson (2014), and Huang et al. (2014), respectively. Our R = 10⁶ line-by-line opacities are resampled at R = 10⁴ to be suitable for retrievals at R = 100 and are available for download at Batalha et al. (2020). For each transmission spectrum, we use PandExo to compute a simulated JWST observation with NIRSpec G395H, which we then bin to R = 100 for the ultimate retrieval. Lastly, we pair PICASO with the open source Nested Sampler dynesty (Speagle 2020), which implements the algorithm developed by Skilling (2004). Mukherjee et al. (2021) outlines the specific hyperparameters used for the sampler, dynesty.

For each planet we compute the posterior probability distributions for five free parameters: (1) the irradiation temperature of the Guillot (2010)–P(T) profile (in K; prior: U(300, 1200)); (2) H₂ bulk abundance (in dex; prior: U(−6, 0)); (3) H₂O/CH₄ abundance ratio (in dex; prior: U(−6, 6)); (4) CO₂/CH₄ abundance ratio (in dex; prior: U(−6, 6)); and (5) xR_p , a scaling factor to the reported radius derived from the Kepler/TESS transit observation that we arbitrarily define at 10 bars (unitless; prior: U(0.5, 1.5)). For this analysis we focus specifically on the retrieval results of the abundance ratios and combine the posteriors for H₂O/CH₄ and CO₂/CH₄ following Equation (2) to compute a posterior for C/O. In total we run 23 unique planetary atmospheres (some of the selected targets are the same across the three samples) for each of the 10 random trials, for a total of 230 retrievals.

We consider three metrics to evaluate how each hypothetical NIRSpec G395H survey did across all ten trials, and whether there is a clear best practice for constructing atmospheric surveys.

The first metric, M_atm, describes the overall precision and accuracy of the individual-planet atmospheric parameter constraints, which in this study is the posterior probability distribution of C/O. This quantity is computed by combining the retrieved posterior distributions of H₂O/CH₄ and CO₂/CH₄ via the C/O approximation defined in Equation (2). For each trial we compute the chi-squared of the injected (C/O_inj,i ) versus retrieved (C/O_ret,i ) log C/O ratio over all 12 planets, while accounting for the 1σ-posterior width of each planet (σ_i ):

$$\begin{eqnarray}&&{M}_{\mathrm{atm}}=\sum _{i=1}^{12}\displaystyle \frac{{\left(C/{O}_{\mathrm{ret},i}-C/{O}_{\mathrm{inj},i}\right)}^{2}}{{\sigma }_{i}^{2}}.\end{eqnarray} \tag{ 5 }$$

Therefore, the lower the M_atm, the better the hypothetical survey did overall in constraining individual-planet atmospheric composition parameters.

The second metric, M_pop,p, describes the overall ability to constrain the population parameters precisely. For this metric, we compute the volume enclosed in the joint, three-dimensional posterior of the three population parameters (m_f , m_r , and b) at 3σ. We also explored using the 1σ or 2σ posterior volume as the metric; the overall conclusions were unchanged. Similar to M_atm, the lower the M_pop,p, the tighter the width of the joint posterior and the better the survey did at precisely constraining the population parameters.

The third metric, M_pop,a describes the overall ability to constrain the population parameters accurately. For this metric, we compute the effective "distance" between the posterior mode and the injected value. We define "distance" as the difference between the injected value and the posterior mode in the three-dimensional hyperparameter space. Similar to previous metrics, the lower the M_pop,a, the higher the accuracy and the better the survey did at accurately constraining the population parameters.

Figure 6 shows the results of the three metrics for all 10 trials and both prior assumptions on population-level parameters. With regard to the atmospheric parameters, QS and 3Bin each had 4 out of 10 trials with the lowest value of M_atm. In the other two trials, ESBE had the lowest value of M_atm. According to the mean M_atm across all ten samples, the QS- and 3Bin-derived samples performed nearly equally (average lines are overlapping in Figure 6 with M_atm = 7.1).

Zoom In Zoom Out Reset image size

With regard to the population-level metrics, considering uniform priors on m_r and m_f , QS had the lowest value of M_pop,p in 6 out of 10 trials, and the lowest value of M_pop,a in 4 out of 10 trials. The ESBE survey had the lowest value of M_pop,p in 3 out of 10 trials, and the lowest value of M_pop,a in 5 out of 10 trials. This suggests that QS and ESBE performed more similarly while 3Bin consistently had the poorest performance. We note that the variability across trials is quite high and that an estimate of the metrics' variances given by the 20th–80th quantile range (approximately the quantiles that would correspond a 1σ uncertainty interval for a normally distributed metric, which these are not) would yield overlapping bands.

Furthermore, we test the robustness of these results against our assumptions of priors on the population-level parameters, m_r and m_f , as described in Section 4. Considering the arctan priors on m_r and m_f , the QS sample achieves the lowest M_pop,p value in 4 out of 10 trials while ESBE achieves the lowest value in 6 out of 10 (for the uniform prior it was 6/10, and 4/10, respectively). For the M_pop,a metric, QS achieves the lowest value 6 out of 10, ESBE achieves the lowest value 3 out of 10, and 3Bin achieves the lowest value once (for the uniform prior it was 4/10, 5/10, and 1/10, respectively). Though the overall conclusions drawn from the metrics are not changed, the choice of prior does affect the results on a trial-by-trial basis.

In Table 1 we show a break down of which sample achieved the lowest metric values for each trial and prior choice. In 3 out of 10 trials, the sample with the lowest M_pop,p metric was changed, and in 2 out of 10 trials the sample with the lowest M_pop,a metric was changed (shown as bolded text in Table 1). In all cases, the arctan prior changed the "winning" sample to either ESBE or QS, not 3Bin. In other words, our overall result that the 3Bin-derived sample achieves the poorest performance is not prior dependent. We verify this by showing in Figure 7 the representative posteriors of Trial 9, where the choice of prior did affect the population-level metric but, qualitatively, the posterior distributions of the population parameters are not drastically affected.

Zoom In Zoom Out Reset image size

Ultimately, 3Bin consistently had the poorest population parameter constraints, while QS and ESBE were relatively tied. According to the mean across 10 trials, with uniform priors on m_r and m_f , QS obtained slightly higher precision constraints, while ESBE obtained slightly higher accuracy constraints, and with arctan priors, QS obtained higher precision and accuracy constraints. However, we reiterate that a mean for the population parameters does not capture the variability across trials, which is significant as is evident from Figure 6. The variability exhibited from trial to trial with these twelve-planet samples illuminates the critical need for samples larger than what can be feasibly done in a single JWST cycle.

Ultimately, we consider QS and ESBE to have done equally well according to the number of times these samples achieved the lowest value of M_pop,p or M_pop,a. This particular result leads to an important conclusion: Even when the 3Bin sample had the "best" precision and accuracy on log(C/O), compared to the QS- and ESBE-samples, 3Bin resulted in population parameters that were worse in terms of both accuracy and precision. This finding showcases why sample selection must be carefully considered for population-level analyses: simply optimizing for the best constraints on individual-planet atmospheric parameters will not necessarily lead to accurately and/or precisely constrained population parameters. Overall, there was no one-to-one mapping between the sample that had best-constrained atmospheric parameters across all 12 planets and the sample that ultimately had the best-constrained population parameters.

It is important to reiterate that five QS sample targets overlap with the ESBE sample, seven 3Bin-sample targets overlap with the ESBE sample, and five QS sample targets overlap with the 3Bin sample. Because only one draw was done per planet per trial, differences between the success of the hypothetical survey (e.g., ESBE versus 3Bin) are driven by the planets that differ between the surveys. In the case of "ESBE versus 3Bin," this would be just the difference between 5 out of 12 planets. These overlaps in the samples strongly motivate increasing the parent sample from which planets for this type of survey could be selected (i.e., those with mass constraints).