Role of Planetary Radius on Atmospheric Escape of Rocky Exoplanets

Large-scale characterization of exoplanetary atmospheres is on the horizon, thereby making it possible in the future to extract their statistical properties. In this context, by using a well-validated model in the solar system, we carry out 3D magnetohydrodynamic simulations to compute nonthermal atmospheric ion escape rates of unmagnetized rocky exoplanets as a function of their radius based on fixed stellar radiation and wind conditions. We find that the atmospheric escape rate is, unexpectedly and strikingly, a nonmonotonic function of the planetary radius R and that it evinces a maximum at R ∼ 0.7 R ⊕. This novel nonmonotonic behavior may arise from an intricate trade-off between the cross-sectional area of a planet (which increases with size, boosting escape rates) and its associated escape velocity (which also increases with size but diminishes escape rates). Our results could guide forthcoming observations because worlds with certain values of R (such as R ∼ 0.7 R ⊕) might exhibit comparatively higher escape rates when all other factors are constant.


INTRODUCTION
The number of confirmed exoplanets has grown explosively in the past decade (Winn & Fabrycky 2015;Perryman 2018;Zhu & Dong 2021), with the current total exceeding 5000. 1 Current and forthcoming surveys are anticipated to raise the number of detected and characterized exoplanets by a substantive degree.As a result of the relatively large sample size, it has become feasible to extract exoplanetary statistics and trends.
In a similar vein, statistics of terrestrial planetsroughly interpreted to possess radii < 2 R ⊕ -situated within or close to the habitable zone (HZ) are likely to become increasingly relevant in the coming decades; the limits of the HZ were computed in Kasting et al. (1993) and Kopparapu et al. (2013Kopparapu et al. ( , 2014)).Rocky planets in the HZ can host liquid water on their surfaces, and such worlds are likely to be detected and characterized with increasing frequency by future telescopes (Fujii et al. 2018;Madhusudhan 2019;Wordsworth & Kreidberg 2022).The atmospheres of these terrestrial exoplanets may yield a wealth of data about their physical, chemical, and geological properties (Seager & Deming 2010;Madhusudhan 2019;Zhang 2020;TRAPPIST-1 JWST Community Initiative et al. 2023), as well as potentially even biological features (Fujii et al. 2018).
Given the centrality of exoplanetary atmospheres, especially those of rocky exoplanets proximal to the HZ as delineated above, it is of crucial importance to understand what statistical trends, if any, may be discernible by future surveys.This era has already commenced, with JWST having obtained the transmission spectrum or detected the thermal emission of some warm rocky exoplanets (Lustig-Yaeger et al. 2023;Greene et al. 2023;Zieba et al. 2023).From the perspective of atmospheric loss -which could sculpt the atmospheres of terrestrial planets (Tian 2015;Owen 2019) -the significance of nonthermal atmospheric ion escape processes driven by the solar (or stellar) wind for this class of planets is well-established in our solar system (Lammer et al. 2008;Lammer 2013;Brain et al. 2016;Dong et al. 2018b;Persson et al. 2020;Ramstad & Barabash 2021;Lichtenegger et al. 2022) and is plausible for exoplanets (as reviewed in Zahnle & Catling 2017;Lingam & Loeb 2019, 2021;Airapetian et al. 2020;Gronoff et al. 2020).
In this work, we carry out a systematic investigation of "nonthermal" atmospheric ion escape -in which the speeds of escaping particles are essentially decoupled from the exobase temperature, and often entail the presence of ions in electromagnetic fields (Tian 2015;Laneuville et al. 2020) -as a function of radius for unmagnetized rocky exoplanets with atmospheric composition akin to Venus.The latter world is recognized as a valuable benchmark for exoplanetary science (Kane et al. 2019;Kane 2022), especially insofar as planets close to the inner edge of the HZ are concerned (Ostberg et al. 2023).The outline of the Letter is as follows.We describe the multi-species MHD model and the setup in Section 2. Next, we present the results and analyze the ensuing implications in Section 3. We round off by summarizing our findings in Section 4.

METHODS AND MODEL SETUP
We leverage the 3-D Block-Adaptive-Tree Solarwind Roe-type Upwind Scheme (BATS-R-US) multispecies magnetohydrodynamic (MS-MHD) model that has been successfully implemented to model Venus-like exoplanets (Dong et al. 2017b(Dong et al. , 2018a(Dong et al. , 2019(Dong et al. , 2020) ) for simulating atmospheric ion escape from exoplanets spanning a range of radii.The MS-MHD model includes four ion species H + , O + , O + 2 , CO + 2 , and the associated ionospheric photochemistry such as photoionization, charge exchange, and electron recombination.For unmagnetized planets, the BATS-R-US MS-MHD model investigates the interactions of magnetized stellar winds with planetary upper atmospheres, and properly accounts for atmospheric escape mechanisms mediated by the stellar wind via mass loading.
In addition to assuming Venusian atmospheric composition (i.e., CO 2 -dominated atmosphere) due to the reason outlined in the last paragraph of Section 1, the model accepts several input parameters to specify planetary and stellar properties, including but not limited to planetary radius and mass, atmospheric profile, as well as stellar winds and radiation.We have chosen to work with the Planet-Star-Orbital (PSO) coordinates, where the X-axis points from the planet toward the star, the Z-axis is perpendicular to the planetary orbital plane, and the Y-axis completes the right-hand system.
We modeled eight exoplanets ranging in size from 0.50 to 2.25 Venus radii (R V ), as illustrated in Given that RV ≈ R⊕, the sizes of these putative planets span the range of around half to twice the size of Earth.(Ehlmann et al. 2016).This upper bound is compatible with empirical and theoretical constraints on the size of the largest worlds with rocky composition (Otegi et al. 2020;Plotnykov & Valencia 2020).Each model planet has a surface pressure of 1 bar; accordingly, larger planets possess more massive atmospheres, which is consistent with intuition.While we could opt to vary the surface pressure, the previous work (Dong et al. 2017b) has demonstrated that the nonthermal atmospheric ion escape rates do not significantly rely on planetary surface pressure; in other words, the choice of surface pressure is not expected to conspicuously affect our calculations on the atmospheric ion escape rates.
For each model planet, we modify the neutral atmospheric scale height (H) by varying the gravitational acceleration (g) since the former is defined as where T is the temperature, k is the Boltzmann constant, and m is the mass of the atmospheric species.Assuming that the planetary radius R is known, the mass of the planet M is derived from the simplified massradius relation for rocky planets (Valencia et al. 2006;Zeng et al. 2016Zeng et al. , 2019; see also Chen & Kipping 2017): where M V and R V are the mass and radius of Venus, respectively.Given M and R, it is found that g obeys In our study, the model exoplanets are assumed to be unmagnetized.This assumption not only matches the current state of Venus, but also might be characteristic of ∼ 50% of all rocky planets in the HZ as per some estimates (McIntyre et al. 2019).However, future work is needed to assess the impact of a planetary magnetic field on nonthermal ion escape (Dong et al. 2017b;Gunell et al. 2018;Egan et al. 2019;Lingam 2019;Dong et al. 2020), as those effects are not yet fully understood.Next, a crucial set of parameters that needs to be specified is that of the stellar wind.In this context, it should be borne in mind that M-dwarfs not only constitute the majority of stars in the Milky Way, but also their (rocky) planets represent promising targets for atmospheric characterization (Shields et al. 2016;Fujii et al. 2018;Lingam & Loeb 2021).Although we do not simulate the conditions experienced by M-dwarf exoplanets (which are diverse) as such, we instead choose to work with ancient solar wind conditions (near Earth) at ∼ 4 Ga consistent with high mass-loss rate of the young Sun (Wood et al. 2021) -with solar wind density n sw = 136.7 cm −3 , solar wind velocity v sw = (-910, 0, 0) km/s, and interplanetary magnetic field B sw = (-15.6,30.2, 0) nT -and a stellar EUV flux that is 12 times that of modern Sun (Boesswetter et al. 2010;Dong et al. 2017a) for our model exoplanets.To an extent, these enhanced stellar wind parameters and EUV flux are reminiscent of those encountered by M-dwarf exoplanets in the HZ (France et al. 2013(France et al. , 2016;;Airapetian et al. 2020).

RESULTS AND DISCUSSION
In this section, we cover the salient results and provide a plausible explanation for the discerned trends.

Results
For the model exoplanets listed in Table 1, we deployed the BATS-R-US MS-MHD model and computed the associated atmospheric ion escape rates; the latter were calculated based on the stellar radiation and wind conditions motivated toward the end of Section 2. The atmospheric ion escape rates are reported in Figure 1 and Table 2 accordingly illustrates the relationship between planetary radius and atmospheric ion escape rate.On inspecting Figure 1, it is apparent that the relatively light ion species O + has the highest ion escape rate and the relatively heavy species CO + 2 has the lowest ion escape rate.This trend is potentially explainable by the fact that lighter ions like O + are typically more abundant at higher altitudes (as reviewed in Mendillo et al. 2018), in the absence of strong (e.g., turbulent) mixing, and are consequently more vulnerable to stellar wind erosion.Hence, given that the stellar wind could strip a relatively higher concentration of lighter ions, it cannot also strip heavier ions at lower altitudes with the same efficiency due to the constraint imposed by the conservation of momentum and energy.
In examining Figure 1, we identify an intriguing nonmonotonic feature in the trend of the ion escape rate as a function of the radius (of rocky planets).For the chosen stellar wind conditions, which are loosely reminiscent of those parameters that may be experienced by temperate M-dwarf exoplanets, total ion escape rate peaks at R = 0.75R V .The maximum total ion escape rate at R = 0.75R V is about eight times higher than the minimum escape rate at R = 2.25R V (see Table 2).The putative rationale for the former feature and its implications are discussed shortly hereafter in Section 3.2.
Figure 2 depicts O + ion density contour plots for all model exoplanets subjected to the same stellar radiation and wind conditions.The black sphere in the middle represents the planetary body.The planetary radius increases from R = 0.50R V (top left) to R = 2.25R V (bottom right) with ∆R = 0.25R V .It is clear that the escaping O + ion density around the planetary body with R = 0.5R V has the highest value, but the escape rate peaks at R = 0.75R V (see Table 2).It is noteworthy that the escape rates in Figure 1 and Table 2 are calculated based on the integration of escape ion flux over a spherical surface around the planet, thus the size of the planet is also important.When planetary radius exceeds R = 1.25RV , one can barely observe the red contours around the planet due to the increasing surface gravity.One of the most striking results manifested is the nonmonotonic variation of the atmospheric escape rate with the planetary radius R. We begin by furnishing a possible qualitative explanation for this behavior, and then elaborate further on this theme.

Analysis
As we are studying atmospheric escape, there are two competing factors at play.Depending on which of the duo is more dominant, the escape rate may be altered accordingly.On the one hand, the cross-sectional area of the planet that interacts with the stellar wind is of importance.In the case of weakly magnetized and wholly unmagnetized planets, this area is roughly proportional to R 2 (cf.Salz et al. 2016), although for magnetized planets the expression is more complex (Blackman & Tarduno 2018;Lingam 2019).Therefore, as the radius is elevated, the cross-sectional area increases, and so does the escape rate because the latter scales with the area.Likewise, it is important to recognize that the area of the atmosphere -the source from which the particles escape the planet -is also proportional to R 2 .
On the other hand, as the radius is increased, the escape velocity of the planet (v e ) is also enhanced: where we have used the mass-radius relationship M ∝ R 3.7 for rocky planets (Zeng et al. 2016).On account of the higher escape velocity (at larger R), it would be harder for particles to escape the planet's gravitational well.In such a scenario, the atmospheric escape rate is presumed to steeply decline with the radius, which runs counter to the trend in the preceding paragraph.Hence, when these two aspects are duly taken into consideration, the nonmonotonic behavior of the atmospheric escape rate as a function of the planetary radius exhibited by our simulations might be explainable.
It is possible to build on this theme.For an unmagnetized planet with a fixed radius, its cross-sectional area πR 2 determines the amount of stellar radiation and wind energy incident on the planet.The larger the crosssectional area is, the more energy the planet ought to receive.The stellar radiation energy can lead to photoionization (the main source of the ionized gases) and heating of the upper atmosphere, and subsequently, the stellar wind could strip those ionized particles from the planet.Given that it is relatively easy for ionized gas to escape from a small unmagnetized planet (with lower gravity), the atmospheric ion escape may become source limited, i.e., almost all ions supplied to the region energized by the stellar wind can gain sufficient energy to escape, and therefore, the source of the ions (i.e., controlled by the ion production rate) could limit the supply and thence the total ion escape flux, rather than the amount of available energy.Therefore, this scenario enables a trend whereby larger unmagnetized planets exhibit higher atmospheric ion escape rates due to bigger reservoirs of ionized gas in their upper atmospheres.
On the other hand, as the planet's size is further increased, the escape velocity is also enhanced, as revealed by (4).Therefore, beyond a certain threshold, it would no longer be easy for atmospheric ions to attain the requisite energy for escaping the planet; to put it differently, the escape rate might transition from a source-limited regime to an energy-limited regime.A brief segue is warranted at this stage.The energy-limited regime is encountered often in hydrodynamic escape (e.g., Watson et al. 1981;Erkaev et al. 2007;Owen & Jackson 2012;Owen & Alvarez 2016;Salz et al. 2016;Kubyshkina et al. 2018;Krenn et al. 2021;Lampón et al. 2021).However, in contrast to hydrodynamic escape on worlds with H/He atmospheres, it is worth recognizing that the adopted EUV flux alone is unlikely to drive atmospheric mass loss on Venus-like planets for two major reasons: (1) the atmosphere is chiefly composed of heavy chemical species; and (2) CO 2 15-µm cooling is very efficient in the upper atmosphere of Venus-like planets (Bougher et al. 1994).In this paper, the effects of EUV flux are, instead, to produce photoionized gas and heat the upper atmosphere (through photoionization and photoelectron heating); subsequently, the stellar wind can strip the ionized gas away (Laneuville et al. 2020).
Circling back to the above theme, the bottleneck on larger rocky planets is now envisioned to be imposed by the escape velocity, which is a monotonically increasing function of R -refer to (4).Hence, we would expect the escape rate to strongly decline with an increase in the radius across this range, as confirmed by Figure 1 and Table 2. On combining these predictions, the escape rate ought to first increase and then decrease with the radius, which is consistent with Figure 1 and Table 2.
To sum up, we have argued that we could witness a change in the regime, going from a limit set by source particles to one enforced by energy.Quite strikingly, such a shift has been empirically documented in the solar system, where nonthermal ion escape on Mars is sourcelimited, whereas the equivalent on Venus is energylimited (Persson et al. 2020;Ramstad & Barabash 2021); note that both these planets are weakly magnetized, allowing for a comparison of similar worlds.This data gives credence to the preceding paragraphs, although we caution that we have offered a simplified analysis, which was necessitated by the complexity of nonthermal escape mechanisms (Cravens et al. 2017;Lingam & Loeb 2019;Airapetian et al. 2020;Gronoff et al. 2020).
In closing, we comment on the ramifications of our modeling.If the nonthermal ion escape rates do indeed peak at R ∼ 0.7 R ⊕ , this datum could have crucial consequences for future observations.Given that the computed escape rate is boosted by a factor of about 2 as we transition from R ∼ 0.5 R ⊕ to R ∼ 0.7 R ⊕ and declines thereafter by nearly the same amount when we move to R > R ⊕ , it is possible that temperate rocky planets with R ∼ 0.7 R ⊕ might harbor relatively thinner atmospheres if all other variables (e.g., outgassing, atmospheric mass) are statistically similar across the range of planetary radii.2If future telescopes such as ARIEL (Tinetti et al. 2018) characterize a sufficiently large sample of terrestrial exoplanets, it may be feasible to test some aspects of our hypothesis empirically.
4. CONCLUSIONS Building on the ongoing characterization of exoplanetary atmospheres by JWST (Gardner et al. 2023)which has led to several interesting discoveries (e.g., Lustig-Yaeger et al. 2023;Greene et al. 2023;Zieba et al. 2023) -and observations by forthcoming ground-based extremely large telescopes, a wealth of data will become available, making it feasible to discern statistical patterns in exoplanetary atmospheres in the future.Nonthermal ion escape facilitated by stellar EUV radiation and driven by the stellar wind represents one of the key regulators of the mass and composition of terrestrial exoplanetary atmospheres (Brain et al. 2016;Zahnle & Catling 2017;Airapetian et al. 2020;Gronoff et al. 2020).Motivated by these facts, we utilized a thoroughly validated multispecies MHD model to investigate the total ion escape rate as a function of the planetary radius R.
After performing the simulations, we unearthed a novel trend for rocky planets with CO 2 -dominated atmospheres: the escape rate is a nonmonotonic function of R, as clearly depicted in Figure 1, with a peak at R ∼ 0.7 R ⊕ .It was found that this characteristic is manifested for intense stellar wind and radiation conditions, which may be associated with M-dwarfs and young Sunlike stars.This nonmonotonic feature runs counter to the naive expectation that smaller rocky planets (e.g., akin to Mars) would automatically exhibit higher escape rates because of their weaker gravity.
The nonmonotonic behavior may arise from a tradeoff between the (cross-sectional) area and the escape velocity of a planet; both of these variables increase for larger planets.On the one hand, when the aforementioned area is boosted, this will enhance the atmospheric ion escape rate since the planet would intercept more of the stellar radiation and wind (due to its greater area), thereupon amplifying the escape rate on account of the enhanced ionized reservoir in the planetary upper atmosphere.On the other hand, when the escape velocity is boosted, this factor will suppress the escape rate because particles require more energy to escape the gravitational well.From a more detailed perspective, we suggested in Section 3.2 that nonthermal atmospheric ion escape may shift from a source-limited regime to an energy-limited one, which could partially explain the simulated results.
If the striking trend displayed by our modeling is correct, it may have tangible observational consequences for future surveys of terrestrial exoplanets.Provided that all other factors are held equal, a higher atmospheric ion escape rate could translate to relatively thinner atmospheres.Hence, if future surveys either discover thinner atmospheres -or, even better, obtain direct evidence of higher nonthermal ion escape rates -at a certain value of R (conceivably R ∼ 0.7 R ⊕ ), such data would serve to corroborate the predictions of our numerical simulations.Moreover, such nonmonotonic behavior would deepen our understanding of how nonthermal atmospheric escape can sculpt the properties of exoplanetary atmospheres of rocky worlds.

Figure 1 .
Figure1.The calculated atmospheric ion escape rates as a function of planetary radius, which manifest a nonmonotonic trend, peaking at R = 0.75RV .

Figure 2 .
Figure 2. Logarithmic scale contour plots of the escaping O + ion density (units of cm −3 ) in the meridional plane for cases with different planetary radii subjected to the same ancient solar conditions.The second plot from the left in the top row (R = 0.75 RV ) has the highest atmospheric ion escape rate.

Table 1 .
The lower bound is approximately the size of Mars, which was capable of retaining an atmosphere at early epochs

Table 1 .
Model exoplanet parameters for rocky exoplanets.

Table 2 .
The calculated atmospheric ion escape rates (of various ion species) in units of 10 26 s −1 as a function of the planetary radius.