Robustness of C/O and mean molecular weight under disequilibrium chemistry

Introduction

Transmission spectroscopy is a primary route to atmospheric characterisation for transiting exoplanets, especially those with extended atmospheres and frequent transits. Wavelength-dependent transit depths arise because starlight passes through the planetary limb and acquires opacity signatures from gases and aerosols [54, 14, 32, 8]. The James Webb Space Telescope (JWST) now provides broad wavelength coverage and high precision, enabling detections of multiple molecules in hot Jupiters, warm gas giants, and sub-Neptunes [53, 6, 64, 39, 7]. Improved data quality does not remove a central limitation: interpretation often compresses diverse spectra into a small set of compositional metrics whose reliability depends on modelling assumptions.

Retrieved abundances of common molecular opacity sources such as $\mathrm{H}_2\mathrm{O}$, CO, $\mathrm{CO}_2$, and $\mathrm{CH}_4$ often carry large fractional uncertainties. Across representative retrievals of hot Jupiters, warm Neptunes, and sub-Neptunes, reliably detected species typically have fractional abundance uncertainties of approximately 50–150%, while marginal detections and upper limits are substantially less constrained [37, 2, 21]. Therefore, atmospheric retrievals commonly report derived quantities such as the carbon-to-oxygen ratio ($C/O$), elemental carbon abundance ($C/H$), and mean molecular weight $\mu$, all of which are calculated from uncertain molecular mixing ratios. Their uncertainties depend on wavelength coverage, spectral resolution, data quality, cloud treatment, and whether individual molecules are strongly required by the data or only moderately favoured in Bayesian model comparison [15, 21]. Typical thresholds for strong and moderate evidence are $\ln\mathcal{B} \gtrsim5$ and $\ln\mathcal{B} \sim2.5$, respectively [40, 3, 63, 37]. Molecular abundance uncertainties do not propagate equally into $C/O$, $C/H$, and $\mu$. Elemental ratios depend on sums and ratios of trace species, whereas $\mu$ weights each species by its molecular mass and may be comparatively insensitive to uncertainty in a single trace molecule when the atmosphere is dominated by $\mathrm{H}_2$. Therefore, the relevant question is not only what value is inferred, but also how strongly each derived metric amplifies or partially cancels the uncertainties in the retrieved molecular abundances.

Atmospheric processing may also affect these metrics differently. Assumptions about clouds, thermal structure, chemistry, reference pressure, and radius normalisation influence the inferred molecular abundances even when the observations are precise. Vertical mixing and photochemistry further modify the abundances at the low pressures probed by transmission spectroscopy. Therefore, a robust compositional metric should satisfy two conditions: it should carry a small propagated uncertainty from the retrieved molecular abundances, and it should remain stable under plausible changes in atmospheric chemistry and transport.

C/O is widely used because it controls the distribution of carbon and oxygen among the dominant molecular carriers in hydrogen-dominated atmospheres and can strongly re-shape the spectrum [38, 44, 23]. C/O is also interpreted as a tracer of planet formation through disk snowline chemistry [19, 50] (Öberg et al. 2011). C/H provides a complementary absolute abundance scale. We define C/H as the elemental carbon abundance relative to hydrogen, thereby avoiding ambiguity with bulk metallicity.

Mean molecular weight, $\mu$, is motivated by the geometry of radiative transfer. In a hydrostatic atmosphere, the pressure scale height satisfies $H \propto T\mu^{-1}g^{-1}$ where $T$ is temperature and $g$ is gravitational acceleration. Therefore, mean molecular weight helps determine the amplitude of spectral features in transmission [36, 20]. This relation motivates the interpretation of muted spectra as evidence for atmospheres with high $\mu$ [8]. Transmission spectra, however, lack an absolute normalisation and contain degeneracies among composition, reference pressure, and transit radius [27, 11, 63]. Aerosols introduce another degeneracy because clouds and hazes can flatten spectral features and mimic the effect of high $\mu$ [22, 5, 42, 9]. Precise observations are required to distinguish clouds at high altitude from atmospheres with high $\mu$ [35].

Sub-Neptunes strengthen the motivation for considering $\mu$ in population studies. Heng et al. (2025) predicted a gradient in $\mu$ across the observed radius valley [24, 49], driven primarily by interior oxygen fugacity rather than atmospheric carbon enrichment. This prediction suggests that $\mu$ may encode interior conditions. It also exposes a modelling challenge: connecting bulk or outgassed composition to an observable atmosphere requires tracking how chemistry, mixing, and irradiation process molecular species at the pressures probed by transmission spectroscopy.

Photochemistry provides a direct test of metric robustness because it alters atmospheric composition at low pressures. Transport and photochemistry modify vertical abundance profiles and shift the observable composition away from local thermochemical equilibrium [51, 17, 61] (Moses et al. 2011; Hu & Seager 2014). The detection of the photochemical product SO$_2$ in WASP-39b demonstrated that disequilibrium chemistry can directly affect observed spectra [53, 55]. Therefore, C/O, C/H, and $\mu$ may respond differently to the same physical perturbation.

Here we assessed the reliability of $\mu$ and C/O as compositional metrics in transmission spectroscopy for sub-Neptunes and hot Jupiters. First, we analytically propagated uncertainties in retrieved volume mixing ratios through a chemical-equilibrium network containing carbon, hydrogen, and oxygen to compare the intrinsic uncertainties in $\mu$, C/O, and C/H. Second, we used VULCAN forward models to quantify how these metrics responded to photochemistry, vertical mixing, metallicity, irradiation, and planetary type (Tsai et al. 2017, 2021). Finally, we applied the same uncertainty propagation framework to molecular abundances from published retrievals. Our goal was not to replace elemental ratios with $\mu$, but to determine when $\mu$ was robust and when it depended strongly on the atmospheric model. By benchmarking $\mu$ against C/O and C/H across hot Jupiter and sub-Neptune regimes, we provided a framework for interpreting JWST transmission spectra in population studies.

Methods

Intrinsic uncertainty propagation in the C–H–O system

We used a hydrogen-dominated C–H–O atmosphere to estimate the intrinsic uncertainty associated with each compositional metric. The calculation used an equilibrium chemistry network containing CH$_4$, H$_2$O, CO, CO$_2$, and C$_2$H$_2$. These species captured the dominant carbon- and oxygen-bearing molecules in hot hydrogen-rich atmospheres and matched the C–H–O system used in previous analytic equilibrium-chemistry studies [29, 28].

We specified the temperature $T$, pressure $P$, elemental ratio C/O, and elemental abundances $n_{\rm O}$ and $n_{\rm C}$. We fixed $n_{\rm O}=5\times10^{-4}$ and sampled $T$, $P$, and C/O over 100 realisations using $T\sim\mathcal{U}(800,2500)\ {\rm K}$, $\log(P/{\rm bar})\sim\mathcal{U}(-3,2)$, and C/O $\sim\mathcal{U}(0.1,1.0)$. We then calculated the carbon abundance $n_{\rm C}$ from the selected C/O ratio according to

The temperature and C/O ranges span the hot-Jupiter and sub-Neptune regimes considered in the forward models. The pressure range was broader than the nominal transmission region used later because this calculation sampled the equilibrium compositions available across a vertically mixed atmosphere, rather than the composition at a single pressure level. In this sense, the sampled pressures approximated the thermochemical reservoir that may be transported into the observable atmosphere by vertical mixing.

For each sampled state, we computed equilibrium constants and molecular mixing ratios $X_i$ by solving the polynomial system for the C–H–O network. The temperature-dependent Gibbs free-energy differences, $\Delta G_i(T)$, were taken from the National Institute of Standards and Technology–Joint Army-Navy-Air Force (NIST–JANAF) Thermochemical Tables and linearly interpolated in temperature, following Heng et al. [29] and Heng & Lyons [28]. The NIST-JANAF tables provide recommended temperature-dependent thermochemical quantities, including Gibbs free energies and equilibrium constants, for gas-phase species (National Institute of Standards and Technology 1998).

For each realisation, equilibrium mixing ratios were computed for H$_2$O, CO, CO$_2$, and CH$_4$, with H$_2$ obtained from number conservation. Elemental budgets were defined as $N_{\rm O}=X_{\rm CO}+2X_{\rm CO_2}+X_{\rm H_2O}$ and $N_{\rm C}=X_{\rm CO}+X_{\rm CO_2}+X_{\rm CH_4}$, so that O/H $=N_{\rm O}/(2X_{\rm H_2})$, C/O $=N_{\rm C}/N_{\rm O}$, and $\mu=2X_{\rm H_2}+18X_{\rm H_2O}+28X_{\rm CO}+44X_{\rm CO_2}+16X_{\rm CH_4}$.

To isolate the intrinsic error amplification of each compositional metric, we assigned every molecular abundance the same fractional uncertainty, $\delta$:

We treated the abundance errors as independent and uncorrelated between species. This assumption isolated how the algebraic form of each metric propagates molecular abundance errors, and was not intended to reproduce the covariance structure of a full atmospheric retrieval. Under this assumption, the fractional uncertainty in a weighted sum $A=\sum_i a_iX_i$ is

and for a ratio $z=x/y$ of independent quantities,

We neglected uncertainty in $X_{\mathrm{H}_2}$, which is not directly measured. Applying these expressions gives the propagated fractional uncertainties

and

The first expression uses only the uncertainty in the O/H numerator because the atmosphere is H$_2$-dominated. The second follows from C/O = $N_{\mathrm{C}}/N_{\mathrm{O}}$, and the coefficients in the third equation are the squared molecular weights, e.g. $18^2=324$.

To assess the robustness of the thermochemical uncertainty framework to disequilibrium effects, we applied an idealised photochemical perturbation to the equilibrium abundances before propagating errors. This perturbation was not intended to reproduce a complete photochemical network. It provided a controlled test of how photochemistry-induced shifts in the abundances of the major observable C–H–O species modified the inferred fractional uncertainties. We approximated the net redistribution among observable species using two simplified stoichiometric proxies representative of the dominant redistribution among the observable C–H–O molecular reservoirs in warm H$_2$-dominated atmospheres. (e.g. Moses et al. 2011; Tsai et al. 2017).

The first proxy represented the net oxidation of CO to CO$_2$ driven by water photolysis. Bolded species are those which appear in the net balanced equation.

which is represented by the net stoichiometric conversion

The second proxy represented oxidation of CH$_4$ to CO through intermediate hydrocarbons and oxygen-bearing species,

The reactions gave the net conversion

In photochemical regions, H$_2$O photolysis supplies H and OH radicals that initiate or facilitate methane destruction. The net conversion treats H$_2$O as the oxygen source incorporated into CO rather than as a unique reactive intermediate. For each Monte Carlo realisation, we first computed the equilibrium abundances of CH$_4$, H$_2$O, CO, CO$_2$, and C$_2$H$_2$. We then converted a fraction $f=\{0.10,0.50,1.00\}$ of the relevant reactant pool, subject to stoichiometric limits that prevent negative abundances. R1 was limited by the available CO and H$_2$O, while R2 was limited by the available CH$_4$ and H$_2$O. We evaluated three perturbation cases: R1 alone, R2 alone, and R2 followed by R1, with the same value of $f$ applied to both reactions in the combined case. The R2 conversion increased the total number of molecules because two reactant molecules produce four product molecules. Therefore, we renormalised all molecular mixing ratios after applying R2. C$_2$H$_2$ did not participate in either proxy reaction, and the H$_2$ mixing ratio was recovered from mixing-ratio closure. The unperturbed equilibrium solution ($f=0$) provided the thermochemical reference case. We then re-evaluated Eqs. (5)–(7) using the perturbed abundances.

This procedure did not introduce an additional photochemical error term, but quantified how plausible photochemistry-induced redistribution among the observable molecular reservoirs modified the propagation of thermochemical uncertainties.

Application to published retrievals

To evaluate the analytic uncertainty framework using representative observational data, we selected four published atmospheric retrievals spanning H$_2$/He-dominated hot Jupiters and a metal-rich sub-Neptune (Table 1). Retrieved molecular abundances of H$_2$O, CO, CO$_2$, and CH$_4$ were compiled directly from the published posterior constraints. These abundances were used to calculate the corresponding O/H, C/O, and mean molecular weight using the framework described above. Uncertainties on the derived quantities were obtained by propagating the reported molecular abundance uncertainties through Eqs. (5)–(7). This comparison provided an observational test of the analytic predictions across a range of atmospheric compositions and retrieval precisions.

Parameter / species	WASP-39b	HD 209458 b	TOI-270 d	HAT-P-26 b
Model	PRISM grey	Clouds + informed CO prior	Fiducial	Best-fit
Retrieval reference	Lueber et al. [37]	Bachmann et al. [2]	Felix et al. [21]	Gressier et al. [25]
$T_{\mathrm{eff},\star}$ (K)	$5485 \pm50$	$6030$	$3506 \pm70$	$5104^{+89}_{-103}$
Mass of planet ($M_{\mathrm{Jup}}$)	$0.281 \pm0.032$	$0.682^{+0.014}_{-0.015}$	$0.0150 \pm0.0014$	$0.0585 \pm0.0069$
Radius of planet ($R_{\mathrm{Jup}}$)	$1.279 \pm0.040$	$1.359^{+0.016}_{-0.019}$	$0.1903 \pm0.0052$	$0.550 \pm0.010$
Stellar reference	Mancini et al. [41]	Rosenthal et al. [52]	Van Eylen et al. [60]	Gressier et al. [25]
Planetary reference	Mancini et al. [41]	Bonomo et al. [12]	Van Eylen et al. [60]	Wakeford et al. [62]
$\log\mathrm{H_2O}$	$-3.10^{+0.20}_{-0.19}$	$-2.80^{+0.42}_{-0.41}$	$-4.1^{+1.7}_{-5.0}$	$-2.07^{+0.39}_{-0.59}$
$\log\mathrm{CO}$	$-2.85^{+0.17}_{-0.28}$	$-4.15^{+0.51}_{-0.30}$	$< -4$	$-5.3^{+2.4}_{-4.1}$
$\log\mathrm{CO_2}$	$-4.34^{+0.24}_{-0.22}$	$-6.16^{+0.38}_{-0.39}$	$-1.29^{+0.17}_{-0.19}$	$-2.91^{+0.56}_{-0.71}$
$\log\mathrm{CH_4}$	$\ldots$	$-9.4^{+2.2}_{-2.3}$	$-1.72^{+0.28}_{-0.29}$	$< -5.49$
O/H	$0.00232^{+0.00078}_{-0.00071}$	$0.0017^{+0.0026}_{-0.0010}$	$0.103^{+0.049}_{-0.036}$	$0.0111^{+0.014}_{-0.0066}$
C/O	$0.63^{+0.36}_{-0.37}$	$0.043^{+0.099}_{-0.071}$	$0.76^{+0.41}_{-0.42}$	$0.11^{+0.32}_{-0.17}$
$\mu$ (amu)	$2.572^{+0.020}_{-0.020}$	$2.538^{+0.051}_{-0.018}$	$5.64^{+1.4}_{-0.91}$	$2.72^{+0.40}_{-0.16}$

Forward-model population synthesis

We generated population-level predictions for C/O and $\mu$ by constructing a suite of atmospheric models for four sub-populations of sub-Neptune- and hot-Jupiter-like exoplanets. The sub-populations varied in planetary class, host-star spectral type, and chemical network (Table 2).

Model group	$T_{\star}$ (K)	Spectral type	Planet type	Chemical network	Mass ($M_{\mathrm{Jup}}$)	Radius ($R_{\mathrm{Jup}}$)	$g_{\mathrm{surface}}$ (cm s$^{-2}$)
A	5772	G	Sub-Neptune	CHO	0.0186	0.196	1200
B	5772	G	Hot Jupiter	CHO	1.03	1.13	2000
C	5772	G	Sub-Neptune	CHONS	0.0186	0.196	1200
D	3473	M	Sub-Neptune	CHO	0.0186	0.196	1200

We computed these models with the one-dimensional photochemical kinetics code VULCAN (Tsai et al. 2017, 2021), which solves the coupled continuity equations for chemical species under the combined effects of thermochemistry, photolysis, and vertical mixing.

For each model, the pressure–temperature profile was calculated using the Guillot [26] formalism, following the analytical solution of Heng & Lyons [28], which allows for scattering. The irradiation temperature was calculated from the stellar temperature and semi-major axis, with the stellar radius fixed to 1 $R_{\odot}$ for all model groups. Although this value is larger than a representative M-dwarf radius, we adjusted $a$ to reproduce each target $T_{\mathrm{eq}}$. Consequently, $a$ scaled linearly with $R_{\star}$, leaving the ratio $R_{\star}/a$, and hence the stellar-flux dilution factor $(R_{\star}/a)^2$, unchanged. Therefore, the assumed stellar radius did not affect the irradiation received by the model atmosphere at fixed $T_{\mathrm{eq}}$; it only produced nominal semi-major axes larger than those expected for a typical M-dwarf radius. The longwave and shortwave opacities were fixed to 0.01 and 0.005 cm$^2$ g$^{-1}$, respectively.

The pressure–temperature profile was coupled to the relevant chemical network, as listed in Table 2, and the model was evolved to a converged steady-state abundance profile. To control the timestep size, we used an adaptive relative tolerance, setting the initial relative tolerance to 1, with lower and upper limits of 0.02 and 2.5, respectively.

Sampled parameters

To explore the effect of chemical composition, we varied the input metallicity, parameterised by O/H, and the elemental C/O ratio. The input values were generated using uniform Monte Carlo sampling. We sampled C/O uniformly between 0.1 and 1.0 and sampled O/H from a log-uniform distribution between $5 \times10^{-5}$ and $5 \times10^{-2}$. We then calculated C/H from the sampled C/O and O/H values. Finally, to incorporate the effect of the pressure–temperature structure, we varied the equilibrium temperature, $T_{\mathrm{eq}}$, from 800 K to 2500 K. Since VULCAN does not use $T_{\mathrm{eq}}$ as a direct input parameter, this variation was imposed by adjusting the semi-major axis, $a$. We calculated $a$ assuming complete heat redistribution and a Bond albedo of $A_{\mathrm{B}} = 0$. The zero-albedo assumption corresponds to complete absorption of the incident stellar radiation and provided a consistent irradiation convention across the generic model grid without introducing an additional poorly constrained parameter.

Using the three variable parameters, O/H, C/O, and $T_{\mathrm{eq}}$, we randomly generated 125 unique parameter combinations. For each combination, we constructed four classes of model: thermochemical+low-mixing, thermochemical+high-mixing, photochemical+low-mixing, and photochemical+high-mixing. The photochemical models included photolysis in the VULCAN chemical solver, whereas the thermochemical models excluded photolysis and evolved the atmosphere under thermal chemistry and vertical transport.

Vertical mixing was parameterised by a constant eddy diffusion coefficient, $K_{\mathrm{zz}}$. Since $K_{\mathrm{zz}}$ remains poorly constrained observationally, it is commonly treated as a free parameter or estimated from simplified transport prescriptions in one-dimensional photochemical models (e.g. Allen et al. [1]; Komacek et al. [34]). We adopted two constant values, $K_{\mathrm{zz}} = 10^{6}$ and $10^{10}$ cm$^2$ s$^{-1}$, which spanned four orders of magnitude and fell within the range of transport efficiencies explored in previous disequilibrium chemistry studies and inferred from global climate model (GCM)-based estimates of vertical mixing (e.g. Baxter et al. [4]; Charnay et al. [16]). The lower value was representative of weak-mixing cases used in sub-Neptune photochemical models, while the upper value was comparable to the strong-mixing cases and upper-atmosphere eddy diffusion coefficients adopted in one-dimensional hot-Jupiter and VULCAN photochemical models (e.g. Miller-Ricci Kempton et al. [43]; Tsai et al. [55]). Therefore, these two cases provided a simple framework for assessing the sensitivity of transmission-region chemistry to vertical transport while minimizing the dimensionality of the parameter space. The four model classes applied to the 125 sampled parameter combinations yielded 500 models per sub-population, and 2000 models in total.

Metric calculation

From each converged model, we calculated $\mu$ and C/O over the pressure range $10^{-3}$–$10^{-1}$ bar, corresponding to the atmospheric region typically probed by transmission spectroscopy.

The mean molecular weight was computed as

while elemental ratios were computed as

where $N_{i,\mathrm{C}}$ and $N_{i,\mathrm{O}}$ are the numbers of carbon and oxygen atoms in species $i$, respectively. These quantities were evaluated at each pressure level and averaged over the selected pressure range.

We quantified robustness by comparing matched model pairs that differed only in chemistry or in $K_{zz}$. For chemistry, we subtracted the $\mu$ of the thermochemical model from the photochemical $\mu$ value and normalised by the thermochemical value. We made this comparison separately for the low- and high-mixing pairs. For vertical mixing, we subtracted the low-mixing value from the high-mixing value and normalised by the low-mixing value, holding the chemistry fixed.

where $\Theta$ is the pressure-averaged transmission-region value of either C/O or $\mu$.

Results

Intrinsic uncertainty is lower for $\mu$ than for elemental ratios

Analytic uncertainty propagation predicts substantially larger fractional uncertainties for elemental abundance ratios than for mean molecular weight (Figure 1). The propagated uncertainties in O/H, C/O, and $\mu$ were evaluated using Eqns. (5)–(7) over the sampled $T$–$P$–C/O grid. Assuming a common per-species abundance uncertainty of $\delta= 1$, the median fractional uncertainty on C/O is approximately 1.26, while O/H exhibits a similarly large uncertainty of $\sigma_{\mathrm{O/H}}/(\mathrm{O/H}) \approx0.7$–1.0. In contrast, the propagated uncertainty on $\mu$ remains more than two orders of magnitude smaller, with a median value of $9 \times10^{-3}$, and varies only weakly across the sampled parameter space.

This difference in error propagation arises from the different roles of the trace molecular species in each derived quantity. Mean molecular weight is dominated by the H$_2$/He background, so uncertainties in the trace C–H–O species perturb only a small fraction of the total atmospheric mass. Elemental abundance ratios are constructed entirely from those same trace species, allowing molecular abundance uncertainties to propagate directly into C/O and O/H. The analytic framework also predicts different sensitivities to atmospheric metallicity. As heavy molecules comprise an increasingly large fraction of the bulk atmosphere, the uncertainty on $\mu$ increases as the dominance of the H$_2$/He background diminishes. In contrast, the fractional uncertainty on C/O is comparatively insensitive to metallicity because both the numerator and denominator scale with the abundances of the same heavy-element carriers, leaving the propagated uncertainty primarily controlled by the relative precision of the retrieved molecular abundances.

The photochemical perturbations preserve the contrast predicted by the analytic calculation (Figure 1). The vertical bars show the maximum decreases and increases relative to the baseline thermochemical calculation when reactions R1 (Eqn. (11)) and R2 (Eqn. (12)) were applied simultaneously, using the same conversion fraction, $f = 0.10$, 0.50, or 1.00, for both reactions. In all cases, the propagated uncertainties in C/O and O/H remain of order unity, whereas those on $\mu$ remain more than two orders of magnitude smaller. The photochemical bars remain within the order of magnitude of the thermochemical calculations, showing that substantial photochemistry-induced redistribution among the observable C–H–O reservoirs does not overturn the contrast between elemental ratios and mean molecular weight.

Derived atmospheric metrics are robust to disequilibrium chemistry

Photochemistry and vertical mixing do not introduce a systematic bias into the transmission-region C/O (Figure 2). Across all model populations, the transmission-region C/O closely follows the input elemental ratio under thermochemical equilibrium, photochemical equilibrium, and both weak- and strong-mixing conditions. All models lie close to the one-to-one relation, indicating that disequilibrium chemistry does not systematically shift the inferred C/O away from its equilibrium value.

The residual changes in transmission-region C/O are also small (Figure 3). Relative to the thermochemical solution, most models differ by less than 1%, irrespective of photochemistry or vertical mixing ($\gtrsim98\%$ of chemistry pairs and $\gtrsim95\%$ of mixing pairs). The few outliers occur at the extremes of the temperature and abundance ranges, where the averaged C/O becomes sensitive to small changes in a limited number of trace molecular reservoirs.

The high-temperature outliers occur in thermochemical mixing pairs with $T_{\mathrm{eq}} \gtrsim2000$ K and low-to-intermediate O/H, approximately O/H $\sim2 \times10^{-4}$–$4 \times10^{-3}$. These models reach $\Delta$C/O $\simeq-3$ to $-4\%$, with high mixing producing lower C/O at pressures probed by transmission spectroscopy than low mixing. In these thermochemical atmospheres where C and O are carried mainly by H$_2$O and CO, increasing $K_{zz}$ modestly increases the H$_2$O abundance and lowers the averaged C/O over the same pressure range. The absolute abundance changes are small, but the elemental ratio is sensitive because most of the carbon and oxygen are carried by only a few species.

A second group appeared in cool, extremely oxygen-poor photochemical model pairs, with $T_{\mathrm{eq}} \sim800$ K and O/H $\sim2 \times10^{-4}$. These models showed the same response as the hot thermochemical outliers: high mixing lowered C/O by approximately 2–3%. However, the controlling chemical reservoirs differed. Under photochemistry and high $K_{zz}$, carbon was redistributed from hydrogen cyanide (HCN), partly into carbon disulfide (CS$_2$) and methane (CH$_4$). The residual CO abundance decreased substantially, whereas the H$_2$O abundance increased only slightly. Ammonia (NH$_3$) increased concurrently as nitrogen followed the same depletion of HCN. Therefore, in this oxygen-poor regime, C/O was sensitive to the redistribution of carbon and oxygen among a small number of chemical reservoirs.

The remaining C/O outliers were confined to cool M-star models at intermediate O/H, with $T_{\mathrm{eq}} \simeq840$–880 K and O/H $\simeq3 \times10^{-3}$–$6 \times10^{-3}$. Unlike the previous groups, high mixing increased the transmission-region C/O by approximately 2–3%. In these models, carbon was redistributed from methane into benzene (C$_6$H$_6$), whereas oxygen remained primarily in H$_2$O. This redistribution increased the carbon inventory relative to the oxygen retained in H$_2$O, producing a mixing response opposite to those of the hot and oxygen-poor outliers.

Mean molecular weight is similarly robust to disequilibrium chemistry for most models (Figure 3). Across the grid, $\mu$ is essentially unchanged for $\gtrsim93\%$ of mixing pairs. The exceptions are concentrated at high oxygen abundance and cool-to-intermediate equilibrium temperature, with O/H $\simeq10^{-2}$–$5 \times10^{-2}$ and $T_{\mathrm{eq}} \simeq1100$–1550 K. In this regime, high mixing raises the transmission-region $\mu$ by approximately 2–7% relative to low mixing.

These outliers occur because the transmission region becomes sensitive to vertical redistribution among the dominant heavy molecular carriers. In low-$K_{zz}$ photochemical models, the $10^{-1}$–$10^{-3}$ bar region remains vertically stratified, with H$_2$O and CH$_4$ enhanced at higher pressures and CO more abundant at lower pressures. Strong eddy mixing homogenises this structure toward an H$_2$O- and CH$_4$-rich, CO-poor composition across the averaged pressure range, increasing the pressure-averaged $\mu$. Similar outliers in the thermochemical mixing comparisons arise from mixing-dependent quenching among the same molecular reservoirs, which also shifts the transmission region toward higher H$_2$O and CH$_4$ and lower CO.

The individual molecular abundances exhibit substantially larger disequilibrium effects than the derived atmospheric metrics. Weakly mixed atmospheres develop larger photochemical departures than strongly mixed atmospheres, yet these molecular abundance changes produce only negligible changes in C/O and $\mu$.

Published atmospheric retrievals support the analytic error-propagation pattern

The benchmark retrievals reproduce the uncertainty hierarchy predicted by the analytic framework and preserved by the disequilibrium models (Table 1). Despite spanning H$_2$/He-dominated hot Jupiters and a metal-rich sub-Neptune, all four retrievals exhibit substantially larger fractional uncertainties in C/O than in $\mu$.

The three H$_2$/He-dominated atmospheres (WASP-39b, HD 209458 b, and HAT-P-26 b) exhibit fractional uncertainties on C/O ranging from approximately 60% to more than 200%, despite molecular abundance uncertainties of tens of percent to order unity. In contrast, the corresponding uncertainties on $\mu$ remain below $\sim10\%$ and are close to 1% for WASP-39b and HD 209458 b.

TOI-270 d departs from this behaviour in the manner predicted by Sect. 3.1. Its substantially higher atmospheric metallicity increases the contribution of heavy molecules to the bulk atmospheric mass, raising the fractional uncertainty on $\mu$ to approximately 20%, while the uncertainty on C/O remains comparable to the other retrievals at approximately 55%.

Discussion

Molecular redistribution under disequilibrium chemistry

Water is the dominant radical source in warm H$_2$-rich atmospheres because it is both abundant and readily photodissociated at ultraviolet wavelengths. The resulting H and OH radicals initiate the oxidation network that controls much of the observable carbon chemistry. OH efficiently oxidises CO to CO$_2$, while H atoms drive methane oxidation through successive hydrogen-abstraction reactions. Although CH$_4$ also photodissociates in the upper atmosphere, radical-driven chemistry dominates throughout most of the transmission region. In contrast, CO remained comparatively resistant to direct photolysis because dissociation of its strong C $\equiv$ O triple bond requires substantially higher photon energies than cleavage of the O − H and C − H bonds in H$_2$O and CH$_4$, respectively: approximately 11.5 eV for C $\equiv$ O, compared with 5.9 eV for O − H and 6.7 eV for C − H.

Introducing nitrogen chemistry expands the set of observable carbon reservoirs rather than fundamentally changing the carbon chemistry. HCN provides an additional sink for carbon under reducing conditions, while N$_2$ remains chemically inert because of its strong triple bond. Sulfur follows a similar pattern, with most sulfur remaining locked in a small number of stable reservoirs (H$_2$S in the deep atmosphere; SH, S$_2$, CS$_2$ and ultimately atomic S at lower pressures under photochemistry), and contributing little to either the observable oxygen chemistry or the bulk atmospheric mass. Therefore, these additional pathways modify molecular partitioning without fundamentally altering the dominant carbon, oxygen, or hydrogen reservoirs.

Mean molecular weight and the radius valley

Our results are particularly relevant to recent proposals that atmospheric $\mu$ may vary systematically across the sub-Neptune radius valley. Heng et al. (2025) predicted such a gradient from geochemical outgassing models, with the trend controlled primarily by interior oxygen fugacity rather than atmospheric carbon enrichment. In this context, the relative robustness of $\mu$ found here is encouraging. If $\mu$ is less sensitive than elemental ratios to molecular abundance uncertainty and transmission-region photochemical redistribution, then population-level trends in $\mu$ may provide a more stable link between observable atmospheres and interior conditions. However, our models do not include interior outgassing or atmospheric escape, so they should be viewed as

testing the atmospheric robustness of $\mu$, not the origin of any radius-valley trend.

The analytic calculation assigns the same fractional uncertainty to every molecular abundance. Molecular constraints in atmospheric retrievals are not uniform: feature strength, wavelength coverage, overlapping opacity, clouds, and model assumptions can cause one gas to be constrained more precisely than another. Therefore, the common fractional uncertainty provides a controlled baseline rather than a realistic retrieval model. It isolates how the mathematical definitions of O/H, C/O, and mean molecular weight propagate a fixed level of molecular abundance uncertainty. Gas-dependent uncertainties could change the numerical contrast between these metrics, particularly when the least-constrained gas dominates an elemental inventory. The published case studies partially relax the common-uncertainty assumption by propagating the gas-dependent abundance uncertainties reported by each retrieval. These cases recover the same qualitative hierarchy between C/O and mean molecular weight. However, propagation of marginal abundance intervals does not capture covariance between retrieved gases.

Uncertainty limits on C/O as a formation tracer

The large propagated uncertainties in C/O introduce an observational caveat to its use as a formation tracer. Disk-chemistry models predicted that the C/O ratios of gas and solids vary across volatile snowlines ((Öberg et al. 2011)), while observations showed that disk C/O can also vary spatially within individual planet-forming systems (e.g. [13]; [33]). Planet-formation calculations further showed that atmospheric C/O depends not only on formation location but also on the relative amounts of gas, pebbles, and planetesimals accreted by the planet. Consequently, planets formed at different orbital radii can acquire similar atmospheric C/O values, preventing a unique formation location from being inferred from C/O alone ([18]; [10]). Our results revealed a complementary limitation: calculating C/O from uncertain molecular abundances can amplify their uncertainties, even when C/O remains chemically meaningful. Therefore, formation interpretations based on retrieved C/O require chemical and disk-evolution modelling together with careful propagation of the molecular abundance uncertainties.

Conclusion

We developed an analytic framework to quantify how uncertainties in retrieved molecular abundances propagate into the derived atmospheric metrics C/O, O/H, and mean molecular weight, $\mu$. We evaluated these uncertainties across a sampled $T$–$P$–C/O$ grid, tested the predictions against photochemical forward models spanning a range of vertical mixing efficiencies, and applied the same approach to representative published atmospheric retrievals of hot Jupiters and sub-Neptunes.

The analysis shows that elemental abundance ratios are substantially more sensitive to molecular abundance errors than mean molecular weight. Photochemical forward models show that this relative ordering is preserved under disequilibrium chemistry. Although photochemistry and vertical mixing substantially redistribute molecular abundances, they produce only minor changes in the derived atmospheric metrics. The same qualitative behaviour is recovered in published atmospheric retrievals.

Our results indicate that the principal limitation on retrieved elemental abundance ratios arises from amplified molecular abundance uncertainties rather than from photochemical redistribution of molecular reservoirs. The method presented here provides a computationally inexpensive way to estimate these uncertainties and establishes a baseline for future investigations incorporating additional atmospheric processes, including condensation, cloud formation, and higher-metallicity atmospheres.

Acknowledgements. The authors thank the anonymous referee for constructive comments that improved the clarity of the manuscript. B.B. acknowledges financial support from the European Research Council (ERC) Geoastronomy Synergy Grant (grant No. 101166936; PIs: Gaillard, Heng, and Mojzsis). The authors acknowledge the use of generative AI tools to assist with language editing and manuscript refinement. All scientific content, analysis, interpretation, and conclusions are the responsibility of the authors.