Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences
You have accessResearch articles

Modelling H3+ in planetary atmospheres: effects of vertical gradients on observed quantities


    Since its detection in the aurorae of Jupiter approximately 30 years ago, the H3+ ion has served as an invaluable probe of giant planet upper atmospheres. However, the vast majority of monitoring of planetary H3+ radiation has followed from observations that rely on deriving parameters from column-integrated paths through the emitting layer. Here, we investigate the effects of density and temperature gradients along such paths on the measured H3+ spectrum and its resulting interpretation. In a non-isothermal atmosphere, H3+ column densities retrieved from such observations are found to represent a lower limit, reduced by 20% or more from the true atmospheric value. Global simulations of Uranus' ionosphere reveal that measured H3+ temperature variations are often attributable to well-understood solar zenith angle effects rather than indications of real atmospheric variability. Finally, based on these insights, a preliminary method of deriving vertical temperature structure is demonstrated at Jupiter using model reproductions of electron density and H3+ measurements. The sheer diversity and uncertainty of conditions in planetary atmospheres prohibits this work from providing blanket quantitative correction factors; nonetheless, we illustrate a few simple ways in which the already formidable utility of H3+ observations in understanding planetary atmospheres can be enhanced.

    This article is part of a discussion meeting issue ‘Advances in hydrogen molecular ions: H3+, H5+ and beyond’.

    1. Introduction

    Hydrogen, the most abundant cosmic element, also dominates the composition of giant planets. Consequently, the most prominent ion species in giant planet atmospheres are the stable H+ and H3+ ions [1]. The proton, H+, is not spectroscopically observable, whereas the spectrum of H3+ is exceptionally rich, particularly the ν2 vibration rotation band in the near-infrared [2], a spectral region accessible from Earth's high-altitude observatories. In fact, the first astronomical spectroscopic detection of H3+, enabled by a confluence of evolving theoretical, laboratory and observational advances, was made in Jupiter's auroral region [3]. Further detections at Saturn and Uranus, and continued ground-based monitoring of H3+ at Jupiter, Saturn and Uranus in the subsequent decades, have demonstrated its remarkable effectiveness as a probe of giant planet upper atmospheres (e.g. [411] and references therein).

    At present, the field of comparative aeronomy—that is, the comparative study of planetary upper atmospheres—relies on a sparse sampling of remote diagnostics, especially for the giant planets. Vertical thermal structure, in particular, is difficult to determine remotely, and yet it plays a primary role in identifying the relevant physical processes at work in planetary upper atmospheres. There are only a handful of temperature profiles obtained for Jupiter [12], Uranus and Neptune, primarily from the Voyager spacecraft [1315], whereas Saturn has now been relatively thoroughly sampled by Cassini [1618].

    The observed exospheric temperatures at all of the giant planets are hundreds of Kelvin hotter than predictions based on solar heating alone, emphasizing a rather fundamental lack of understanding in the energy balance in giant planet atmospheres, and highlighting the need for more thorough spatio-temporal thermospheric temperature constraints. Modellers are actively seeking an explanation for this energy discrepancy, which may simply involve redistribution of auroral energy inputs [1921], or perhaps alternative energy sources, such as wave-driven heating from below [22]. In the meantime, measurements of H3+ temperatures are vital for bridging this knowledge gap, as H3+ is thought to be in quasi-LTE with the surrounding neutral atmosphere [1,23,24], and valuable insights have already been provided by H3+ observations to-date. However, derived H3+ temperatures also suffer from a key ambiguity: the vast majority of ground-based observations are column integrations through the entire ionosphere, from top-to-bottom, and therefore measured emissions result from convolution of the vertical structures in both H3+ density and temperature.

    Here, we investigate how giant planet atmospheric models can help supplement interpretation of H3+ spectroscopic observations. These calculations are concentrated on unravelling the column-integrated density and temperature degeneracy behind the observed spectra, and on minimizing—or at least understanding—the effect of thermospheric gradients on analysis of H3+ datasets. First, in §2, we briefly describe the modelling approach as well as the observational constraints used. Next, in §3, we consider a series of increasingly realistic ionospheric models and examine the complications that can arise in interpreting column-averaged observations. Finally, we combine the insights from these sections in order to demonstrate a preliminary method of retrieving altitude profiles of H3+ temperature from nadir viewing geometry observations.

    2. Methods

    (a) Modelling overview

    The majority of H3+ ions in giant planet ionospheres are in photochemical equilibrium (PCE), as the H3+ chemical lifetime is much shorter than the transport time scale at most altitudes, and therefore the ion continuity equation simplifies to equating local production and loss (i.e. Ps = Ls) [25,26]. While transport processes are still relevant—and highly so for H+, especially at high altitude—the dominance of chemical loss at low altitudes justifies the use of one-dimensional (1D) simulations over those regions, which offer the advantages of simplicity and computational freedom over more dynamically comprehensive three-dimensional simulations.

    Ionization fractions at the giant planets are roughly of order 10−6 [27], indicating that ion chemistry and dynamics are largely inconsequential for the underlying neutral atmosphere. However, this generality is less true at auroral latitudes, where there appears to be a correspondence between auroral emission signatures with complex hydrocarbons and stratospheric hazes [28,29], and where H3+ can act as a thermostat, helping to maintain a cooler thermosphere [30] and limit atmospheric escape.

    Two models are adopted in the present work, which is focused on non-auroral latitudes at Jupiter and Uranus. The simulations are conducted in one dimension, owing to the prevalence of PCE for H3+ distributions, and there are separate neutral and plasma modules in order to enable a more computationally efficient exploration of ion chemistry.

    The neutral module is described in detail by Moses & Poppe [31], and is actually a combination of a meteoroid ablation code [32,33] with the Caltech/JPL 1-D KINETICS photochemical model [34,35]. KINETICS solves the coupled mass-continuity equations as a function of pressure, and (for the neutral species) includes molecular and eddy diffusion transport terms. It has proved to be effective and highly adaptable, having been applied to all of the giant planets, and currently treats 70 hydrocarbon and oxygen species that interact via approximately 500 recently updated chemical reactions [31,36]. Input for the meteoroid ablation code follows from revised constraints on interplanetary dust fluxes in the outer Solar System based on in situ spacecraft data [37]. The resulting oxygenated and hydrocarbon mixing ratios are in agreement with a wide range of observational constraints [31]. Therefore, after adjusting the KINETICS simulations for the solar and geometric conditions explored here, the resulting neutral atmospheres serve as an excellent background for exploring realistic ion-neutral photochemistry at the giant planets.

    Plasma densities and temperatures follow from another 1D model called BU1DIM (the Boston University 1-D Ionosphere Model). BU1DIM was originally developed for Saturn [38,39], though has since been applied to Earth [40] and Mars [41]. Its most recent iteration has been generalized for application to any planetary atmosphere, and includes significantly expanded chemistry [25]. BU1DIM describes the time- and altitude-dependent structure of an ionosphere by solving the coupled continuity, momentum and energy equations for all ion species of interest. Jupiter's magnetic field is specified using results from the Juno spacecraft [42]. At Uranus, however, magnetic field measurements are limited to a single flyby [43]. The primary effect of magnetic fields on 1D ionospheric calculations is to constrain the plasma motion (e.g. introducing a sin2 I term into the expression for vertical ion drift velocity, where I is the magnetic dip angle [38,44]). Therefore—partly due to incomplete knowledge of Uranus' magnetic field, and partly due to the predominance of PCE at H3+ altitudes—magnetic field lines at Uranus are considered to be vertical here in order to focus investigations on the effect of vertical thermospheric gradients on derived H3+ parameters. Modelled ion production rates follow from the attenuation of solar Extreme UltraViolet (EUV; 10–121 nm) and soft X-ray photons (combined, the XUV) [45], which are extrapolated to Jupiter and Uranus based on measurements from the Thermosphere Ionosphere Mesosphere Energetics and Dynamics Solar EUV Experiment (TIMED/SEE) [46]. In addition, secondary ionization and thermal electron heating rates are specified using parameterizations derived from coupled electron transport calculations at Saturn [47]. Aside from solar XUV radiation, no other sources of energy input are considered here (e.g. energetic particle precipitation).

    Early theoretical models of giant planet ionospheres predicted electron densities that were up to an order of magnitude too large based on later spacecraft measurements, with Saturn exhibiting the most extreme discrepancy [27,48]. One commonly adopted mechanism for reducing modelled electron densities in order to better match observations was to convert H+ into a molecular ion via the reaction

    Without the introduction of some form of ion-neutral charge-exchange reaction, such as (2.1), modelled H+—and hence electron density, ne—is unrealistically large, as the radiative recombination rate coefficient for H+ is extremely slow (approx. 10−12 cm3 s−1 for typical giant planet thermospheric electron temperatures) [49]. The (2.1) reaction rate is thought to be near its maximum kinetic value [50,51], however, the fraction of molecular hydrogen in the fourth or higher vibrational state is not constrained by observations at present. For Jupiter, we adopt the vibrational density results from calculations by Majeed et al. [52], which lead to an effective H2 vibrational rate coefficient in combination with [51]. For Uranus, as two of the dominant sources of vibrationally excited H2 have been shown to be photon-induced fluorescence and dissociative recombination of H3+ ions [52]—two solar-driven processes—we scale the fractional H2 vibrational populations for Jupiter by (1/r2) to account for the diminution of solar photons with distance. Thus adjusted, the Majeed et al. results are then interpolated onto the appropriate Uranus pressure grid. Further model inputs and specific settings are discussed in relation to their corresponding results in §3.

    (b) Observations and data reduction

    While this is primarily a modelling study, there are two primary sources of data used to constrain the model results at Jupiter. First, the Galileo G0N radio occultation, obtained on 8 December 1995, sampled Jupiter's dusk ionosphere near 24° S latitude and 292° E longitude [53,54]. This measurement provides a representative ionospheric electron density profile suitable for demonstrating the effect of vertical temperature gradients on retrieved H3+ parameters, at least when combined with model simulations that reproduce both the Galileo electron density profile and subsequent H3+ column density observations.

    Such H3+ column densities are the second source of data in this study: ground-based, high-resolution spectroscopic observations of Jupiter in the L telluric window. Jupiter high-resolution (R∼25 000) spectroscopic data spanning 3.26–4 μm were obtained over four nights in April 2016 (14, 16, 20 and 23) using the Near InfraRed Spectrograph (NIRSPEC [55]) on the 10 m Keck II telescope, as described in Moore et al. [56]. Combined, these observations yielded global coverage of Jupiter, with overlapping data from 2+ nights for more than half of the planet. For the current work, a spectrum corresponding to the G0N latitude and local time was extracted and reduced as described in Moore et al. [56]. Briefly, using the line list of Neale et al. [57] and the partition function and total emission formulation of Miller et al. [24], assuming conditions of q-LTE, Gaussian line-fitting techniques [58] are used to fit the modelled spectra to observed H3+ R- and Q-branch lines in order to retrieve column-averaged vibrational temperatures and column-integrated densities [59].

    3. Results and discussion

    The following subsections are dedicated to investigating the errors introduced in retrieved H3+ densities and temperatures due to realistic vertical atmospheric gradients. These errors are expected because standard reduction of H3+ observations, which typically have nadir viewing geometries, starts by assuming a uniform layer of constant density and temperature, whereas we know that neither density nor temperature are constant within giant planet H3+ layers. First, in §3a, we examine a range of realistic temperature and density gradients in order to quantify the degree to which retrieved H3+ column densities are underestimated. Next, in §3b, we produce global simulations of Uranus' ionosphere to demonstrate the effect of gradients on retrieved H3+ temperatures. Finally, in §3c, we combine H3+ spectral observations with electron density constraints to produce a realistic H3+ density profile at Jupiter, and we use the insights from §3a,b in order to extract a preliminary H3+ temperature profile.

    (a) Effect of vertical atmospheric gradients: H3+ density

    We first examine the effect of vertical gradients using a series of simplified, synthetic slab atmospheres. A column-integrated line-of-sight will transect the H3+ layer in a planetary ionosphere at some oblique angle—or along the zenith for observations from directly overhead. Giant planet H3+ emissions are optically thin [60], and lie in a spectral region with strong methane absorption (near 3.4 μm) [61]. The intensity of the H3+ spectrum in each slab scales linearly with density and exponentially with temperature [30]. Furthermore, no self-absorption is considered between slabs. Thus, the observed spectrum follows from the net sum of each individual ‘slab’ emission element within the ionosphere, where the slab emission depends on the column density and temperature of a slab. The retrieved H3+ ‘temperature’ will only represent the ‘true’ atmospheric temperature for an infinitely thin or isothermal H3+ layer. In all other, more realistic cases, the derived temperature and density are thus necessarily a result of the convolution of the altitude structure in both density and temperature.

    Based on previous modelling results [62,63], we begin by approximating Jupiter's H3+ layer to have an equivalent slab width of approximately 1000 km in altitude, as the calculated H3+ density is a nearly constant ∼5000 cm−3 over this altitude range [62]. This leads to an implied column density of ∼0.5 × 1016 m−2. Over this same region, Jupiter's temperature increases from—very approximately—525 K to 775 K [12]. The temperature gradient is thus roughly 0.25 K km−1, and the mean temperature of the H3+ layer, Tmean, is approximately 650 K. We then divide this synthetic H3+ layer into an arbitrary number n slabs, each with column density Nslab = 0.5 × 1016/n, a temperature Tslab based on the temperature gradient, and a corresponding slab emission. We choose n = 10 for this example, as fewer slabs are better-represented graphically, but discuss other variations below. This general atmospheric structure is illustrated in figure 1.

    Figure 1.

    Figure 1. A simplified representation of Jupiter's H3+ layer based on Maurellis & Cravens [62]. (Online version in colour.)

    The composite H3+ spectrum based on the structure from figure 1 will clearly be dominated by higher altitude (redder, hotter) slab spectra. By fitting this spectrum as described in §2b, we retrieve a column density and column-integrated temperature of 0.4 × 1016 m−2 and 695 K, respectively. These values demonstrate an important complication with using observed H3+ spectra to derive realistic column densities, as the retrieved density from the simulated spectral fit is only 80% of the true density. This discrepancy is due to the exponential relation between H3+ temperature and emission, combined with the standard assumption that the ‘H3+ layer’ is isothermal, and this weighting is also reflected in the retrieved temperature, 695 K, which is approximately 7% higher than the true mean temperature represented in figure 1 (650 K).

    A more realistic treatment of the simplified situation depicted above would complicate matters further due to the problems common to all astronomical observations, such as contaminations due to other emissions and absorptions, detector noise, and so forth. Before stressing this initial (Nfit/Ntrue) = 0.8 result too strongly, however, it is important to investigate the various sensitivities that lead to such a discrepancy in the column density. For example, while the conditions in figure 1 have been chosen to roughly represent Jupiter's ionosphere, the mean temperature and the temperature gradient will be different at other locations at Jupiter and at other planets. Similarly, we might question what effect a more realistic density profile would have. The rest of this subsection is therefore devoted to outlining how different choices in representing the atmosphere and in generating the synthetic spectra affect the retrieved H3+ column density.

    First, we examine the choice of the number of slabs n on the derived column density. As demonstrated in figure 2, increasing the number of slabs leads to an asymptotic approach towards (Nfit/Ntrue) = 0.822, which is a approximately 3% increase over the 10 slab representation imagined in figure 1. There is a corresponding inverted trend in the total H3+ emission with slab number, as representations with fewer slabs associate a wider range of the H3+ layer with higher temperatures. In the example considered for figure 2, the net emission decreases by approximately 5% towards an asymptotic value of 21.35 μW m2 sr−1. For comparison, typical giant planet models blanket the ionosphere with fewer than 100 grid points. The derived column density will also depend on the mean temperature in the atmosphere, however. Figure 3 explores this effect while keeping the temperature gradient and the number of slabs fixed (i.e. the temperature gradient is as shown in figure 1, 0.25 K km−1, and the mean temperature is varied).

    Figure 2.

    Figure 2. Total H3+ emission (red squares) and retrievals of H3+ column density Nfit, relative to the true H3+ column density Ntrue (black circles), based on slab spectra for the conditions introduced in figure 1. Dashed lines represents asymptotic values for large slab numbers. (Online version in colour.)

    Figure 3.

    Figure 3. Retrievals of H3+ column density Nfit, relative to the true H3+ column density Ntrue. The number of slabs n is 36, and the mean temperature is varied while holding the vertical temperature gradient constant (i.e. as shown in figure 1, 0.25 K km−1).

    Based on figures 2 and 3, derived H3+ column density is highly dependent upon the atmospheric temperature profile in the H3+ layer. Before exploring a wider range of temperature gradients, such as might be more widely representative of H3+ in giant planet ionospheres, we investigate the additional impact of density gradients. (It should be noted that the Nfit/Ntrue ratio is not sensitive to the true slab column density, which is perhaps not surprising given the linear relation between density and emission and the fractional error introduced from temperature gradients.) Whereas the temperature profile is generally monotonically increasing in the lower thermosphere, there will be both positive and negative gradients in H3+ density.

    Figure 4 presents the combined effects of both density and temperature gradients in the H3+ layer on simulated retrieved column densities, Nfit. These results follow from an atmospheric layer divided into n = 36 slabs, with a true slab column density Ntrue = 0.5 × 1016 m−2 and mean temperatures Tmean of (a) 450 K, (b) 650 K and (c) 850 K. These calculations follow the approach outlined in figure 1, except the absolute H3+ temperatures are allowed to range only between 150 and 1250 K, temperatures observed to be appropriate for giant planet upper atmospheres. This requirement means that, for the largest temperature gradients explored, there are regions of the H3+ layer that are isothermal (i.e. at either 150 K or 1250 K). As expected due to the linear relationship between H3+ emission and density, there is no discrepancy between modelled and ‘observed’ column density when there is no temperature gradient. The generally vertical contours in figure 4 indicate that (Nfit/Ntrue) is most sensitive to temperature gradients, though density gradients begin to play an important role when the temperature gradient is larger than approximately 0.5 K km−1. For realistic temperature gradients at the giant planets, approximately 0.4–2 K km−1 [12,18], the retrieved H3+ column density ranges from approximately 20–90% of the true value, depending primarily on the mean temperature within the H3+ layer. Observed uncertainties in Nfit are often comparable, so the effect illustrated in figure 4 might be absorbed into experimental errors for low S/N data. For instance, in a study of the anticorrelation between H3+ density and temperature, Melin et al. [59] find that a S/N of 12 is required to achieve a 10% column density uncertainty for a column-averaged H3+ temperature of 600 K.

    Figure 4.

    Figure 4. Contours of retrieved H3+ column density Nfit, relative to the true H3+ column density Ntrue, as a function of temperature and density gradients within the H3+ layer. These results follow the approach outlined in figure 1. Specifically, an H3+ layer is divided into n = 36 slabs, with a true slab-column density of 0.5 × 1016 m−2 and with mean temperatures of (a) 450 K, (b) 650 K and (c) 850 K. Within those constraints, synthetic model spectra are then generated based on imposed density and temperature gradients, and Nfit is derived from the composite synthetic spectrum. The dashed line in (b) is described below in §3c structure. (Online version in colour.)

    A more realistic H3+ layer would experience variations in both temperature and density gradients with altitude, and so could not be represented as simply as in figure 4. Careful examination of the evolution of dT/dz and dn/dz with altitude could give some idea of the net error induced in Nfit based on those gradients, though such an analysis would also rely on prior knowledge of the atmosphere, and so significantly reduce the value of added observation. Therefore, while the preceding figures establish that, unless the H3+ layer is in an isothermal atmosphere, the H3+ column density retrieved from observations will represent a lower limit, it is not practical at this stage to examine an all-inclusive range of possible atmospheric structures in order to assign definitive quantitative values to those lower limits. Instead, these results serve as further motivation for investigating other similar complications of interpreting H3+ observations, and for outlining ‘toy’ model parameters that would be relevant for development of a full H3+ retrieval model (e.g. such as in [64]).

    (b) Effect of vertical atmospheric gradients: H3+ temperature (Uranus)

    While density and temperature structures within a planet's H3+ layer affect the retrieved column density, they also affect the interpretation of the retrieved column-averaged temperature. In order to demonstrate this effect, we model the global distribution of H3+ at Uranus. First, the background atmosphere is based on Moses & Poppe [31], appropriate for globally averaged conditions with dust-derived oxygen influxes of 1.2 × 105 H2O molecules cm−2 s−1, 2.5 × 105 CO molecules cm−2 s−1 and 3.0 × 103 CO2 molecules cm−2 s−1, consistent with H2O, CO and CO2 observations [6567]. This 1D atmosphere is applied uniformly at Uranus. While clearly not fully realistic, using a fixed neutral atmosphere, where ion and neutral chemistry is not fully coupled, allows for clearer elucidation of the effects of a varying H3+ layer on retrieved temperatures, as will be demonstrated below. Next, a simulation date of 15 September 2017 is chosen. This choice is largely arbitrary for the purposes of demonstrating the effects of gradients on retrieved H3+ temperatures; however, there do happen to be contemporaneous H3+ observations from September 2017 [68], for which Uranus' sub-solar latitude was −38°. Finally, 1D ionospheric calculations are performed globally as described in §2a, with a 1° latitude resolution. Profiles of background neutral density, temperature and ion density at 30° S latitude are shown in figure 5. Electron temperatures (not shown) are calculated to diverge from the neutral temperature around 2000 km altitude and reach approximately 800 K at 30° S latitude, 12 solar local time (SLT). (Note that altitude levels throughout the text are referenced to the 1 bar pressure level.) Calculated H3+ temperatures are found to be equal to the background neutral temperature.

    Figure 5.

    Figure 5. Model profiles for Uranus. Representative (a) background neutral parameters, which come from Moses & Poppe [31] and are held fixed at all latitudes, and (b) ion density profiles for 30° S latitude, 12 SLT. (Online version in colour.)

    Figure 6 presents global ionospheric density results at Uranus, plotted versus SLT and planetocentric latitude. Note that, as H3+ temperatures were found to be identical to the neutral temperature, at least for the conditions in figure 5, these global simulations do not include plasma temperature calculations, and therefore any H3+ temperature variations are a reflection of the vertical distribution of H3+ and the background temperature profile. Following the preceding approach, Nfit and Tfit follow from generating and fitting synthetic H3+ spectra based on modelled H3+ distributions. These simulations use 138 H3+ slabs (i.e. grid points in altitude), and therefore Nfit is within approximately 0.2% of its asymptotic value based on figure 2.

    Figure 6.

    Figure 6. Global ionospheric model results for Uranus. Contours of (a) the true modelled H3+ column density, Ntrue, (b) the H3+ density retrieved from a fit of the modelled spectrum, Nfit, (c) the Nfit/Ntrue ratio, (d) the column-averaged H3+ temperature from a fit to the modelled spectrum, (e) the altitude of the peak of the H3+ density and (f ) the column-averaged H3+ lifetimes, weighted by H3+ density. In addition, white/black dashed lines in (a,b,f ) indicate the minimum or maximum of each parameter in SLT versus latitude. Solar zenith angle contours, relevant for all panels, are shown in (f ). (Online version in colour.)

    While direct reproduction wasn't the goal of these simulations—and would be at least partially coincidental anyway due to the homogeneous neutral atmosphere—calculated column densities are broadly consistent with H3+ observations [6,68]. As expected from §3a, the true modelled H3+ column density, figure 6a, is slightly larger than the retrieved value, figure 6b, with a global mean column density ratio, figure 6c, of 0.87. The column density ratio is nearest to 1.0 at dawn, and at high northern latitudes, indicating that more of the H3+ layer there is in the isothermal region of the model atmosphere above 2800 km altitude, as would be expected based on solar zenith angle (SZA) effects and nightside recombination chemistry, which depletes lower-altitude layers near the electron density peak more rapidly [69]. This is more evident from figure 6d, which reveals that the column-averaged H3+ temperature is higher in those regions, a direct response of the peak altitude of the H3+ layer being shifted towards higher altitudes, as seen in figure 6e. This is primarily an SZA effect, as more oblique slant paths through the atmosphere generate higher-altitude photoionization; however, there is also a slight offset post-noon due partly to conversion of H+ to H2+ (and thus, H3+, following reaction with H2) via reaction (2.1). The H3+ lifetime at the H3+ peak is on the order of ∼1 h for most of the Uranus day, as would be expected based on the modelled peak density (approx. 3000 cm−3) and dissociative recombination rate of H3+ with electrons (approx. 10−7 cm3 s−1 [70]). Column-averaged H3+ lifetimes, weighted by H3+ density, are slightly higher than at the peak altitude, but still typically less than 3 h (figure 6f ).

    Owing to the simplified nature of the Uranus simulations, and owing to the sparse thermospheric constraints at present, a detailed model-data comparison is not warranted here. Instead, we emphasize the temperature variations shown in figure 6d. Despite the fact that the thermospheric temperature profile is identical at all latitudes, retrieved H3+ temperatures vary by greater than 35 K, or roughly 5% of the mean temperature. This is simply understood as primarily a SZA effect: at low SZAs the H3+ layer is lower in the ionosphere, probing the lower temperatures there, whereas the reverse is true for high SZAs. This result is, again, not surprising; however, it highlights two important points: (1) observed H3+ temperature variations do not necessarily imply anything about the energetics of the thermosphere, and (2) these SZA effects should be accounted for when interpreting measured temperature variabilities. Finally, in reviewing the results of figure 6, it is important to also emphasize the elements that are missing from the simulations presented there. In particular, we have neglected energetic particle precipitation, which would be expected to lead to increased ionization and enhanced thermospheric temperatures, mainly at high magnetic latitudes. At Uranus, the magnetic polar regions are at mid- and low-latitude as a result of the tilted magnetic dipole axis [71], though their exact longitudes are unknown at present due to uncertainty in Uranus' rotation period. Thus, the chief value of figure 6 is in demonstrating the qualitative effects of H3+ density and temperature gradients on retrieved parameters. More realistic global variations of ionization and heating at Uranus would be expected to lead to different quantitative structures.

    (c) Vertical structure of H3+ density and temperature (Jupiter)

    One of the most sought-after observables for planetary atmospheres is the thermal structure, as so much of the rest of planetary dynamics and energetics depend upon it. Furthermore, as established in §3a,b, atmospheric gradients in the H3+ layer can confuse measurements of H3+ density and temperature, meaning the very parameter we want to constrain, T(z), is itself limiting observational insight. We now introduce a potential method of unravelling this confusion by combining column-integrated H3+ measurements with forward modelling.

    For this proof-of-concept, we use the Galileo G0N radio occultation, obtained on 8 December 1995, which sampled Jupiter's dusk ionosphere near 24° S planetocentric latitude and 292° E longitude [54]. Figure 7 presents a model reproduction of the G0N occultation, along with corresponding background neutral atmospheric parameters and modelled ion densities. The model simulation is for 24° S latitude, with a forced vertical drift WD of 50 cm s−1, and plasma density comparisons are extracted for 18 SLT in accordance with the dusk terminator measurement of G0N. As described in §2a, the solar flux is specified using extrapolated TIMED/SEE measurements, the secondary ionization and photoelectron heating rate are parameterized [47], and the effective (2.1) reaction rate comes from Majeed et al. [52], with the rest of the chemistry as specified in Moore et al. [25].

    Figure 7.

    Figure 7. From (ac), Jupiter neutral density and temperature structure; corresponding modelled ion densities; and the Galileo G0N radio occultation electron density profile (black) [54] compared with the model reproduction (red). (Online version in colour.)

    The good model-data agreement in figure 7 gives some confidence that the broad electron density features are well represented (that is, aside from the noisy densities in the topside and the narrow structures below the peak, which might be attributed to gravity waves [72]). As an additional test of the model simulation, we turn to H3+ observations obtained in April 2016 using Keck/NIRSPEC (§2b). First, a series of spectra obtained with the NIRSPEC slit oriented E-W near Jupiter's Great Red Spot are combined and reduced in order to obtain a calibrated H3+ spectrum at 24° S latitude, 17 SLT. Second, we derive a spectral fit to the data, and compare the fit parameters to the H3+ column density from the ionospheric simulation shown in figure 7. An extracted portion of the H3+ spectrum, its corresponding spectral fit, and the modelled column densities are shown in figure 8. Retrieved parameters are Tfit = 774 ± 64 K and Nfit = (2.31 ± 0.88) × 1015 m−2, and the latter is over plotted on the modelled column densities with horizontal error bars that identify the range of local times that contributed to the observed spectrum.

    Figure 8.

    Figure 8. (a) Extracted, calibrated spectral regions from Keck/NIRSPEC observations at Jupiter (grey circles), along with the H3+ spectral fit (orange). (b) Diurnal variation of modelled H3+ column density (black), along with the corresponding value retrieved from the Keck observations (red circle). Vertical error bars come from the spectral fit; horizontal error bars indicate the SLT region over which the spectra were obtained. (Online version in colour.)

    There is good model-data agreement in figure 8b, which lends additional confidence that the model is an accurate representation of both the observed electron density (figure 7) and H3+ column density (figure 8b). Before moving on, however, there are a couple of important caveats to emphasize. First, and perhaps most important, the ionospheric simulation is for 20 April 2016 during solar minimum, in accordance with the Keck H3+ observations, whereas the Galileo radio occultation was obtained on 8 December 1995, nearly 21 years prior, also during solar minimum. It would be far more surprising if Jupiter's ionosphere had not changed in those intervening years than if it had, especially given the variability present in ionospheric radio occultations at giant planets [53,73]. Therefore, the fact that the model agreement is good in both figures 7 and 8 is most likely a coincidence rather than an indication of atmospheric stability, though the fact that the approximately 21 year separation between the two datasets is nearly 2 full solar cycles allows some minimum of hope for a happy coincidence to be maintained. Second, the Galileo radio occultation was at 68° W (System III) longitude, whereas the Keck observations were centred at 308° W longitude, a slightly different magnetic environment. The majority of the modelled H3+ layer is still in photochemical equilibrium, meaning that this variation should have minimal effect, at least if solar photons are the main ionospheric driver as assumed here, though high altitude ion drifts would be altered. Nevertheless, given that the model is able to reproduce both available datasets, and given that there are no better combinations of ionospheric constraints available, we shall progress forward under the assumption that the ionospheric model simulation provides as accurate a representation of the H3+ density structure as possible at present.

    To re-state the problem: the observed H3+ spectrum is a function of the integrated density N(z) and temperature T(z) profiles from the emitting H3+ layer. Thus, the IR spectrum I(λ) also contains information about both of them. Essentially, there are three unknowns, so if either N(z) or T(z) can be convincingly constrained then the other can in principle be derived when combined with the observed spectrum. The spectrum in this example is known from Keck/NIRSPEC observations. Based on the good model-data agreement for the electron density profile and the H3+ column density, we proceed under the assumption that N(z) is appropriately constrained. Therefore, we should also be able to reconstruct at least some limited representation of T(z) from the above two inputs.

    First, as in §3a, the modelled H3+ layer is idealized as n slabs, each of column density Nslab = Ntrue/n. Second, a temperature is assigned to each slab, randomly selected from 150 to 1250 K, and a synthetic H3+ spectrum (i.e. the sum of the slab spectra) is computed. In principle, the specified temperatures of each slab could be independent of each other, but for the present case a monotonically increasing (or isothermal) temperature profile is enforced, consistent with observation [53,74]. This step is repeated until the modelled and observed spectra converge to within some pre-defined tolerance. In this case, the tolerance is set to a maximum of 0.03% disagreement between the slab-modelled spectrum and the spectrum obtained from the H3+ fit. Once converged, it is useful to also provide some estimate of the sensitivity of the result. For this purpose, n-1 slab temperatures are held fixed to their converged values while the other is freely varied until the modelled column-averaged temperature exceeds the measured temperature uncertainty. The derived temperature uncertainties thus represent ∼8% errors in this case, as Tfit = 774 ± 64 K. Finally, based on the converged slab temperatures weighted by their uncertainties, an analytical temperature profile is derived following [75]:

    where z is the altitude element, Texo the exospheric temperature and Ai are constants of the fit.

    Results from the above process with n = 4 are shown in figure 9. Each slab is represented by a shaded column, with varying vertical extent due to enforced equal slab column densities and variable slab number densities. The original modelled H3+ density is given by the red curve. These variable slab widths enable higher altitude resolution near the H3+ density peak and allow for more reasonable error bars than significantly smaller slab widths would. Corresponding slab temperatures are shown as black circles on the right, along with the derived temperature profile in red. Horizontal grey error bars indicate the 1σ uncertainty in slab temperature, and vertical grey error bars demarcate the vertical extent of each slab. The best-fit parameters for the temperature profile in figure 9 are Texo = 788 K, A1 = 136 K, A2 = 159 km, and A3 = 121 km2 K−1, and are only representative of the altitude range with significant H3+ density (e.g. between ∼300 and 2500 km).

    Figure 9.

    Figure 9. (a) The modelled H3+ density (red), represented as 4 slabs of equal column density (and thus varying vertical extent and number density; grey shading). (b) Best-fit slab temperatures (black circles), along with estimated uncertainties (grey lines), and the derived temperature profile (black). (c) Corresponding density (red) and temperature (black) gradients. See text for description of methods. (Online version in colour.)

    Combining the results of figure 9 with those from §3a, it appears that the Keck/NIRSPEC observations were minimally affected by gradients in the H3+ layer. The derived thermal gradient is effectively zero above 1050 km (i.e. less than or equal to 0.05 K km−1), less than 0.2 K km−1 for altitudes between 750 and 1050 km, approximately 0.3 − 1.4 between 550 and 750 km, and approximately 1.8 K km−1 at the bottom side. In situ measurements at Jupiter [12] and Saturn [18] find approximately 2 K km−1 near 400 km and 0.4 K km−1 at the base of the thermosphere, respectively. Based on figure 4, this implies that retrieved H3+ column densities are less than 80% of the modelled values only for altitudes less than 750 km, a range that encompasses 12% of the total column density. Meanwhile, absolute modelled density gradients are less than 8 cm−3 km−1 everywhere, and only approximately 3 cm−3 km−1 where there is also a substantial temperature gradient (near 400 km altitude). Progression of the impact of the derived density and temperature gradients on retrieved column density is shown by the dashed gray curve in figure 4b. In total, the calculated error in column density, based on the combination of figure 4b and the density structure from figure 9, is approximately 10%, all associated with gradients at altitudes below 800 km. This error is comparable to the observational uncertainties (figure 8).

    4. Conclusion

    This study has investigated the effect of atmospheric H3+ density and temperature gradients on the interpretation of observations of the composite spectrum. Overall atmospheric structure is found to cause observed H3+ column densities, Nfit, to represent a lower limit. The degree to which Nfit constrains the true ionospheric column density, Ntrue, depends primarily on the magnitude of any temperature gradients and secondarily on the mean temperature within the H3+ layer. Density gradients can act to reduce the retrieved Nfit/Ntrue ratio even further, provided there is also a temperature gradient present.

    Giant planets generally exhibit strong positive temperature gradients in the lower thermosphere, and consequently, low altitude H3+ is most significantly underestimated. This is also the region where a majority of H3+ is produced, and where the atmosphere is most electrically conductive. Therefore, one immediate caution based on the above results is that nearly all of the error in retrieved H3+ densities due to atmospheric density and temperature gradients is associated with this low altitude region. The total error in H3+ column density from nadir column-integrated observations may be small (e.g. 10%), but the local error in H3+ number density can be large (typically 50% or more, figure 4), and this should be considered when estimating ionospheric electrical conductivities associated with H3+.

    Based on thermal structure in the atmosphere, derived H3+ temperatures are found to represent primarily the temperature at the H3+ density peak. This result, while not surprising, does serve to emphasize that observed H3+ temperature variations may not always represent any inherent evolution in the thermosphere, and may instead be attributed to simple photochemical effects. For example, an X-ray flare would lead to a brief burst of high energy photons, producing a low altitude ionization layer and weighting derived H3+ temperatures towards the (generally) lower temperatures there. Similarly, high SZA regions will absorb ionizing radiation higher in the atmosphere and therefore exhibit higher average H3+ temperatures. This same effect is expected post-sunset, as dissociative recombination with electrons more quickly depletes the lower altitude H3+ near the electron density peak, shifting the effective H3+ layer towards higher and higher altitude throughout the night.

    While this work gives some guidance on the degree to which atmospheric gradients induce additional uncertainty in retrieved column-integrated H3+ densities and temperatures, real-world atmospheric variations far-outstrip those considered here. Therefore, these calculations cannot represent a definitive quantitative manual. Instead, their primary value is in providing qualitative insight, and in demonstrating the potential for enhancing the scientific impact of H3+ observations through complementary modelling studies.

    Finally, by combining two ionospheric datasets with a model simulation that accounts for the above effects, a method for deriving a temperature profile from an overhead H3+ observation is presented. This method relies on the relatively straightforward nature of H3+ solar-driven photochemistry at non-auroral latitudes, and furthermore requires at least some knowledge of the atmospheric structure, so is not easily applicable everywhere. Nevertheless, given the abundance of H3+ observations already obtained, especially at Jupiter and Uranus, it offers potential for improving constraints on global temperature variations at those planets, and it serves as a first-step towards developing a complete H3+ retrieval tool. Ideally, an independent, more traditional method of deriving thermal profiles could first be used to validate this approach. The H3+ limb profiles obtained by the JIRAM instrument [76] on-board the Juno spacecraft at Jupiter may represent the perfect opportunity for such a validation.

    Data accessibility

    The data used herein are available in the public archives or from L.M. upon request.

    Author's contributions

    L.M. led the project, contributed to collection, reduction and analysis of ground-based data, performed ionospheric modelling simulations, produced the figures and wrote the paper. H.M. contributed to collection, reduction and analysis of ground-based data, provided routines for generating synthetic H3+ spectra, and took part in detailed discussions. J.O'D. wrote the Keck observing proposal, contributed to collection, reduction and analysis of ground-based data, provided inputs and advice for interpreting and generating H3+ spectra, and took part in detailed discussions. T.S. contributed to collection, reduction and analysis of ground-based data, including specialized spectral fitting techniques, and took part in detailed discussions. J.M. conducted neutral atmospheric modelling simulations, and provided detailed guidance in incorporating those results in ionospheric calculations. M.G. contributed to development of the secondary ionization and thermal electron heating parameterizations in the modelling. S.M. took part in detailed discussion, and provided expert guidance on the history, chemistry and spectral behaviour of H3+. C.S. contributed to discussion of atmospheric retrieval methods, sensitivities and visualization of results. All authors reviewed and edited the manuscript.

    Competing interests

    We declare we have no competing interests.


    L.M. was supported by the National Aeronautics and Space Administration (NASA) under Grant NNX17AF14G issued through the SSO Planetary Astronomy Program and grant no. 80NSSC19K0546 issued through the Solar System Workings Program. J.M. acknowledges support from NASA Solar System Workings grant no. NNX16AG10G and 80NSSC19K0546. M.G. acknowledges support from STFC of UK under grant ST/N000692/1.


    Some of the data presented herein were obtained at the W. M. Keck Observatory, which is operated as a scientific partnership among the California Institute of Technology, the University of California and the National Aeronautics and Space Administration. The Observatory was made possible by the generous financial support of the W. M. Keck Foundation. We are grateful to the TIMED/SEE PI, Tom Woods, and his team for providing us with the solar flux dataset and associated routines for extrapolation to planets.


    One contribution of 18 to a discussion meeting issue ‘Advances in hydrogen molecular ions: H3+, H5+ and beyond’.

    Published by the Royal Society. All rights reserved.