Inferring Shallow Surfaces on sub-Neptune Exoplanets with JWST

Planets smaller than Neptune and larger than Earth make up the majority of the discovered exoplanets. Those with H$_2$-rich atmospheres are prime targets for atmospheric characterization. The transition between the two main classes, super-Earths and sub-Neptunes, is not clearly understood as the rocky surface is likely not accessible to observations. Tracking several trace gases (specifically the loss of ammonia (NH$_3$) and hydrogen cyanide (HCN)) has been proposed as a proxy for the presence of a shallow surface. In this work, we revisit the proposed mechanism of nitrogen conversion in detail and find its timescale on the order of a million years. NH$_3$ exhibits dual paths converting to N$_2$ or HCN, depending on the UV radiation of the star and the stage of the system. In addition, methanol (CH$_3$OH) is identified as a robust and complementary proxy for a shallow surface. We follow the fiducial example of K2-18b with a 2D photochemical model (VULCAN) on an equatorial plane. We find a fairly uniform composition distribution below 0.1 mbar controlled by the dayside, as a result of slow chemical evolution. NH$_3$ and CH$_3$OH are concluded to be the most unambiguous proxies to infer surfaces on sub-Neptunes in the era of the James Webb Space Telescope (JWST).


INTRODUCTION
Sub-Neptune-sized planets (R p ∼ 1.6-3.5 R ⊕ ) constitute the main population of exoplanets we have discovered (Hsu et al. 2019). Yet their formation (Bean et al. 2021), atmospheric composition (Moses et al. 2013), and interior structure Aguichine et al. 2021) are not well understood. A recurring question is whether these planets are close to scaled-up terrestrial planets or scaled-down ice giants. The massradius relation, when available, provides the first assessment. However, the bulk density can be highly degenerate from different mass fractions of the H 2 /He envelope and iron/rocky core. Particularly the low-density (0.25 ρ ⊕ -0.75 ρ ⊕ ) planets can either consist of a dense core enclosed by a massive H 2 /He envelope or a lighter core (dominated by rocky material and/or water) with a shallow H 2 -rich atmosphere.
Atmospheric observations provide potential diagnostics to distinguish the above two classes. Yu et al. (2021) suggest that atmospheric chemistry can be utilized to infer the pressure level of the atmosphere-interior interface. When the atmosphere is thin (e.g., < 10 bar), photochemically unstable gases (e.g., ammonia (NH 3 ) and methane (CH 4 )) would decline without being reformed through thermochemistry in the deep atmosphere. Hu et al. (2021) focus on the distinctive features of an ocean planet with a thin atmosphere where solubility equi-librium is expected. The temperate sub-Neptune K2-18b with a recent water detection (Benneke et al. 2019;Tsiaras et al. 2019) is applied in both studies as a fiducial example. While Yu et al. (2021) do not consider atmosphere-surface exchanges, Hu et al. (2021) assume a plausible range of CO 2 concentrations in the initial atmosphere consistent with a massive water ocean. Nevertheless, both conclude that CH 4 in a 1-bar atmosphere is about 10-100 times less abundant than that in a massive atmosphere, independent of the amount of CO 2 . Yu et al. (2021) show that the slow chemical recycling makes NH 3 the most sensitive species and a prominent proxy for the presence of surfaces. However, the amount of NH 3 in a 1-bar atmosphere predicted in both studies differ by orders of magnitude due to different assumptions of initial composition. Lastly, previous studies have commonly relied on 1D models and overlooked the interaction between dayside and nightside.
In this work, we first revisit the destruction mechanism driven by photochemistry that has been suggested to serve as an observational discriminator between shallow and deep atmospheres. We elucidate the timescale and duality of the nitrogen conversion with a quiet and an active M star, and for a range of atmospheric metallicity and vertical mixing. Finally, a key addition to previous work is that we apply a 2D model including day-night transport to reevaluate the viability of utiliz-ing atmospheric chemistry as a proxy for the presence of surfaces.

1D radiative transfer model
To facilitate comparison with previous work, we take K2-18b as a test case for identifying the surface underneath a small H 2 atmosphere. The pressure-temperature (P -T ) profiles of K2-18b in dry radiative-convective equilibrium are computed by the radiative transfer model HELIOS (Malik et al. 2019a;Malik et al. 2019b). We assume a moderate 100 times solar metallicity (same as Yu et al. (2021)) and an internal heat with T int = 30 K (Lopez & Fortney 2013) with uniform heat redistribution. The P -T profile is not self-consistently evolving with chemistry but simply fixed as input to the photochemical model. We discuss the radiative effects of the surface in Section 3.1 and the radiative feedback of disequilibrium chemistry in Section 6.

1D photochemical model
The atmospheric composition is computed using the photochemical kinetics model VULCAN (Tsai et al. 2017;Tsai et al. 2021). The model treats photolysis, thermochemical kinetics, and vertical mixing. VUL-CAN has been applied to a variety of planetary atmospheres, including hot-Jupiters, sub-Neptunes, and Earth (e.g. Zilinskas et al. 2020;Tsai et al. 2021). We consider species containing N, C, O, H for this study. The large-scale mixing process is parameterized by the eddy diffusion using the expression as a function of pressure in bar (P bar ) (Tsai et al. 2021) in our nominal model We further test the sensitive to eddy diffusion and the effects of surface sinks of ammonia, as discussed in Section 3.5 and 6.

2D photochemical model
The day-night circulation is expected to be crucial in regulating the compositional variation across a tidallylocked planet (e.g. Agúndez et al. 2014;Drummond et al. 2020;Feng et al. 2020;Wardenier et al. 2021). We run the 2D photochemical model VULCAN on a pressurelongitude grid to account for the day-night transport on a meridionally averaged equatorial plane. The 2D model solves the continuity equations including horizontal transport where n i is the number density (cm −3 ) of species i, P i and L i are the chemical production and loss rates (cm −3 s −1 ) of species i, and φ i,z , φ i,x are the vertical and horizontal transport flux, respectively. The construction and benchmarks of 2D VULCAN are discussed in detail in a follow-up paper (in prep

Transmission spectrum modeling
We use the radiative transfer and retrieval framework NEMESIS (Irwin et al. 2008) to produce the transmission spectra. The forward models were computed using the correlated-k technique, with the k-tables being computed using the methodology of Chubb et al. (2021). Specifically, the sources of the opacity data for the molecules of interest are: NH 3 (Coles et al. 2019), CO (Li et al. 2015), CO 2 (Yurchenko et al. 2020), CH 3 OH 1 (Gordon & et al. 2017).
3. 1D RESULTS: SLOW LOSS OF AMMONIA 3.1. Minimal surface effects on the thermal structure Before looking into the compositional diagnostics, we first discuss how the existence of surfaces might impact the thermal structure. The radiative effect of the surface is neglected in Yu et al. (2021) where the temperature profile is simply truncated at different pressure levels. In practice, the surface partially absorbs stellar irradiation and re-emits back to the atmosphere. A canonical example is that of Earth, where roughly half of the solar radiation is absorbed by the surface but the atmosphere is mostly opaque to infrared radiation. This potentially makes the temperature higher than that of a "surfaceless" atmosphere at the same pressure level, which can in turn alter the composition and the ability to recycle atmospheric species.
We test the thermal effects of the surface using 1D radiative-convective calculations with a generic surface placed at different pressure levels. We assume a surface albedo of 0.1 but find no differences in the range of 0.1-0.3 (corresponding to the albedo range of land areas on Earth) 2 . The resulting temperature profiles for 100 times solar metallicity are depicted in Figure 1. The convective zone extends from a few bar to 0.1 bar and the temperature is set by the dry adiabat. The presence of a surface turns out to have negligible thermal effects in this region since the atmosphere absorbs most of the stellar irradiation. Only for surface pressures around 0.01 bar, the greenhouse warming near the surface starts to be notable, because the atmosphere is optically thin toward stellar radiation while being thermally opaque in this region. In all cases, our test validates the approach of directly truncating the model for shallow atmospheres with surface pressure down to about 1 bar (also see May & Rauscher (2020) for the global heat transport).

Nitrogen conversion for a quiet M star and an active M star
Photochemically active gases, including NH 3 , HCN, and CH 4 , are lost in the absence of a deep atmosphere where themochemistry operates fast and fully recycles them back. Comparing the results of a surface at 1 bar to the deep/no surface case, Yu et al. (2021) find that NH 3 is decreased by about 5 orders of magnitude while HCN is decreased by about 3 orders of magnitude. In this section, we will address the timescale of nitrogen conversion and overall compositional evolution driven by photochemistry. It is relevant to consider the host star at different stages when the timescale of chemical destruction is comparable to the timescale of stellar evolution. We apply the synthetic M2 star at the age of 45 Myr and 5 Gyr from HAZMAT (Peacock et al. 2020) for the stellar UV spectra, to represent a young, active and an old, quiet M star 3 . Figure 2 shows snapshots of NH 3 and HCN abundances at different times with a 1-bar surface compared with a deep, 1000-bar surface, for both a quiet M star and an active M star 4 . All models start from initial abundances in chemical equilibrium. We first confirm that while recycling is enabled in the deep surface scenario, NH 3 and HCN settle to a steady state within 1000 years. In the 1-bar atmosphere, NH 3 mostly drops off between 10 kyr and 1 Myr to a final mixing ratio several orders of magnitudes lower, qualitatively matching the decrease of NH 3 in Yu et al. (2021). However, as NH 3 is gradually lost, HCN does not readily follow NH 3 . With a quiet M star, HCN actually increases to a uniform abundance close to 1% in the period between 0.1 and 10 Myr before gradually declining after 10 Myr.
The timescale of nitrogen conversion given by the evolution of major nitrogen species is further captured in Figure 3 (a). Overall, we find that nitrogen conversion takes place between 10 kyr and 1 Myr for an active M star and is much slower for a quiet M star. NH 3 is evidently converted to N 2 for an active M star, but HCN can serve as a transient nitrogen pool in the case of a quiet M star.
3.2.1. The two conversion paths of NH3 HCN is produced in the upper atmosphere initiated by NH 3 photolysis: While HCN is subject to photodissociation into cyanide radical (CN) and H, CN can rapidly recycle back to HCN with CN + H 2 −−→ HCN + H in an H 2 atmosphere. In fact, the major effective sink of HCN is the OH radical produced by water photolysis. HCN is consumed by OH as C irreversibly goes into CO and N converted into N 2 . The primary scheme for the loss of NH 3 and HCN 10 3 10 4 10 5 10 6 10 7 10 8 Time ( The critical difference between a quiet and an active M star is that an active M star produces orders-ofmagnitude more OH radicals and hence HCN is scavenged more efficiently after building up. On the other hand, HCN is gradually consumed by OH only after 10 Myr with a quiet M star. We can estimate the overall nitrogen conversion timescale of scheme (4) with its rate-limiting step (Tsai et al. 2018) where k is the rate coefficient for the reaction NH 2 + NO −−→ N 2 + H 2 O. Taking the initial abundances at 0.1 mbar and 1000 years, the timescale of NH 3 by (5) is ∼ 0.1 Myr for an active M star and ∼ 10 Gyr for a quiet M star, consistent with the kinetics results. This is not to be confused with the timescale of NH 3 photodissociation, which operates much faster, about a few hours at the same pressure level. We conclude this section by reiterating that although the ammonia loss in the shallow-surface case is driven by photochemistry, the process of nitrogen conversion occurs over a timescale of a million years or longer for quiet M stars. It resembles denitrification processes on a geological timescale but only takes place in the atmosphere rather than shifting between reservoirs.

Carbon conversion
Most of the carbon is initially bound in the photochemically unstable CH 4 and converted into CO and CO 2 over time. The same evolution of carbon species is shown in Figure 3 (b). First, carbon conversion takes about the same timescale as the nitrogen conversion, O(Myr). Second, the conversion of CH 4 to CO and CO 2 is initiated by the CH 4 photodissociation channel: CH 4 hν −−→ CH + H 2 + H. CH reacts with H 2 O to form H 2 CO, which is readily photodissociated into CO. This conversion proceeds faster with an active M star but the final abundances of CO and CO 2 remain close independent of the host star. Third, CH 4 is still continuously evolving after millions of years with a quiet M star, making it ambiguous to compare with the deep-atmosphere abundance as a proxy for surface pressure. Last and most interestingly, methanol (CH 3 OH) is produced as a by-product of carbon conversion. A small part of H 2 CO, instead of being photodissociated, can thermally react with hydrogen to form CH 3 O and CH 3 OH. The mixing ratio of CH 3 OH is increased to ∼ 10 −6 , compared to 10 −9 and 10 −8 in the deep surface scenario for a quiet and an active M star, respectively, where CH 3 OH is transported and destroyed in the deep atmosphere.

Sensitivity to metallicity
Sub-Neptunes can in general have a wide range of atmospheric metallicity (Moses et al. 2013;Fortney et al. 2013). The equilibrium abundance of NH 3 decreases with increasing metallicity as a result of hydrogen deficiency. To test the sensitivity of NH 3 conversion to metallicity, we further explore lower metallicities with the same P -T structure. Figure 3 (c) illustrates the evolution of major nitrogen species for 1×, 10×, and 100× solar metallicity, with a 1-bar surface and an active M star, where the abundances are normalized to those in equilibrium for a clearer comparison between different metallicities. The final [NH 3 ]/[NH 3 ] EQ slightly de-creases with smaller metallicity while [HCN]/[HCN] EQ is almost identical for all metallicities. We find that the trends of NH 3 and HCN destruction remain robust. The O(Myr) timescale of nitrogen conversion holds across different metallicities as well.

Sensitivity to vertical mixing
The eddy diffusion profile in our nominal model prescribed by Equation (1) is close to that in Yu et al. (2021). The inverse square root of pressure dependence of K zz is commonly assumed to reflect the turbulent mixing due to gravity wave breaking in stratified atmospheres (e.g., Lindzen 1981; Moses et al. 2016). Since the eddy diffusion cannot be derived from first principles and is often treated as a free parameter, we further test the sensitivity to the strength and structure of vertical mixing by running the 1-bar surface case for a range of uniform K zz profiles from 10 5 to 10 9 cm 2 /s. This range of eddy diffusion coefficient is consistent with that derived from the average vertical wind in the GCM of K2-18b (see Figure 7) and other sub-Neptunes in a similar dynamical regime (e.g., Charnay et al. (2015)).
The mixing ratios of NH 3 , HCN, and CH 3 OH for different assumptions of eddy diffusion coefficient are reported in Figure 3 (d). NH 3 in the deep model retains a uniform abundance ( Figure 2) and does not depend on the eddy diffusion, whereas NH 3 decrease in the shallow model generally correlates with smaller eddy diffusion, since the lower well-mixed region leads to deeper UV penetration and photochemical destruction. Overall, NH 3 always has lower abundances in the shallow model than those in the deep model but the exact ratio, [NH 3 ] 1bar /[NH 3 ] 1kbar , depends on the strength and shape of K zz . On the other hand, CH 3 OH consistently shows little sensitivity to eddy diffusion and about two orders of magnitude increase with respect to the deep model. Conversely, HCN highly depends on eddy diffusion as a result of the balance between photochemical production/destruction in the upper atmosphere and vertical transport. HCN can either raise or lower in the presence of a shallow surface, making it unsuitable as a proxy for the surface level.

2D RESULTS: GLOBALLY UNIFORM DISTRIBUTION CONTROLLED BY THE DAYSIDE
We have established the long, O(Myr) timescale of the nitrogen conversion in Section 3. On a tidally locked exoplanet, the conversion driven by photochemistry only occurs on the dayside but completely shuts down on the permanent nightside. Without photolysis, NH 3 is able to maintain a high abundance on the nightside. To answer the question how the chemical evolution with a shallow surface would be regulated by the global circulation between the dayside and nightside, we apply a 2D photochemical model on a meridionally averaged equatorial plane to address the effects of day-night transport. The temperature and wind for the 2D photochemical model are adopted from the 3D GCM output of K2-18b, whose global structure can be found in Figure 7. We employ four quarters (Figure 7 (d)) in our 2D model for clarity, while the full 2D chemical output is included in the supplementary figure. To isolate the effects of horizontal transport, Figure 4 compares the steady-state abundances with a 1-bar surface and an active M star to those without including the zonal wind.
It is evident that without a recycling mechanism, the horizontal transport is able to exhaust NH 3 on the nightside. Even the modest zonal wind (∼ 1 m/s) yields a horizontal transport timescale less than a year, orders of magnitude faster than the chemical evolution. As a result, most species are homogenized below 1 mbar and exhibit abundances close to those predicted by the 1D model without day-night transport. The model for a quiet M star shows qualitatively consistent results (Figure 6), with reduced photochemical destruction of NH 3 and higher abundance of HCN. The results we find for the fiducial example of K2-18b can be generalized to most temperate, tidally-locked planets with a chemically inactive surface, since the equatorial jet is a robust dynamical property (e.g. Pierrehumbert & Hammond 2019).

SPECTRAL INDICATOR
Since the 1D model for the dayside is verified to represent the globally uniform distribution, we compare the transmission spectra from our deep-atmosphere and shallow-surface models in Figure 5. The most distinctive differences are CO 2 absorption features at 4.2 -5 µm and the lack of NH 3 at 3, 8.8 -12 µm for the shallowsurface case. Interestingly, CH 3 OH absorbs between 9 and 10 µm, well within the NH 3 band. To quantify the ability of detecting these molecules with JWST, we used Pandexo (Batalha et al. 2017) to generate the required observations to disentangle between our two models. We find that with it is possible to identify the NH 3 feature at 3 µm with only 3 transit observations using NIRSpec. When observing with MIRI, it is possible to further discriminate CH 3 OH from NH 3 with around 20 transits. Conversely, HCN is not abundant enough in the region sensitive to transmission to show observable features. (surface) boundary for all species. Regarding NH 3 , agricultural ammonia on Earth is absorbed by surfaces with a fairly short residence time (∼ days, e.g., Seinfeld & Pandis 2016;Jia et al. 2016). We have tested our model with Earth surface conditions assuming biotic processes are involved, and even the slowest dry deposition of NH 3 on a desert (v dep = 0.0002 cm/s, Jia et al. 2016) acts as a more efficient sink than the photochemical sink (Section 3). However, without the biosphere transferring ammonia to different organic nitrogen, it is conceivable that ammonia returns to the atmosphere at once and makes effectively net zero deposition. In general, various geological processes can be crucial in controlling trace gases in the atmosphere, especially determining the redox state of the atmosphere. For example, Wogan et al. (2020) find volcanic outgassing unlikely to be reducing, while Lichtenberg (2021) demonstrate that the interior of sub-Neptunes can remain reducing due to vigorous internal convection and Zahnle et al. (2020) suggest iron-rich impacts can also generate reducing atmospheres that favor NH 3 and CH 4 . In the case of an atmosphere with scattering clouds/hazes or less irradiated than K2-18b Blain, D. et al. 2021), the surface temperature can be brought down to suppress water evaporation (Scheucher et al. 2020) and allow liquid water oceans, which can participate in regulating the inventory of soluble gases (e.g., NH 3 , CO 2 , and CH 3 OH). With mobile lid tectonics (Tackley et al. 2013;Meier et al. 2021), the long-term evolution of CO 2 in the shallow atmosphere case is expected to be governed by outgassing and weathering (e.g. Pierrehumbert 2010).
Since CO 2 generally has other geological sources and HCN is unsuitable as a surface proxy for its strong dependence on the stellar type and vertical mixing, we propose including CH 3 OH as a complementary indicator along with NH 3 . On Earth, CH 3 OH is mainly produced by plants and deposited to the ocean (Yang et al. 2013), whereas in a H 2 -dominated atmosphere with a shallow surface, CH 3 OH is produced during the process of CH 4 oxidation. Detection of NH 3 without CH 3 OH is consis-tent with the deep/no-surface scenario, while detection of CH 3 OH but without NH 3 indicates the presence of a surface. Non-detection of both NH 3 and CH 3 OH can imply a global water-ocean world (Hu et al. 2021) as they are highly soluble in water. We summarize the flowchart of inferring the surface property of sub-Neptunes with H 2 -dominated atmospheres using JWST in the bottom panel of Figure 5.

SUMMARY
1D modeling of sub-Neptune exoplanets suggests that NH 3 is depleted below detectable level on the dayside in the presence of a shallow surface, while HCN can either increase or decrease depending on vertical mixing and quiet/active M stars. CH 3 OH is found to be consistently enhanced for all shallow-surface simulations. We construct a 2D photochemical framework to account for the day-night circulation and find the global abundance is overall quenched from the dayside as the chemical conversion takes O(Myr). Our results suggest that HCN is not applicable to determine the pressure level of the surface. Instead, the shallow-surface scenario can be ruled out by NH 3 detection, whereas positive detection of CH 3 OH with negative detection of NH 3 indicate shallow and dry surfaces.

ACKNOWLEDGMENTS
S.-M.T. acknowledges support from the European community through the ERC advanced grant EXO-CONDENSE (#740963; PI: R.T. Pierrehumbert). T.L. has been supported by the Simons Foundation (SCOL award #611576). This work has made use of the synthetic spectra from the HAZMAT program; doi: 10.17909/t9-j6bz-5g89. We thank J. Moses for thoughtful discussion.

A. RADIATIVE FEEDBACK FROM DISEQUILIBRIUM CHEMISTRY
Since the temperature structure is fixed by thermochemical equilibrium, we have also checked the radiative feedback as a result of disequilibrium chemistry. We re-run the 1D radiative-convective calculations with the final disequilibrium composition from Section 3. We find the resulting temperature from disequilibrium chemistry can be about 100 K lower at most in the convective zone above 1 bar, mainly from the decrease of opacity-predominating water. We have performed sensitivity tests and found the temperature difference does not alter the presented results (also see the thermal structure variance explored in Yu et al. (2021)). We will leave the effects of water condensation for future work.

B. RUNNING 2D PHOTOCHEMICAL MODELS FOR MILLIONS YEARS
The horizontal transport flux by the zonal wind in Equation (2) yields a hyperbolic partial differential equation. The time step is ought to follow the same format as Courant-Friedrichs-Lewy (CFL) condition when numerically integrating the system. In other words, the time step must not exceed the time that zonal flow travels across adjacent vertical columns. For K2-18b, the time step limit is ∆x u ∼ 10 5 s, which means more than 10 8 integration steps are required to integrate the system to a million year for the long-term chemical evolution. To circumvent this computational load, we arbitrarily slow down the zonal wind (e.g., by 1000 times) such that a larger time step can be adopted. Once the chemistry has evolved after the Myr integration, we recover the correct wind velocity and run the model to final steady state. This seemingly risky approach can be justified by the timescale argument: The horizontal transport effectively interacts with vertical mixing and chemical evolution , which manifest drastically different timescales. The timescale of vertical mixing is within hours and the chemical evolution takes ∼ Myr, whereas the nominal horizontal transport has a timescale of ∼ days. Therefore, the long-term compositional evolution is expected to be qualitatively unaffected for any horizontal transport much slower than vertical mixing but orders of magnitude faster than the chemical evolution. That is, the tuned-down zonal wind still plays the role of passing on the chemical evolution from dayside to nightside. Finally, we have confirmed this approach by comparing the simulation at about 5000 years to that with 1000 times slower zonal wind and found no differences. 3) temperature (black), zonal wind (solid), and eddy diffusion profiles (dashed) for the four quarters of K2-18b. All temperature profiles in black lines and overlapped because the temperature is globally uniform, whereas the zonal wind and eddy diffusion (derived from the vertical wind) profiles are color coded for the dayside, evening, nightside, and morning quarter.