The Impact of Cometary "Impacts" on the Chemistry, Climate, and Spectra of Hot Jupiter Atmospheres

For the sake of computational efficiency and broader model comparability, our cometary ablation and breakup model is highly parameterized. We assume that the comet encounters the atmosphere at a zero angle of incidence (i.e., $\cos (\theta )=1$) and remains undeformed (i.e., spherical) until breakup. Further, while the impact itself is modeled in a time-dependent manner, tracing the journey of the comet through the atmosphere, the resulting changes to the atmospheric chemistry via mass deposition are all applied to the atmospheric model concurrently.

The impact itself can be split into two phases: in the upper atmosphere, where the pressure is low, the comet slows due to atmospheric drag, drag whose friction drives surface ablation. However, as the density increases, so too does the drag, leading to an increase in the ram pressure and eventual breakup when the ram pressure exceeds the tensile strength of the comet (Passey & Melosh 1980; Mordasini et al. 2016). During the ablation phase, at each time step, the velocity of the impactor is calculated using the velocity evolution equation of Passey & Melosh (1980; and references therein):

$$\begin{eqnarray}&&\displaystyle \frac{{dV}}{{dt}}=g-\displaystyle \frac{{C}_{D}{\rho }_{a}{{AV}}^{2}}{M},\end{eqnarray} \tag{ 1 }$$

where

$$\begin{eqnarray}&&A={S}_{F}{\left(\displaystyle \frac{M}{{\rho }_{c}}\right)}^{\tfrac{2}{3}}\end{eqnarray} \tag{ 2 }$$

is the effective cross-sectional area of the comet, g is the gravitational acceleration of the planet, C_D is the drag coefficient, ρ_a is the atmospheric density, S_F is the shape factor of the comet (S_F = 1.3 for a sphere), and ρ_c , V, and M are the cometary density, velocity, and mass, respectively.

Similarly, the ablation-driven mass deposition is given by

$$\begin{eqnarray}&&\displaystyle \frac{{dM}}{{dt}}=-\displaystyle \frac{{C}_{H}{\rho }_{a}{{AV}}^{2}}{2Q}\left(\displaystyle \frac{{V}^{2}-{V}_{\mathrm{cr}}^{2}}{{V}^{2}}\right),\end{eqnarray} \tag{ 3 }$$

where C_H is the heat transfer coefficient, Q is the heat of ablation of the cometary material (see Table 1), and V_cr = 3 km s⁻¹ is the critical velocity below which no ablation occurs due to the reduced deposition of frictional energy (Passey & Melosh 1980).

At the same time, the ram pressure (P_ram) is also calculated and checked against the tensile strength (σ_T ) of the comet, with breakup occurring when

$$\begin{eqnarray}&&{P}_{\mathrm{ram}}\gt {\sigma }_{T},\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&{P}_{\mathrm{ram}}={C}_{D}{\rho }_{A}{V}^{2}.\end{eqnarray} \tag{ 5 }$$

At this point, the impact is terminated and any remaining mass is distributed deeper into the atmosphere over a pressure scale height using an exponentially decaying function in order to simulate the rapid breakup (i.e., shredding/crushing) of the comet and the resulting mass distribution due to the inertia of the incoming material (Korycansky et al. 2006).

While relatively simple, when properly configured (see Section 2.3), this approach reproduces both the observed breakup location of comet Shoemaker–Levy 9 with Jupiter in 1994 (around 1 bar; see, for example, Korycansky et al. 2006), as well as breakup locations calculated using more complex physical processes that are believed to underlie cometary destruction. This includes cometary deformation (i.e., pancaking; e.g., Mac Low & Zahnle 1994; Field & Ferrara 1995; Mordasini et al. 2016) as well as instabilities such as the Rayleigh–Taylor (Field & Ferrara 1995; Alibert et al. 2005c) or Kelvin–Helmholtz (Korycansky et al. 2002) instabilities. The latter of these instabilities drive breakup by inducing surface waves that lead to the shedding of droplets of surface material, and they eventually grow to be on the order of the size of the comet, disrupting it and leading to breakup. Next, once the mass-distribution profile for the impact has been calculated, it is then converted into a mass-density profile by evenly distributing the mass at each vertical level over a 150 km sided box. This box size was chosen to ensure that comets of all radii of interest can be modeled: bigger box sizes lead to smaller comets having a negligible impact on the atmospheric composition, whereas smaller box sizes lead to instabilities in the chemical solver, as the most-massive comets completely overwhelm the unperturbed initial atmosphere. This occurs when any constituent component of the atmosphere becomes more fractionally abundant than hydrogen, even temporarily (e.g., during initial material deposition); this can easily be avoided by ensuring that the incoming mass is distributed/mixed over a sufficiently large area. We note that, due to the 1D nature of our atmospheric models, changing the box size is degenerate with changing the comet's mass (increasing the box size means distributing the mass over a larger area, resulting in local density changes that are equivalent to a smaller comet spread over the original, smaller, area). Hence, we chose to keep the box size fixed for our models, varying the comet radius/density as a more intuitive mechanism by which to adjust the amount of cometary material delivered to the atmosphere. Finally, the cometary composition (Table 2) is used to convert this mass density into a number density for each of the constituent species; this is combined with our unperturbed hot Jupiter model to generate an initial molecular abundance profile that will be evolved by ATMO.

The impact of a comet/planetesimal might also lead to atmospheric heating due to at least some of the kinetic energy lost by the impacting comet being transferred to the surrounding atmosphere, with the exact ratio of kinetic energy lost to ablation and the surrounding atmosphere being controlled by the heat transfer coefficient C_H (Alibert et al. 2005a; Ragossnig et al. 2018). For our models, we set C_H = 0.5 (which is the upper limit of the range of potential values given by Svetsov et al. (1995) and Alibert et al. (2005a)), evenly splitting the lost kinetic energy between the comet (and hence thermal ablation) and atmospheric heating (which takes the form $0.5\,{C}_{H}\left({m}_{t}{v}_{t}^{2}-{m}_{t+1}{v}_{t+1}^{2}\right)\,=\tfrac{3}{2}{k}_{B}T$). Here, we found that the long-term effect of such a heating term was minimal. While it does lead to some initial heating of the upper atmosphere, the short radiative timescales of this region, paired with the strong dayside irradiation caused by the very close host star, mean that, after a few days/weeks of model time, models with and without kinetic energy driven heating are nearly indistinguishable. This remains true on the dark nightside, when redistribution is efficient and hence strong day–night heating can play a similar role (Showman et al. 2020). However, for the sake of energy conservation, we retain this initial heating term for all models shown here.

ATMO is a 1D "radiative-convective" disequilibrium chemistry model for simulating the structure and composition of hydrogen-dominated planetary atmospheres, such as hot Jupiters and brown dwarfs (Amundsen et al. 2014; Tremblin et al. 2015, 2016; Drummond et al. 2016). It has been regularly applied to interpret transit observations, both as a forward and retrieval (i.e., recovering observations) model, including the analysis of early-release JWST observations (e.g., Fu et al. 2022). In fully consistent modeling mode, ATMO iteratively solves for both the atmospheric pressure–temperature profile and chemical abundances while including radiative transport (Amundsen et al. 2014, 2016) and chemical kinetics (Drummond et al. 2016), including vertical mixing and photochemistry. The inclusion of chemical kinetics is key for a number of reasons, ranging from the strong vertical mixing associated with molecular abundance gradients that result from the impact, to the significant fraction of a comet's mass that is made up of water, a significant opacity source (Seager & Sasselov 2000; Fortney et al. 2008; Burrows et al. 2010) that is susceptible to photodissociation in the presence of strong UV radiation (Fleury et al. 2019). Further, a number of recent studies, such as Moses (2014) and Molaverdikhani et al. (2020) have shown that the inclusion of chemical kinetics, turbulent and molecular diffusion, photochemical processes, and free election interactions can all significantly influence the resulting steady-state structure.

For this study, we consider the archetypal hot Jupiter HD209458b as our test bed and initial unperturbed atmosphere. HD209458b has a radius of 1.35 R_J, surface gravity $\mathrm{log}(g)=2.98$, and orbital radius of 0.048 au (Charbonneau et al. 2000; Southworth 2010). The host star is Sun-like (a G0 star versus the G5 Sun) with R = 1.114 R_⊙, hence the short-wavelength stellar stellar spectrum is a modification of the Sun-like spectrum of Kurucz (1970) ¹ (as well as Kurucz (1991) and Castelli & Kurucz (2003)), while the XUV spectrum (used for 34 reactions included in the photochemistry scheme) is purely Sun-like because HD209458 exhibits Sun-like activity levels (Louden et al. 2017). This irradiation intercepts our 1D atmospheric model at an angle of 45°, chosen because such a configuration typically provides a good fit to observations (Irwin et al. 2020).

Vertically, our models extend between 10⁻⁵ and 10² bar, allowing us to model both the upper atmosphere, which is probed by transmission spectra, as well as the deep, optically thick atmosphere in which the cometary breakup and resulting inertia-driven mass deposition occurs. These regions are linked via vertical mixing, which is parameterized using a constant eddy diffusion coefficient K_zz = 10⁹ cm² s⁻¹, a value that was chosen to be broadly compatible with prior studies using ATMO (Drummond et al. 2016; Venot et al. 2020) and of hot Jupiter atmospheres (Moses et al. 2011; Miguel & Kaltenegger 2014; Molaverdikhani et al. 2020; Miles et al. 2023; Samra et al. 2023). We note that small changes in the adopted value of the mixing strength (K_zz = 10⁸ → 10¹⁰) did not significantly alter our results, other than the time taken for the model to settle.

The radiative transfer scheme takes a correlated-k approach (Goody et al. 1989; Amundsen et al. 2014) to computing the total opacity, with the model including 20 opacity sources: H₂–H₂ collision-induced absorption or CIA, H₂–He CIA, H₂O, CO, CO₂, CH₄, NH₃, Na, K, L, Rb, Cs, FeH, PH₃, H₂S, HCN, C₂H₂, SO2, Fe, and H⁻ free–free and bound–free; see Table 1 of Phillips et al. (2020) for details regarding these opacity sources. Under this approach, the incoming radiation is split in a number of bands in which opacities are more easily parameterized, 32 for our actual atmospheric models and 500 when calculating transmission-spectra/optical depths, uniformly distributed in wavenumber space between 31 cm⁻¹ and 50,000 cm⁻¹.

Finally, the chemical-kinetics scheme used is the C₀–C₂ chemical network of Venot et al. (2012), which includes 957 reversible and 6 irreversible reactions involving 105 species, plus helium as a potential third body in some reactions.

Cometary parameters (Table 1) have been selected to be broadly compatible with observations of comets within our own Solar System, as well as recent observations of interstellar ices, presuming that said ices are preserved during the star and planet formation process. Where a choice existed in the exact value of a parameter, such as tensile strength, either an aggregate of reported values was selected, or we ran a series of tests in order to find the value that best reproduced the breakup of comet Shoemaker–Levy 9 (SL9).

To that end, we consider:

• 1.  
  A cometary radius range of between 2.5 km and 30 km (see the review by A'Hearn 2011, including references therein, as well as Schambeau et al. 2015).
• 2.  
  Two cometary densities, 1 g cm⁻³ (slightly compacted or "dirty" H₂O ice) and 2 g cm⁻³ (a denser compacted ice mix), both of which are broadly compatible with models/observations (see the review of Weissman & Lowry 2008) under the assumption that the atmospheric pressure has resulted in compressive waves that enhance sub 1 g cm⁻³ densities to ∼1 g cm⁻³ (Field & Ferrara 1995).
• 3.  
  An initial cometary velocity of 20 km s⁻¹, in agreement with the majority of the models of Korycansky et al. (2002).
• 4.  
  A drag coefficient equal to that of a sphere in a moderately turbulent flow, C_D = 0.5.
• 5.  
  A latent heat of ablation equivalent to that of pure water, i.e., Q = 2.5 × 10¹⁰ erg g⁻¹ (Field & Ferrara 1995; Mordasini et al. 2016).

As for the tensile strength, estimates for comets vary widely, with studies of both different bodies and different locations on the same body reporting significantly different values. For example, Spohn et al. (2015; as well as Basilevsky et al. 2016 and Möhlmann et al. 2018) report values for comet 67P/Churyumov–Gerasimenko, based upon Rosetta mission observations, ranging from 6.5 × 10⁴ erg cm⁻² at the Philae lander impact site to 4 × 10⁷ erg cm⁻² within the Abydos valley region of the comet. As such, during initial model development, we tested cometary impact models with tensile strengths between that of loosely packed snow and solid water ice, 10³ and 10⁷ erg cm⁻² respectively (Petrovic 2003), exploring which value best reproduced the breakup pressure of comet SL9 (0.1 → 1 bar). We found that a tensile strength of between 10⁶ and 10⁷ erg cm⁻² was optimal, and as such, we settled on σ_T = 4 × 10⁶ erg cm⁻², as used by Mordasini et al. (2016).

Finally, we come to the cometary composition (Table 2), which we treat a bit differently. Rather than converging on a single composition based upon a broad mean of literature values, we instead investigate nine different composition models. This allows us to study how differences in the composition of the incoming material are reflected in the final atmospheric temperature and chemistry. The compositions we consider include:

• 1.  
  A model of the solar nebula that lead to the formation of our Solar System (Lodders 2003, Table 11).
• 2.  
  A Solar System cometary composition model taken from Bockelée-Morvan et al. (2004; referred to as comet ver. A).
• 3.  
  A model of the composition of comet 67P/Churyumov–Gerasimenko based upon Rosetta mission data taken from Le Roy et al. (2015;  referred to as comet ver. B).
• 4.  
  A modification of the aforementioned 67P cometary composition in which the CO₂ abundance has been significantly enhanced (it should be noted that we assume all other cometary parameters, such as Q, are unchanged).
• 5.  
  A model of the composition of the protoplanetary disk CRBR 2422.8-3423, based upon Spitzer observations, from Pontoppidan et al. (2005).
• 6.  
  A pair of planetesimal models, based upon the observationally driven planet formation models of Alibert et al. (2005b), for which two different N₂/NH₃ ratios where considered (N₂/NH₃=1 in ver. A and N₂/NH₃=10 in ver. B).
• 7.  
  A pair of interstellar ice models based upon JWST observations of two sightlines (ver. A and ver. B) probing the very outer envelope of the class 0 protostar Cha MMS1, taken from McClure et al. (2023).

We note that, in general, our composition models are carbon depleted relative to the Sun (C/O = 0.55; Asplund et al. 2009), suggesting that the observed objects/regions either have a lower C/O ratio than our Solar System, formed interior to the snowlines of key carbon bearing species (which significantly affects cometary C/O ratios; Öberg et al. 2011), or contain additional refractory carbon that is missing from the listed volatile inventories.

We start by exploring the initial mass-deposition profiles (Figure 1) that are common to all of our cometary models regardless of composition. Here, the two distinct regimes that make up the cometary impact are clear: in the upper atmosphere, where the atmospheric density/pressure is low, the incoming comet retains its global integrity, with only surface material being ablated and deposited into the atmosphere at a rate that depends upon the size of the comet (Equation (3))—larger comets exhibit higher ablation rates. However, as the comet continues its downward journey, the atmospheric density/pressure increases, thus increasing both the rate of ablation as well as the stresses on the comet. Eventually, the density becomes high enough that the ram pressure exceeds the tensile strength of the comet (Equation (5)), leading to breakup and the final distribution of any remaining cometary material. As can be seen in Figure 1, the exact location at which the breakup occurs is dependent upon the cometary mass, with larger comets experiencing weaker breaking due to atmospheric drag (Equation (2)) and hence moving faster, increasing the ram pressure for an otherwise identical atmospheric density.

Zoom In Zoom Out Reset image size

Examples of the initial fractional molecular abundance (i.e., volume-mixing ratio) profiles that result from processing the mass-deposition profile with a cometary composition model (Section 2) and combining it with a unperturbed HD209458b-like atmospheric model can be found on the top row of Figure 2, with the left profile resulting from an R = 5 km, ρ = 1 g cm⁻³, cometary impact (one of the smaller comets considered here) and the right from an R = 22.5 km, ρ = 1 g cm⁻³, cometary impact (one of the larger comets).

Zoom In Zoom Out Reset image size

Here, we can see a number of features that will influence the steady-state atmospheric composition. For example, because the abundances of CO₂, CH₄, and NH₃ in the unperturbed atmosphere are low, even the smallest amount of deposited cometary material will lead to a large increase in the initial abundance, with the relative increase being much larger than that found for H₂O or CO. However, this is not to say that the initial enhancement in H₂O and CO abundances is small: The R = 22.5 km comet impact model reveals that, at around the breakup pressure, H₂O has become the second-most abundant component of the atmosphere, and CO is more abundant than He. However, all of this is temporary: vertical mixing leads to this strong local enhancement rapidly becoming more uniformly spread throughout the atmosphere, with H₂O/CO becoming less abundant than He at all pressures in around three months of simulation time.

Moving on to the steady-state fractional abundances, which are shown on the bottom row of Figure 2, we find that the impact of the smaller comet (left) has had a significant effect on the atmospheric chemistry. For example, the water content of the atmosphere is over an order of magnitude higher than that found in the unperturbed atmosphere, with this effect growing to a factor of ∼1000 as the impactor size increases (right).

Starting with the smaller R = 5 km cometary impact, which reaches steady state after approximately 25 yr of simulation time, we find that H₂O, CO, and CO₂ have been uniformly enhanced, with H₂O overtaking CO in overall abundance and CO₂ showing the strongest relative enhancement: 80 times versus four and two times for H₂O/CO, respectively. This can be linked to the change in C/O ratio, and more generally to the overall oxygen content, of the atmosphere—for example, the unperturbed HD209458b atmosphere's C/O ratio was 0.58 (i.e., close to solar), whereas this atmosphere's C/O ratio is 0.48. Consequently, we have seen a shift in the carriers of carbon away from hydrocarbons, such as CH₄ or C₂H₂ (reduced by, on average, respective factors of 10⁴ and 10⁹) and toward oxygen-rich atmospheric components, such as CO₂. This occurs despite the fact that additional CH₄ was introduced to the atmosphere by the cometary impact. Looking into CH₄ in more detail, we find that the largest reduction in abundance (5 orders of magnitude) occurs in the upper atmosphere, at pressures less than 10⁻² bar, suggesting that the primary driving force behind this change is photolysis, specifically the photolysis of water. During and after the cometary impact, ablation and vertical mixing result in a high abundance of water in the very upper atmosphere. Here, the strong stellar irradiation causes water photodissociation, freeing hydrogen, oxygen, and hydroxide (OH) that go on to form CO, CO₂, O₂, NO, NO₂, etc. (see the extended Section 3.5). This explanation is reinforced by the abundance of excited oxygen, O(¹D), which is massively enhanced at early times (average abundance of 1.3 × 10⁻⁸ versus an unperturbed abundance of 4.4 × 10⁻¹³, for our 22.5 km impact case) when H₂O photolysis is at its peak, and remains enhanced by around 2 orders of magnitude at steady state (4.02 × 10⁻¹¹). It may also be responsible for the factor ∼10 enrichment of atomic hydrogen (H) found in the middle to deep atmosphere (although the temperature of the deep atmosphere may also play a role; see Section 3.3). It should be noted that we have confirmed that the primary driving forces behind this CH₄ destruction are reactions with H and OH, both of which are created by water photolysis in the upper atmosphere.

Moving on to the R = 22.5 km cometary impact, which reaches a steady state after approximately 45 yr of simulation time, we find an atmosphere that is highly influenced by the low C/O ratio and hence high oxygen content of the incoming material. Once again, we find that H₂O, CO, and CO₂ have been uniformly enhanced, by 14, 22, and ∼1000 times, respectively. However, unlike in the R = 5 km cometary impact, we no longer find that H₂O is more abundant than CO. Despite the fact that more H₂O is delivered to the atmosphere than CO by the comet, the final CO/H₂O ratio tilts even more in favor of CO than in the unperturbed model (see Section 3.4). However, just because the oxygen is not in water does not mean that it is not present, as the C/O ratio will attest: here we find C/O = 0.42, suggesting that the carbon-to-oxygen ratio of the atmosphere is converging toward that of the incoming comet (this can be seen in Section 3.6). The oxygen-rich nature of the atmosphere is also reflected in the abundances of hydrocarbons, such as CH₄, whose global fractional abundance has dropped by a factor of ∼10⁴ and very upper atmosphere abundance has dropped by more than 10 orders of magnitude. We once again attribute this effect to water photolysis, which, as we discuss in Section 3.4, has grown nonlinearly with incoming cometary mass.

In the above analysis, we have not discussed in detail how nitrogen-bearing molecules behave. This is because, unlike the molecules discussed above, we find neither a uniform increase or decrease in abundance for NH₃, NO, HCN, or even oxygen-rich NO₂. Instead, we find that the change in abundance is pressure dependent. For example, at lower impactor masses, we find that NH₃ is enhanced in the upper and lower atmosphere while being diminished in the middle atmosphere (between ∼10⁻⁴ and ∼10⁻² bar). Then, as the impactor mass increases, we find a general increase in NH₃ except in the very upper atmosphere (P < ∼ 10⁻⁴ bar) where the trend is reversed and NH₃ is further diminished. Yet, despite this increase, for all the models considered here, the abundance of NH₃ in the middle atmosphere never exceeds that found in our unperturbed HD209458b model.

As for why this trend occurs, the general enhancement in NH₃, and other nitrogen-bearing molecules, is due to a combination of the influx of nitrogen from the impacting comet as well as the increase in atmospheric metallicity favoring the formation of heavier molecules (see Section 3.6 for more details). At the same time, the upper atmosphere suppression can be linked to photolysis, both directly, via the breakdown of NH₃, and indirectly, via the presence of additional "free" oxygen linked to water photolysis. Together, these effects also lead to the formation of more oxygen-rich nitrogen-bearing molecules, such as NO and NO₂, both of which are enhanced by over 5 orders of magnitude, and both of which are heavy and hence prone to sinking out of the upper atmosphere.

Our cometary impacts affect more than just the atmospheric chemistry. As shown in Figure 3, the temperature of the atmosphere has increased at all pressures, with the strongest heating occurring in the middle/deep atmosphere (P>10⁻¹ bar). Here, we find that, except for the largest impactors, the temperature of the deep atmosphere increases with impactor size (we discuss why this relationship breaks down for large impactors in Section 3.4). This change is driven by the presence of water: essentially, water acts as a greenhouse gas, reducing the amount of heat radiated away from the deep atmosphere and hence increasing the temperature. As for where the energy to heat the deep atmosphere comes from, ATMO includes an artificial internal energy source that is designed to heat the deep atmosphere in order to account for the radius inflation of hot Jupiters, something that 1D "radiative-convective" convective models cannot account for a priori (this is the so-called radius-inflation problem). However, the artificial nature of this heating in ATMO does not preclude such an effect from occurring in a real hot Jupiter that experiences a similar H₂O enrichment. As discussed in Tremblin et al. (2017) and Sainsbury-Martinez et al. (2019, 2021, 2023), the most likely solution to the radius-inflation problem is the vertical transport of heat from the upper to the deep atmosphere, where it heats the internal adiabat (and also causes the internal adiabat to form at lower pressures, increasing the apparent radius). This heating would be similarly affected by an enhanced water greenhouse effect, because the input energy source does not rely on radiative transport, and deep cooling—radiative loss—would be similarly suppressed. As such, given the link between the temperature of the deep atmosphere and the impactor mass, and hence atmospheric C/O ratio (and metallicity), we propose that said C/O ratio might act as a secondary effect controlling the exact level of radius inflation observed for a hot Jupiter.

Zoom In Zoom Out Reset image size

We note that the change in temperature of the lower/deep atmosphere has consequences for the amount of free/atomic hydrogen in the deep atmosphere, with heating causing increasing amounts of molecular hydrogen to dissociate into atomic hydrogen.

As discussed above and shown in Figure 4, we find that the abundance of H₂O departs from a linear cometary-mass/H₂O-abundance relationship at very high impactor masses. More specifically, we find that there is a threshold for water deposition above which the steady-state water abundance is either constant or lower than that found with a lower impactor mass. It should be noted that the exact location of this threshold is composition dependent (i.e., cometary water fraction dependent): for example, we find that the water enhancement starts to be suppressed when R > 20 km for a 1 g cm⁻³ comet with JWST Ice ver. A composition.

Zoom In Zoom Out Reset image size

Our analysis suggests that the cause of this H₂O deficiency at high impactor mass is a combination of vertical mixing and photolysis. Massive impacts cause a correspondingly massive initial influx of H₂O, primarily in the lower atmosphere, which is then rapidly mixed into the upper atmosphere due to the resulting density/molecular gradient, where it then undergoes photolysis. While this mixing/photolysis occurs for all impacts, as the impactor mass grows, so too does the amount of water that ends up in the upper atmosphere at early times. This water is then exposed to significant UV irradiation, leading to strong initial photolysis and hence an overall reduction in the total total atmospheric H₂O by the time that the atmosphere has reached a steady state. This conclusion is reinforced by both our alternate (higher/lower) K_zz tests, which reveal that more efficient diffusion of water out of the upper atmosphere mutes this enhanced photolysis, and our multi-impact tests (in which cometary material is spread over 10 smaller impacts), which reveal slightly slower photolysis and hence higher steady-state water abundances (due to the slower delivery of water to the atmosphere).

Interestingly, we do not find such a break in the cometary-mass/abundance relation for other molecules. However, theoretically, it should occur for any other molecule that undergoes photolysis, when its initial abundance is large enough. For example, hints of this trend are found for our high-CO₂ cometary composition models, where we find that the rate of CO₂ abundance increase with impactor mass appears to slow at large impactor sizes. However, because of the numerous other paths by which CO₂ can become depleted in these models, the effect is rather muted.

We next explore how the abundances of other molecules change with impactor size (mass/radius). Here, we find that our models typically fall into one of a number of regimes, examples of which are shown in Figure 5. We note that the trends below are discussed in relation to an R = 2.5 km, ρ = 1 g cm⁻³, cometary impact model, not the unperturbed HD209458b-like atmosphere. Consequently, a trend may be referred to as an enhancement with impactor mass/radius even when the molecular abundance is less than that found in the unperturbed model—it is just less by a smaller amount.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

View figure set (42 panels)

Tarball of figure set images

The first regime is a simple, uniform, increase in abundance with impactor mass, as seen for H₂O below the "saturation" limit (Figure 4). Other molecules that behave like this include other oxygen-bearing molecules, such as CO, CO₂, and O₂.

Next, we have molecules such as NO (Figure 5(D)), O(³P), H, and OH, which follow a similar trend of increasing abundance with impactor mass, but with a pressure dependence in the strength of said enhancement: generally, we find that the relative enhancement increases with pressure. We attribute this pressure dependence to both the effects of photochemistry, which is confined to the upper atmosphere, as well as molecular mass.

We also find a large number of molecules that are increasingly diminished in abundance with impactor mass in the upper atmosphere, but that show much weaker impactor-mass dependence in the lower atmosphere. An example of this trend can be seen for CH₄ (Figure 5(A)), which, as we discuss in Section 3.2, shows a relatively weak decrease in abundance with impactor mass in the deep atmosphere, and a much stronger dependence of abundance on impactor mass in the upper atmosphere, where free oxygen and hydroxide are abundant, leading to oxygen-rich molecules being chemically preferred. A similar story holds true for other carbon/nitrogen-rich and oxygen-poor molecules, such as NH₃ (Figure 5(C)—although here the deep enhancement is stronger due to overall nitrogen enrichment by the comet), CH, CH₃, and C₂H₄.

Finally, there are a number of molecules for which no clear trend is evident. For example, HCN (Figure 5(B)) exhibits a rather interesting impactor-size/abundance dependence in the upper atmosphere: when the impactor is relatively small, we find that the abundance of HCN decreases with increasing impactor size; however, when the impactor size is relatively large, the trend is reversed, and the abundance instead increases slightly with impactor size. The underlying cause of this relationship is likely to be the nonlinear interaction of a number of competing effects, including the overall enrichment in atmospheric nitrogen caused by larger impactor, a more general increase in atmospheric metallicity with increasing cometary-mass deposition, and photochemistry in the upper atmosphere changing the local atmospheric composition. For HCN, in the lower atmosphere, the primary controller of the abundance is the total nitrogen content of the atmosphere, hence we find little to no enhancement in HCN for any cometary impacts that do not carry nitrogen-bearing molecules, such as the Comet ver. A composition. Furthermore, even for those impacts that do introduce additional nitrogen to the atmosphere, the upper atmosphere HCN abundance still remains less than that found in our unperturbed model.

Examples of additional molecules, which may be abundant enough in the upper atmosphere to impact observations, are shown in extended Figure 5. Together, these profiles reveal the complex ways in which molecular abundances change as we adjust the composition of the atmosphere, reinforcing the need for robust models when interpreting observations of planetary atmospheres.

In order to study the full ensemble of more than 400 cometary impact models considered here, we next calculated the "global" atmospheric C/O ratio for each model and the upper atmosphere (P < 10⁻² bar) mean fractional abundance for each molecule in each model. We focus on the upper atmosphere because it is this region that is potentially observable with transmission spectra (i.e., it is here that the opacity typically reaches unity—see Section 3.7). We also limit ourselves to molecules with upper atmosphere abundances >10⁻¹⁶, because lower-abundance molecules are not only highly sensitive to slight perturbations in the chemistry (driven by small number statistics) but are also highly unlikely to be detectable. Abundance–C/O maps and trends for all of the molecules with upper atmosphere abundances that fulfill this criterion are shown in the extended version of Figure 6 and given in Table 3, respectively.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

View figure set (46 panels)

Tarball of figure set images

It is important to note that our ensemble of models only samples a limited region of the C/O parameter space. For the nine different cometary compositions under consideration, there are eight different cometary ice C/O ratios, because the planetesimal models differ only in nitrogen content. Furthermore, as the size of the impactor increases, the relative effect that it has on the atmospheric composition/density also grows, leading to the overall atmospheric C/O ratio slowly converging toward that of the incoming material. This can be seen in Figure 6, where we find eight distinct abundance–C/O ratio curves. In a way, the steps along each curve are indicative of how ongoing bombardment by material with very similar compositions (or C/O ratios) might change the global chemistry of a hot Jupiter.

Here, as shown in Figure 6, we focus on four molecules that are exemplary of most features found for our ensemble of cometary impact models. Generally, these features are either an increase or decrease in abundance with increasing atmospheric C/O ratio and/or impactor mass (and hence atmospheric metallicity). However, there are a number of molecules whose behavior falls outside these bounds. For example, the upper atmosphere abundance of CO₂ is essentially independent of the C/O ratio, only being affected by the increase in metallicity associated with cometary impacts. We also find that, as discussed in Section 3.2 and Section 3.5, the behavior of nitrogen-bearing species follows its own trends.

We now look at these trends in more detail, starting with oxygen-rich molecules. As shown in Figure 6(A) for H₂O, we find that the abundance of oxygen-rich molecules increases with decreasing atmospheric C/O ratio (and increasing metallicity). This occurs because, at lower C/O ratios, the amount of free oxygen available rises, leading to reactions that form oxygen-rich molecules being preferred over the formation/presence of hydrocarbons. At the same time, regardless of cometary composition, we find an increase in the abundance of most oxygen-bearing molecules with impactor size. This effect is especially strong for water because, apart from one theoretical composition, water is the primary constituent of our cometary ices and hence is directly delivered, in large quantities, during an impact. In brief summary, this occurs because cometary impacts change not only the C/O ratio of the atmosphere, but also the metallicity, especially the C/H, O/H and N/H number ratios. This is because these ratios are orders of magnitude higher in cometary ice than in an unpolluted, H₂-dominated, hot Jupiter atmosphere.

Metallicity-driven effects also explain the trend found for CO₂, as shown in Figure 6(C). Here, we find little to no change in upper atmosphere abundance with C/O ratio, but strong abundance dependence on the impactor mass, and hence the overall carbon/oxygen content of the atmosphere.

Next, we have a number of molecules, typically hydrocarbons such as CH₃, as shown in Figure 6(B), that increase in abundance with C/O ratio but decrease in abundance with impactor size. The former effect is simply the result of additional atmospheric carbon leading to more hydrocarbons forming, while the latter effect occurs because, regardless of composition, the cometary C/O ratio is always <1, hence more oxygen is always deposited than carbon.

Finally, we find that the upper atmosphere abundances of most nitrogen-bearing molecules such as NO, as shown in Figure 6(D), behave rather differently from the above trends. At low impactor masses, we find either a simple increase or decrease in abundance with impactor size. However, as this size grows, we eventually cross a threshold and the trend reverses. In order to understand why this occurs, it is important to note that, as seen in Figure 5(C/D), such reversals are limited to the upper atmosphere, vanishing when we calculate the global mean. Thus, the local nature of this reversal, combined with the universality of this trend (i.e., it occurs even when comets do not introduce nitrogen to the atmosphere), suggests that the underlying cause is photochemistry. Not only is the upper atmosphere very highly UV irradiated, but all of the carriers of nitrogen in our model are either photochemically active or directly form from molecules that result from the photodissociation of other abundant molecules, such as water. For example, at early times, the strong photolysis of incoming water leads to a strong initial switch from NH₃ toward NO and NO₂ in the upper atmosphere. This is then vertically mixed, resulting in the lower-atmosphere enhancements seen in Figure 5(C/D). However, as the amount of water in the upper atmosphere grows, so too does the strength of this initial photolysis (see Section 3.4), leading to increasingly strong initial NH₃ depletion and NO/NO₂ enhancement. In turn, the compositional gradients for these molecules between the upper and lower atmospheres also grows, strengthening vertical mixing and hence further enhancing the deep NO and NO₂ enrichment. The relatively high mass of NO and NO₂ molecules then bakes this deep enrichment in. Because both NO and NO₂ are heavy, they tend to stay in (or sink into) the lower atmosphere, even when the compositional gradients that drove the initial mixing/diffusion have disappeared. On the other hand, NH₃ is much lighter and hence is more freely mixed between the upper and lower atmospheres as the simulation progresses. Thus, the switch in NH₃ abundance in the upper atmosphere can also be linked to the nonlinear change in water abundance with impactor mass. For the most massive impactors, the water content of the atmosphere is actually reduced, hence more NH₃, which was sequestrated in the middle/lower atmosphere at earlier times, can be mixed back into the upper atmosphere and survive, leading to the reversal/enhancement seen in Figure 5(C). Such a scenario also explains the difference in reversal threshold impactor masses for different compositions (Figure 6(D)). As the cometary composition changes, so too does the water (oxygen) and nitrogen content of the upper atmosphere, impacting reaction (production) and photolysis rates and hence the point at which we switch between enhancement and depletion. This conclusion is reinforced by our analysis (not shown) of the change in abundance of nitrogen-bearing molecules with N/O ratio. At "low" N/O ratios, oxygen-bearing nitrogen species are preferred, whereas at "high" N/O ratio, hydronitrogens become increasingly abundant. This is true in both the upper atmospheres as well as when taking a "global" average.

We note that a similar effect, vis-à-vis photolysis and heavier molecules sinking out of the upper atmosphere, also explains the abundance trends of many oxygen-bearing organic molecules, such as CH₃OH, which show depletion/enhancement with impactor mass when we look at the upper or "global" atmospheric averages, respectively.

To understand how the effects of cometary impacts might be detected, it is important to understand how the aforementioned vertical differences in atmospheric chemistry and temperature impact observable features. To that end, we now focus on a single cometary impact, with R = 10 km, ρ = 2 g cm⁻³, and JWST Ice ver. A composition, exploring the time-dependent changes to synthetic opacity maps (optical depth—Figure 7) and transmission spectra (Figure 8) and comparing them to an unperturbed HD209458b atmospheric model. Due to our models' 1D nature, which means that cometary material remains in the model instead of being mixed with the unimpacted surroundings (a process that occurs over tens of years, as seen for SL9), we treat different points in time as being representative of different scenarios.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Early points in time, e.g., 1 yr and 5 yr after impact, are considered to be an optimistic example of the highly localized changes that might occur shortly after impact, before any dissipation/mixing occurs. As such, observations similar to what we discuss for these points in time represent a best-case scenario in which the cometary impact is recent, occurs close to or on the terminator, and alters the terminator region chemistry enough to be observable. This is similar to the short-lived changes observed in the 2.3 μm methane feature observed after the impact of comet SL9 with Jupiter, which Flagg et al. (2016) linked to local (particulate) material delivery and a resulting change in the opacity and albedo of Jupiter's upper atmosphere. Such changes effectively covered the atmosphere at the 45° latitude of Jupiter (Flagg et al. 2016), and we assume a similar scenario occurs for the terminator region of our impacts models. Consequently, our 1 year and 5 year impact models likely overestimate the effect that a single impact has on the terminator region. Yet we still discuss them here, because the effects observed have important implications for the effect of cometary material on planetary atmospheres, reinforcing the need to pair observations with a robust atmospheric model.

On the other hand, we consider the steady state (25 yr after impact for this model) to be representative of an atmosphere that has undergone continuous and/or ongoing bombardment by comets with similar compositions, leading to strong changes in the global composition/ chemistry of the atmosphere. Such a scenario might occur, for example, during and after the migration of a hot Jupiter from its formation location to its current tidally locked orbit.

We start, in Figure 7, by exploring how our cometary impact changes the atmospheres optical depth (τ). Once the model has settled down slightly from the initial perturbation, i.e., after the first few months of simulation time, we find a prolonged period in which the atmospheres opacity is significantly (Figure 7(A)) enhanced with respect to our unperturbed model, with τ = 1 being reached at an average pressure of ∼10⁻³ bar versus ∼10⁻¹ bar in the unperturbed model. The underlying cause of this spike in opacity appears to be the strong H₂O (and to a lesser extent CO and CO₂) enrichment that occurs during the cometary impact. This material acts as a strong opacity source, absorbing incoming radiation and causing local heating of the upper atmosphere, an effect that is also slightly enhanced by the initial heating associated with the cometary impact. In turn, the hotter upper atmosphere strengthens H₂–H₂ and H₂–He collision-induced absorption (CIA), which is highly temperature dependent (Borysow et al. 2001), shaping the increase in optical depth seen in Figure 7(A/B/C) for almost all wavelengths except those associated with K/Na absorption.

Over time, as the initial upwelling of water settles or is photodissociated, the opacity of the upper atmosphere drops, leading to cooling and hence a weakening of CIA. This can be seen in Figure 7(B) as an increase in the pressure at which the atmosphere becomes optically thick (τ = 1).

This process continues as the atmosphere evolves toward a steady state (as shown in Figure 7(C)), with the depth at which τ = 1 reaching an average of ∼10⁻² bar. Yet this is still much shallower than the average τ = 1 pressure, ∼10⁻¹ bar, in an unperturbed HD209458b-like atmosphere, reflective of the enhanced water content (Section 3.2) and upper-atmosphere temperature (Section 3.3) driven by the cometary impact.

However, it should be noted that, as previously alluded to, this enhancement in opacity is not uniform, with different components of the atmosphere affecting the opacity at different wavelengths. For example, between 4 and 5 μm, all of our atmospheric models exhibit an almost order-of-magnitude spike in optical thickness, a spike that is absent for HD209458b. This is caused by CO and CO₂, which have strong absorptions at 4.6 μm and 4.2 → 4.4 μm, respectively, and which are both significantly enhanced (Section 3.2) in abundance when compared with our unperturbed atmosphere.

To explore these differences in more detail, we now turn to synthetic transmission spectra (Figure 8).

One year after impact, we find a transmission spectrum spectrum (Figure 8(A)) that remains significantly impacted by the still-settling cometary material. In particular, we find that, while the apparent planetary radius, R_P /R_S , is enhanced at almost every wavelength, the strength of most (other than the 4 → 5 μm CO/CO₂ or the 16 μm CO₂ features) absorption features is reduced, with respect to our unperturbed HD209458b atmosphere. This combination suggests that the planet itself is not larger/inflated, but that instead we are probing lower-pressure regions of the atmosphere, where the densities of most molecules causing the observed absorption features are reduced. This is exactly what the opacity maps discussed above (Figure 7(A/D)) reveal: the τ = 1 region, which our synthetic transmission spectrum probes, occurs at a much lower pressure for our impacted model than for our unperturbed model. Furthermore, the shape of this change to the transmission spectrum (and opacity maps) is highly reminiscent of the high-temperature H₂–H₂ CIA absorption curve seen in Figure 1 of Borysow et al. (2001). Thus, the absorption features seen in Figure 8(A) are a combination feature of CIA and H₂O opacities, and the enhanced CIA opacity means that fewer water molecules are needed before the atmosphere becomes optically thick, hence reducing the relative strength of the water-absorption feature. This conclusion is reinforced by the strengthening of the CO₂/CO feature, which occurs because the relative enhancement of these molecules is so boosted by the impacting comet (peaking at a factor of 10⁶ for CO₂) that any reduction in probed atmospheric density is more than balanced out by the increase in fractional abundance.

Five years after impact, although the initial upper atmosphere H₂O enhancement has started to significantly settle/photodissociate, and the atmosphere has become optically thinner (Figure 7(B)), the apparent radius of the planet has continued to increase (Figure 8(B)), and so too has the strength of the absorption features, in particular the H₂O, CO and CO₂ features. This occurs because the very outer atmosphere has cooled slightly (as water settles), weakening the very low pressure H₂–H₂ CIA opacity, and thus enabling us to probe higher-pressure/density/temperature regions of the atmosphere, leading to enhancements in both molecular absorption and CIA effects. This is reflected in both the peak-to-trough feature strength, which has a maximum around twice that found one year after impact, as well as the "continuum" level, which tracks the strength of the deeper/hotter H₂–H₂ CIA rather than just being driven by the altitude at which the atmosphere becomes optically thick. The subtlety of this effect reinforces the value of exploring synthetic transmission spectra over just looking at abundances and opacities. It is the best way to understand how vertical differences in composition and temperature impacts observations.

We note that the above increase in apparent radius, and hence the asymmetry in observed planetary radii between eastern and western terminators, is one of the likeliest ways that a very recent, near-terminator impact might be identified. And such a scenario is not unlikely: as discussed in Shoemaker & Wolfe (1982) and Zahnle et al. (2001), the cratering (impact) rate at the apex of motion (i.e., the eastern terminator) can reach 2.6 times the planetary average, while also dropping by a factor of 15 on the trailing edge (i.e., the western terminator).

While the aftermath of a single impact is highly unlikely to be captured, ongoing, or a period of significant, bombardment by comets with similar compositions might lead to observable changes in the global atmospheric composition/chemistry. We explore such a scenario using a steady-state atmospheric model whose synthetic transmission spectra are shown in Figure 8(C). These spectra are typical of those found throughout our ensemble of steady-state atmospheric models, all of which increase the water (and hence oxygen) content of the atmosphere. Similar to what was found at earlier times, the influx of high-opacity material, such as H₂O, CO, and CO₂ has, via atmospheric heating and CIA effects, resulted in a general enhancement in the continuum level with respect to our unperturbed atmosphere, while also weakening many atmospheric features, such as H₂O, K, and Na, due to our observations probing a lower-density region of the atmosphere (Figure 7(C/D)).

We note that the above enhancement of the continuum level is highly correlated with the atmospheric water abundance and hence impactor mass, leading to some of our higher-cometary-mass models showing significant enhancements in the apparent radius of the planet, at least when integrating over the wavelength range considered here. This is important because some observatories, including JWST, observe exclusively in the near/mid-infrared, and as such, an exaggerated radius might be reported (Mordasini et al. 2016). There are two solutions to this problem: The first is to take additional observations with either a UV or optical telescope, both of which are less sensitive to these water- and CIA-driven effects. The second possible solution is to pair these observations with a robust atmospheric model that includes a complete treatment of atmospheric opacity sources. This works because a number of atmospheric features help to break the degeneracy between the observed "continuum" level and the planetary radius.

Examples of these features have been labeled in Figure 8(C) and include:

• 1.  
  enhanced CO₂ absorption at 2.6 μm 2.8 → 3.9 μm, 4.2 → 4.4 μm, and 15.1 → 15.3 μm;
• 2.  
  enhanced CO absorption at 4.6 → 4.7 μm;
• 3.  
  possible weak NH₃ absorption at 9.0 → 9.1 μm;
• 4.  
  possible OH absorption at >10 μm;
• 5.  
  reduction in CH₄ absorption at 1.9 μm, 2.8 → 3.8 μm, and 7.7 μm;
• 6.  
  slightly suppressed H₂O absorption effects, due to probing lower-density regions of the atmosphere, at 1.3 → 1.6 μm, 1.7 → 2.1 μm, and ∼3 μm, etc.

Many of these features fall into the region in which JWST is expected to be most sensitive, which is the fixed-slit mode of NIRSPEC, which has R → 2600 between 0.97 and 5.27 μm (Jakobsen et al. 2022). The mid-infrared CO₂ feature at ∼15.2 μm is also potentially detectable with MIRI.

Further, by comparing these features with potential efractory/metallic elements, such as K (0.77μm) and Na (0.58μm), which should be unaffected by cometary impacts, it is possible that we might be able to identify objects in which the oxygen or carbon metallicity has been enhanced by cometary material delivery. Of course, this would require us to know/measure the abundance of atmospheric K/Na in the first place, but a even a low estimated refractory-to-volatile ratio might be a good indicator that an object is worthy of further study via both observations and atmospheric modeling.

Such features, if identified and attributed to changes in atmospheric C/O ratio and metallicity, provide the most likely pathway to identifying hot Jupiters that have either undergone significant cometary bombardment in their past or are still being bombarded to this day.

It should be noted that we have focused the above discussion on spaced based observations, such as those using JWST, which are are currently the best way to constrain/observe the presence of water on distant planets. However, future ground-based telescopes, such as ELT, might/will be able to break the veil of telluric contamination and constrain the atmospheric C/O ratio and metallicity via high spectral resolution observations paired with robust atmospheric modeling, such as that available with ATMO.

We reiterate that the above discussion of observations represents a best-case scenario, either capturing a near-terminator impact that very recently occurred or observing a hot Jupiter atmosphere that has been significantly polluted by massive or ongoing bombardment. Consequently, rather than being treated as examples of cometary impacts that observers should go out and hunt for, they should be taken as optimistic examples of how material delivery by cometary impacts might reshape hot Jupiter atmospheres and break the link between current atmospheric composition and formation location.

Also, all of the results presented here are for single impacts and 1D atmospheric models. This adds some inherent uncertainty to some of our results and their applicability. To alleviate these concerns, in addition to the models discussed throughout this work, we ran two additional models in which cometary material delivery was spread over ten smaller impacts. For our 5 km equivalent model, we found that the overall atmospheric abundances were very similar, just with a slight enhancement in CH₄ due to its continuous delivery. This effect became more pronounced as we move to the 22.5 km equivalent model, where we find that continuous material delivery has resulted in an approximately sevenfold increase in overall CH₄ and NH₃ abundance, as well as a 15% increase in atmospheric H₂O. The latter change suggests that the enhanced water photolysis discussed in Section 3.4 is indeed likely to be muted in a real hot Jupiter atmosphere, because horizontal mixing will have an effect similar to that of multiple weaker impacts, spreading the water out horizontally instead of temporally. Thus, combining both multiple smaller impacts (the continuous/ongoing bombardment case) with 3D mixing means that the enhanced H2O dissociation seen here is likely to be observed in only the most water-enriched atmosphere. To truly understand these effects, in the future, we must turn to 3D atmospheric models that also include vertical mixing, which may have an effect opposite to that of horizontal/temporal spreading, delivering material from the deep atmosphere to the very upper atmosphere, where it photodissociates.

We note that a similar situation might also impact the reversal in the abundance of nitrogen-bearing molecules in the upper atmosphere, although the fact that said reversal occurs at smaller cometary sizes in our multi-impact model and appears to be linked with molecular masses and vertical mixing means that we are much more confident in the robustness of this feature.