Role of Planetary Obliquity in Regulating Atmospheric Escape: G-dwarf versus M-dwarf Earth-like Exoplanets

Over the past decade, much attention has been directed toward understanding what factors contribute to exoplanetary habitability (Cockell et al. 2016). In particular, it is widely accepted that orbital parameters play a major role in governing habitability (Shields et al. 2016). One of the chief orbital parameters is the obliquity (axial tilt). The fact that Earth's obliquity is subject to only mild fluctuations is believed to play a vital role in maintaining its stable climate.

As a consequence, numerous studies have analyzed how obliquity affects planetary climate (Spiegel et al. 2009; Ferreira et al. 2014; Rose et al. 2017; Kang 2019). There is broad consensus that climate is sensitive to changes (even minimal ones) in obliquity as it can transition from one stable state to another (Williams & Pollard 2003; Linsenmeier et al. 2015; Kilic et al. 2017; Colose et al. 2019). A number of observational techniques based on amplitude and frequency modulations in light curves arising from thermal emission and scattered light have been proposed for inferring planetary obliquity (Gaidos & Williams 2004; Schwartz et al. 2016; Kane & Torres 2017; Rauscher 2017). Thermal phase curves of Hot Jupiters have already yielded constraints on planetary obliquity; for example, CoRoT-2b has an inferred obliquity of 458 ± 14 (Adams et al. 2019).

Exoplanets around M-dwarfs are typically anticipated to have very low (or zero) obliquities due to rapid tidal energy dissipation (e.g., Heller et al. 2011). This effect may be particularly pronounced for the inner planets of multi-planet systems such as the Kepler-186 system (Bolmont et al. 2014). Nevertheless, there are mechanisms that permit the existence of high-obliquity M-dwarf exoplanets. Perhaps the most famous among them are the "Cassini states" that involve the precession of the planet's spin and orbital angular momenta at the same rate (Colombo 1966; Peale 1969; Winn & Holman 2005); the Moon has an non-zero obliquity of 67 due to this reason. Hence, it is feasible for certain M-dwarf exoplanets to have high obliquities (Dobrovolskis 2009; Wang et al. 2016; Shan & Li 2018; see also Millholland & Laughlin 2019 and Millholland & Batygin 2019).

Another factor that plays a vital role in regulating surficial habitability is the presence of an atmosphere. Moreover, an atmosphere also permits the detection of biosignature gases (e.g., molecular oxygen) via spectroscopy (Kaltenegger 2017; Schwieterman et al. 2018; Madhusudhan 2019). Recent numerical and theoretical studies indicate that both magnetized and unmagnetized planets around M-dwarfs might be particularly susceptible to the depletion of ∼1 bar planetary atmospheres over sub-Gyr timescales due to the high ultraviolet fluxes and intense stellar winds that they experience (Dong et al. 2017b, 2018b, 2018c; Garcia-Sage et al. 2017; Lingam & Loeb 2018, 2019; Airapetian et al. 2019).

In view of the preceding discussion, it is worthwhile to ask the following question: how does obliquity regulate atmospheric escape rates? It is, however, important to note that the rotational and magnetic axes of Earth are separated only by ∼10°. Hence, it is plausible that these two axes are potentially aligned for Earth-like planets as well, viz., the angle between the two might be fairly small. In other words, the magnetic obliquity, i.e., the angle between the magnetic axis and orbital axis, may approximately coincide with the conventional planetary obliquity. Furthermore, as the magnetic obliquity determines the orientation of the planetary magnetic field, it can influence magnetospheric properties and the resultant escape of atmospheric ions, which forms the subject of this Letter.

Thus, in this study, we opt to perform a parametric analysis of how magnetic obliquity affects the atmospheric ion loss from magnetized exoplanets. We focus on two distinct examples due to their astrobiological relevance: an Earth-like planet around a solar-type star, and an Earth-like planet around a late-type M-dwarf using TRAPPIST-1 as a proxy.

A multi-fluid MHD model is utilized herein to simulate the interaction between stellar winds and exoplanets and the concomitant atmospheric ion loss from Earth-like worlds. On such worlds, the major neutral and ion species in planetary upper atmosphere are atomic oxygen and O⁺, respectively (Mendillo et al. 2018), therefore we focus on O⁺—the dominant escaping species—in our model. The multi-fluid MHD model that we develop and employ herein possesses separate continuity, momentum, and pressure equations for each species (see the Appendix for details). Despite the higher computational expense incurred by the multi-fluid MHD approach, they are more realistic and accurate than the traditional single-fluid approach (Tóth et al. 2012; Dong et al. 2014, 2018a).

To summarize, (1)–(14) in the Appendix permit us to simulate the O⁺ ions, stellar wind protons (H⁺), and electrons (e⁻) with complete self-consistency; all essential interactions between fluids are accounted for, as well as their individual evolution. More specifically, both elastic and inelastic collisions and the associated heating and cooling terms are explicitly present in our model. It is important to recognize that in addition to photoelectron heating, our model also incorporates Joule heating, as seen from manipulating the second term in the square brackets on the rhs of (4) along the following lines.

$$\begin{eqnarray*}&&\sum _{t={\rm{O}},{{\rm{O}}}^{+}}\displaystyle \frac{{\rho }_{e}{\nu }_{{et}}}{{m}_{e}+{m}_{t}}\left[\displaystyle \frac{2}{3}{m}_{t}{\left({{\boldsymbol{u}}}_{{\rm{t}}}-{u}_{e}\right)}^{2}\right]\propto \displaystyle \frac{{{\boldsymbol{J}}}^{2}}{{\sigma }_{e}}={\boldsymbol{J}}\cdot {\boldsymbol{E}}\end{eqnarray*}$$

where m_t/(m_e + m_t) ≈ 1 in view of the fact that the mass of oxygen is much larger than the electron mass.

A number of chemical reactions such as photoionization, electron impact ionization, charge exchange between neutrals and ions, and recombination are included in the form of source terms delineated in the Appendix. Additional details concerning the reactions are outlined in Table 1. For a brief discussion of the utility of multi-fluid MHD models in our solar system, we refer the reader to Section 2 of Dong et al. (2017a).

Table 1.  The Elastic Collision Rates and Chemical Reaction Rates Used in the Multi-fluid MHD Code

Elastic Collision	Rate (s−1)	References
ion–ion (νst)	$1.27\times \tfrac{{Z}_{s}^{2}{Z}_{t}^{2}\sqrt{{m}_{{st}}}}{{m}_{s}}\tfrac{{n}_{t}}{{T}_{{st}}^{3/2}}$	Schunk & Nagy (2009)
ion–neutral (νsn)	Csnnna	Schunk & Nagy (2009)
electron–ion (νes)	$54.5\times \tfrac{{n}_{s}{Z}_{s}^{2}}{{T}_{e}^{3/2}}$	Schunk & Nagy (2009)
ion–electron (νse)	$1.27\times \tfrac{\sqrt{{m}_{e}}}{{m}_{s}}\tfrac{{n}_{e}{Z}_{s}^{2}}{{T}_{e}^{3/2}}$	Schunk & Nagy (2009)
electron–neutral (νen)	$8.9\times {10}^{-11}{n}_{n}(1+5.7\times {10}^{-4}{T}_{e}){T}_{e}^{1/2}$	Schunk & Nagy (2009)
Chemical Reaction
Primary Photolysis and Particle Impact	Rate (s−1)
${\rm{O}}+h\nu \to {{\rm{O}}}^{+}+{e}^{-}$	See table footnoteb	Schunk & Nagy (2009)
${{\rm{e}}}^{-}+{\rm{O}}\to {{\rm{e}}}^{-}+{{\rm{O}}}^{+}+{e}^{-}$	See reference	Cravens et al. (1987)
Ion–Neutral Chemistry	Rate (cm3 s−1)
${{\rm{H}}}^{+}+{\rm{O}}\to {\rm{H}}+{{\rm{O}}}^{+}$	3.75 × 10−10	Schunk & Nagy (2009)
Electron Recombination Chemistry	Rate (cm3s−1)
${{\rm{O}}}^{+}+{e}^{-}\to {\rm{O}}$	$3.7\times {10}^{-12}{\left(\tfrac{250}{{T}_{e}}\right)}^{0.7}$,	Schunk & Nagy (2009)
${{\rm{H}}}^{+}+{e}^{-}\to {\rm{H}}$	$4.8\times {10}^{-12}{\left(\tfrac{250}{{T}_{e}}\right)}^{0.7}$,	Schunk & Nagy (2009)

Notes. For elastic collisions, Z_s, m_s, n_s, and T_s denote the charge state, mass (in amu), number density (in cm⁻³), and temperature (in K) for a given species. Note that ${m}_{{st}}=\tfrac{{m}_{s}{m}_{t}}{{m}_{s}+{m}_{t}}$ and ${T}_{{st}}=\tfrac{{m}_{s}{T}_{t}+{m}_{t}{T}_{s}}{{m}_{s}+{m}_{t}}$ represent the reduced mass and temperature.

^aBoth ${\nu }_{{{\rm{H}}}^{+}{\rm{O}}}$ and ${\nu }_{{{\rm{O}}}^{+}{\rm{O}}}$ are for resonant ion–neutral interactions, where ${\nu }_{{{\rm{H}}}^{+}{\rm{O}}}=6.61\times {10}^{-11}{n}_{{\rm{O}}}{T}_{{{\rm{H}}}^{+}}^{1/2}{(1-0.047{\mathrm{log}}_{10}{T}_{{{\rm{H}}}^{+}})}^{2}$ for ${T}_{r}\gt 300\,{\rm{K}}$, and ${\nu }_{{{\rm{O}}}^{+}{\rm{O}}}\,=3.67\times {10}^{-11}{n}_{{\rm{O}}}{T}_{r}^{1/2}{(1-0.064{\mathrm{log}}_{10}{T}_{r})}^{2}$ for T_r > 235 K. T_r = (T_i+T_n)/2, where T_i and T_n are the ion and neutral temperatures, respectively. Atomic oxygen number density n_O has units of cm⁻³. ^bThe photoionization rate has been appropriately rescaled based on the extreme ultraviolet (EUV) values at Earth and TRAPPIST-1g. The photoionization rate equals 4.13 × 10⁻⁷ s⁻¹ for Earth, and 32.96 × 10⁻⁷ s⁻¹ for TRAPPIST-1g.

Download table as:  ASCIITypeset image

As mentioned in Section 1, we investigate two different Earth-like planets by changing the host star that they orbit; this endeavor is important because stellar parameters influence numerous aspects of habitability (Lingam & Loeb 2019). The first involves a G-type star with the Sun constituting the proxy. The second is based on TRAPPIST-1 because it is a well-known example of late M-dwarfs.

To investigate Earth-like planets, we make use of Earth's thermospheric temperature profile (Schunk & Nagy 2009) and specify a fiducial surface pressure of 1 bar and magnetic dipole moment equal to the Earth's. It is worth noting that Earth's thermosphere is much hotter (∼1000 K) than the surface and lower atmosphere. Hence, a surface temperature change of ∼20 K induced by obliquity as per climate models (Kang 2019) is not likely to affect the upper atmosphere significantly. As noted previously in Section 2, O⁺ is the primary ion undergoing atmospheric escape in our model because we are interested in exoplanets that resemble Earth. It is feasible to select alternative atmospheres (e.g., resembling Venus), along the lines of Dong et al. (2017b, 2018b), and study the effects of varying magnetic obliquity but this is left for future publications.

Our simulation domain commences at 100 km, where the density of O⁺ ions is roughly in photochemical equilibrium. In our numerical model, float boundary conditions are employed for the velocity u and the magnetic field B. The simulation box is extended up to 200 planetary radii by means of a non-uniform spherical grid. In the lower ionosphere and thermosphere, the lowest spatial resolution of ∼10 km—several times smaller than the thermospheric scale height—is attained to capture fine-scale variations in the upper atmosphere. The angular (i.e., horizontal) resolution is 3^◦ × 3°. The equations in the Appendix are solved by means of an upwind finite-volume scheme (Tóth et al. 2012); see Dong et al. (2017a) for further details.

Because we are studying Earth-like planets around solar-type stars and late M-dwarfs, we require the appropriate stellar wind parameters to compute self-consistent ion escape rates. For the Sun, current solar wind parameters are adopted from Schunk & Nagy (2009). Given that we use TRAPPIST-1 as a proxy, we use the simulated stellar wind parameters from Dong et al. (2018b) at TRAPPIST-1g. The parameters for the two cases are presented in Table 2. The Planet-Star-Orbital coordinates are used herein. In this system, the X-axis is directed from the planet toward the star, the Z-axis is normal to the planet's orbital plane, and the Y-axis is perpendicular to the X- and Z-axes.

Recall that our primary intention is to determine how planetary obliquity regulates the atmospheric ion loss from a magnetic perspective. The final results are depicted in Table 2. There are several interesting features that stand out in both cases, viz., solar-type stars and late M-dwarfs.

First, the atmospheric ion escape rate is ${ \mathcal O }({10}^{26}{s}^{-1})$ in the case of G-dwarf Earth-like planets, whereas the escape rate increases to ${ \mathcal O }({10}^{28}\,{{\rm{s}}}^{-1})$ for M-dwarf planets due to the extreme stellar wind conditions and high-energy radiation in the close-in habitable zones (HZs). In other words, a ∼1 bar atmosphere of an Earth-like planet would take ${ \mathcal O }({10}^{10})$ yr to be depleted for a G-type star and ${ \mathcal O }({10}^{8})$ yr for an M-dwarf based on normal stellar wind conditions.

Second, the variation in the atmospheric ion escape rates is virtually independent of the magnetic obliquity for an Earth-analog around a solar-type star. We find that the total variation is less than 10%. In contrast, when we consider an Earth-like planet around a late M-dwarf, we determine that the variation is modest (but non-negligible); in quantitative terms, the maximum escape rate is more than twice (or 200%) the minimum value. The chief reason why the obliquity plays a weak role in determining the escape rate for solar-type stars stems from the temperate stellar wind and radiation in HZs.

As shown in Figure 1, the magnetosphere of the G-dwarf planet is larger than that of the M-dwarf planet; therefore, regardless of magnetic obliquity's value, the ionosphere does not experience much difference. On the other hand, for an Earth-analog around TRAPPIST-1, the dual effect of the compressed magnetosphere and high-energy radiation makes the ion sources (e.g., electron impact ionization and charge exchange) more sensitive to the magnetic configuration.

Zoom In Zoom Out Reset image size

Third, we see that the maximal ion escape rate is attained at a magnetic obliquity of 90°, whereas the minimum occurs at 0° or 180° (Figure 2). While the atmospheric escape rates at 0° and 180° are nearly the same, there is a clear distinction between 90° and 270°. The reason behind the latter behavior has to do with the relative orientation of the interplanetary magnetic field (IMF) and the planetary magnetic field. At the magnetic obliquity of 90°, the IMF can directly connect to the dayside planetary surface due to the field polarity, whereas the IMF can only connect to the nightside surface at the magnetic obliquity of 270°; see the third column of Figure 1. Therefore, stellar wind particles, especially electrons (with relatively low energy) can be transported along the field lines and ionize the atomic oxygen in the upper atmosphere via impacts as shown in Figure 3.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 1 depicts the contour plots of the O⁺ ion density for an Earth-like planet around a solar-type star (top two panels) and a late M-dwarf (bottom two panels) at different values of the magnetic obliquity. Consistent with the results in Table 2, more O⁺ ions escape from the M-dwarf planet due to the extreme stellar wind and radiations. Figure 2 shows the corresponding atmospheric oxygen ion escape rates, where two peaks clearly occur at magnetic obliquities of 90° and 270°; the former is slightly higher than the latter.

In the 90° or 270° cases, the cusp region (i.e., the region filled with open magnetic field lines) directly faces the star and thus the stellar wind. Therefore, the stellar wind particles can impact the dayside upper atmosphere relatively easily and deposit their energy. Compared to the 270° case, the 90° case possesses a better magnetic connection with the planet (based on our model set-up), thereby achieving a slightly higher escape rate primarily due to the high electron impact ionization.

In order to demonstrate this point, the oxygen ion density and velocity (top row) and different ionization rates (second to fourth rows from the top) are illustrated for an Earth-like planet around a late M-dwarf in Figure 3. The ion outflow of O⁺ (namely, the polar wind) driven primarily by the electron pressure gradient, ∇p_e, is rendered in the top row of Figure 3. The bottom row of Figure 3 shows the difference of the total ionization rate between 90° and 270° (left column), resembling the difference in electron impact ionization rate (right column); in other words, the difference in total ionization rate is mainly regulated by electron impact ionization.

Lastly, over an extended period of time, due to the polarity reversals of stellar and planetary magnetic fields (Glatzmaier 2013), variations in the escape rate with obliquity may get averaged out.

In this Letter, we have described a sophisticated numerical code for simulating the escape of atmospheric ions from Earth-like (exo)planets. We have incorporated the appropriate upper-atmospheric chemistry and evolve each species separately via the multi-fluid MHD equations of Section 2. This model was applied to study two different exoplanets: the first around a solar-type star and the second orbiting a late M-dwarf, for which TRAPPIST-1 was used as a proxy. Our goal was to determine how the atmospheric ion escape rates vary with the planetary obliquity from a magnetic perspective.

We found that the maximum escape rates arose at obliquities of 90° or 270° (depending on field polarities), whereas the minimum rates were attained at 0° or 180°. The reason is that the cusp (comprising open field lines) faces the stellar wind at obliquities of 90° or 270°, and allows the stellar wind particles to deposit their energy in the planetary upper atmosphere. For Earth-like planets around solar-type stars, it is found that the escape rate is virtually independent of the obliquity. On the other hand, for late M-dwarfs, we determined that the escape rate varies by more than a factor of ∼2.

From our simulations, we found that the timescale required to deplete a ∼1 bar Earth-like atmosphere is ${ \mathcal O }({10}^{10})$ and ${ \mathcal O }({10}^{8})$ yr, for solar-type stars and late M-dwarfs, respectively. If we assume that the source of atmospheric O is water from oceans, we find that the mass of Earth's oceans (M_oc) cannot be depleted over the main-sequence lifetime of a solar-type star. In contrast, for a late M-dwarf we determine that M_oc could be depleted over a timescale of ${ \mathcal O }({10}^{10})$ yr, which is shorter than the star's lifetime.

There are two conclusions to be drawn from this finding. When studying Earth-like planets around solar-type stars, at least insofar as the atmospheric ion escape rates are concerned, the effects of obliquity is ostensibly minimal. In contrast, when it comes to exoplanets around late M-dwarfs, obliquity might play an important role. There are, however, two different cases to consider for M-dwarf exoplanets. In the first case, if the atmosphere is completely depleted over a timescale that is orders of magnitude smaller than 1 Gyr, changing this value by a factor of ∼2 will probably not have major implications for the origin and evolution of life.

However, let us consider the second case, wherein the M-dwarf exoplanet under question has a sufficiently massive atmosphere or rapid outgassing to permit the retention of an atmosphere over a few Gyr. The habitable timescale for biological evolution will be roughly halved as one moves from an obliquity of 0° to an obliquity of 90°. It is therefore instructive to carry out a thought experiment. Suppose that an Earth-like planet can retain an atmosphere for 4 Gyr for an obliquity of 0° and that biological evolution unfolds in a similar fashion as on Earth;⁵ at an obliquity of 90°, the depletion timescale is 1.9 Gyr. For this specific hypothetical planet, at 0° obliquity, enough time might exist for the emergence of complex multicellular organisms, whereas an obliquity of 90° may not suffice for the evolution of eukaryotic analogs and complex multicellularity (Knoll 2015).

Some caveats regarding our treatment are worth emphasizing here. We have chosen to vary only the stellar parameters, thus leaving planetary parameters (e.g., size and magnetic field strength) fixed. In actuality, habitable exoplanets are probably very diverse and the escape rates will change accordingly; however, it is plausible that the qualitative trends described herein are still partly valid. Moreover, we have incorporated only O⁺ as it constitutes the dominant ion species in the Earth's ionosphere, but subsequent treatments should incorporate additional minor species. As we utilize a multi-fluid MHD model, kinetic effects contributing to atmospheric ion escape are not included in our model (e.g., see Strangeway et al. 2005). Finally, the issue of atmospheric depletion is difficult to address comprehensively because it also necessitates knowledge of other pertinent issues including outgassing, bolide impacts, (super)flares, and associated phenomena (e.g., coronal mass ejections).

To summarize, planetary magnetic obliquity does not appear to affect atmospheric ion escape rates for habitable planets around solar-type stars, whereas it has a weak-to-moderate influence on the escape rates for late M-dwarf exoplanets.

The authors acknowledge fruitful discussions with Yutong Shan, Anthony Del Genio, Michael Way, Fei Dai, Gongjie Li, and Joshua Winn. C.D. was supported by NASA grant 80NSSC18K0288. M.L. was supported by the Institute for Theory and Computation at Harvard University. Resources for this work were provided by the NASA High-End Computing (HEC) Program through the NASA Advanced Supercomputing (NAS) Division at Ames Research Center. The Space Weather Modeling Framework that contains the BATS-R-US code used in this study is publicly available from http://csem.engin.umich.edu/tools/swmf. For distribution of the model results used in this study, please contact the corresponding author.

In this section, we describe the multi-fluid MHD model, endowed with the electron pressure equation, that is used to simulate the oxygen ion loss from Earth-like exoplanets.

The multi-fluid MHD model comprises three fluids, two of which are ionic O⁺ and H⁺ (with subscript s) and the last is the electron fluid with subscript e. For the background neutral species, we employ the subscript n. In the multi-fluid MHD equations, ρ, u, p, $\overleftrightarrow{I}$, k_B, and γ = 5/3 represent the mass density, velocity vector, pressure, identity matrix, Boltzmann constant, and specific heat ratio, respectively. An extended description of the multi-fluid MHD equations can be found in Rubin et al. (2014), Huang et al. (2016), and Dong et al. (2017a):

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{\partial {\rho }_{s}}{\partial t}+{\rm{\nabla }}\cdot ({\rho }_{s}{{\boldsymbol{u}}}_{{\boldsymbol{s}}})={{ \mathcal S }}_{s}-{{ \mathcal L }}_{s}\end{array}\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{\partial \left({\rho }_{s}{{\boldsymbol{u}}}_{{\boldsymbol{s}}}\right)}{\partial t}+{\rm{\nabla }}\cdot \left({\rho }_{s}{{\boldsymbol{u}}}_{{\boldsymbol{s}}}{{\boldsymbol{u}}}_{{\boldsymbol{s}}}+{p}_{s}\overleftrightarrow{I}\right)={n}_{s}{q}_{s}\left({{\boldsymbol{u}}}_{{\boldsymbol{s}}}-{{\boldsymbol{u}}}_{+}\right)\times {\boldsymbol{B}}\\ \,+\,\displaystyle \frac{{q}_{s}{n}_{s}}{{{en}}_{e}}\left({\boldsymbol{J}}\times {\boldsymbol{B}}-{\rm{\nabla }}{p}_{e}\right)+{\rho }_{s}{\boldsymbol{G}}\\ \,+\,{\rho }_{s}\sum _{t=\mathrm{all}}{\nu }_{{st}}({{\boldsymbol{u}}}_{{\boldsymbol{t}}}-{{\boldsymbol{u}}}_{{\boldsymbol{s}}})+{{ \mathcal S }}_{s}{{\boldsymbol{u}}}_{{\boldsymbol{n}}}-{{ \mathcal L }}_{s}{{\boldsymbol{u}}}_{{\boldsymbol{s}}}\end{array}\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}\begin{array}{rcl} & & \displaystyle \frac{\partial {p}_{s}}{\partial t}+\left({{\boldsymbol{u}}}_{s}\cdot {\rm{\nabla }}\right){p}_{s}=-{\gamma }_{s}{p}_{s}\left({\rm{\nabla }}\cdot {{\boldsymbol{u}}}_{{\boldsymbol{s}}}\right)\\ & & \quad +\sum _{t=\mathrm{all}}\displaystyle \frac{{\rho }_{s}{\nu }_{{st}}}{{m}_{s}+{m}_{t}}\left[2{k}_{B}\left({T}_{t}-{T}_{s}\right)+\displaystyle \frac{2}{3}{m}_{t}{\left({{\boldsymbol{u}}}_{t}-{{\boldsymbol{u}}}_{{\boldsymbol{s}}}\right)}^{2}\right]\\ & & \quad +{k}_{B}\displaystyle \frac{{{ \mathcal S }}_{s}{T}_{n}-{{ \mathcal L }}_{s}{T}_{s}}{{m}_{s}}+\displaystyle \frac{1}{3}{S}_{s}{\left({{\boldsymbol{u}}}_{n}-{{\boldsymbol{u}}}_{{\boldsymbol{s}}}\right)}^{2}\end{array}\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{\partial {p}_{e}}{\partial t}+\left({{\boldsymbol{u}}}_{e}\cdot {\rm{\nabla }}\right){p}_{e}=-{\gamma }_{e}{p}_{e}\left({\rm{\nabla }}\cdot {{\boldsymbol{u}}}_{e}\right)\\ +\,\sum _{t={\rm{s}},{\rm{n}}}\displaystyle \frac{{\rho }_{e}{\nu }_{{et}}}{{m}_{e}+{m}_{t}}\left[2{k}_{B}\left({T}_{t}-{T}_{e}\right)+\displaystyle \frac{2}{3}{m}_{t}{\left({{\boldsymbol{u}}}_{{\boldsymbol{t}}}-{{\boldsymbol{u}}}_{{\boldsymbol{e}}}\right)}^{2}\right]\\ -\,{k}_{B}\displaystyle \frac{{{ \mathcal L }}_{e}{T}_{e}}{{m}_{e}}+\displaystyle \frac{2}{3}{n}_{n}({\nu }_{\mathrm{ph},n}{{ \mathcal E }}_{\mathrm{ns}}^{\mathrm{exc}}-{\nu }_{\mathrm{imp},n}{{ \mathcal E }}_{\mathrm{ns}}^{\mathrm{pot}})\\ -\,\displaystyle \frac{2}{3}{n}_{e}{n}_{n}{{ \mathcal R }}_{{en}}^{\mathrm{inelastic}}+\displaystyle \frac{1}{3}{S}_{e}{\left({{\boldsymbol{u}}}_{{\boldsymbol{n}}}-{{\boldsymbol{u}}}_{{\boldsymbol{e}}}\right)}^{2}\end{array}\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\boldsymbol{B}}}{\partial t}={\rm{\nabla }}\times ({{\boldsymbol{u}}}_{+}\times {\boldsymbol{B}}-\eta {\boldsymbol{J}})\,\end{eqnarray} \tag{ 5 }$$

where ν signifies the collision frequency between two species and ${{\boldsymbol{u}}}_{+}$ refers to the charge-averaged velocity

$$\begin{eqnarray}&&{{\boldsymbol{u}}}_{+}=\sum _{s=\mathrm{ions}}\displaystyle \frac{{q}_{s}{n}_{s}{u}_{s}}{{{en}}_{e}},\end{eqnarray} \tag{ 6 }$$

while η is the magnetic diffusivity

$$\begin{eqnarray}&&\eta =\displaystyle \frac{1}{{\mu }_{0}{\sigma }_{e}}=\displaystyle \frac{1}{{\mu }_{0}}\left(\displaystyle \frac{1}{{\sigma }_{{en}}}+\displaystyle \frac{1}{{\sigma }_{{ei}}}\right),\end{eqnarray} \tag{ 7 }$$

where σ_e is the electron conductivity and is composed of both electron–neutral (${\sigma }_{{en}}={e}^{2}{n}_{e}/{{\rm{\Sigma }}}_{n^{\prime} }{\nu }_{{en}^{\prime} }{m}_{e}$) and electron–ion (${\sigma }_{{ei}}={e}^{2}{n}_{e}/{{\rm{\Sigma }}}_{s^{\prime} }{\nu }_{{es}^{\prime} }{m}_{e}$) collisions.

The preceding set of equations consists of source (${ \mathcal S }$) and loss (${ \mathcal L }$) terms:

$$\begin{eqnarray}&&{{ \mathcal S }}_{s}={m}_{s}{n}_{s^{\prime} }\left({\nu }_{\mathrm{ph},s^{\prime} }+{\nu }_{\mathrm{imp},s^{\prime} }+\sum _{i=\mathrm{ions}}{k}_{{is}^{\prime} }{n}_{i}\right)\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&{{ \mathcal L }}_{s}={m}_{s}{n}_{s}\left({\alpha }_{R,s}{n}_{e}+\sum _{t^{\prime} =\mathrm{neutrals}}{k}_{{st}^{\prime} }{n}_{t^{\prime} }\right)\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&{{ \mathcal S }}_{e}={m}_{e}\sum _{s^{\prime} }\,{n}_{s^{\prime} }({\nu }_{\mathrm{ph},s^{\prime} }+{\nu }_{\mathrm{imp},s^{\prime} })\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}&&{{ \mathcal L }}_{e}={m}_{e}{n}_{e}\sum _{s=\mathrm{ions}}{\alpha }_{R,s}{n}_{s}.\end{eqnarray} \tag{ 11 }$$

Note that the source and loss terms consist of photoionization (${\nu }_{\mathrm{ph},s^{\prime} }$), charge exchange (${k}_{{is}^{\prime} }$), electron impact ionization (${\nu }_{\mathrm{imp},s^{\prime} }$), and recombination (α_R,s).

In addition, inelastic collisions between electrons and O (neutral oxygen) are effective at cooling the former in the lower thermosphere, where collisions occur regularly. Therefore, we incorporate the cooling rate coefficient ${{ \mathcal R }}_{{en}}^{\mathrm{inelastic}}$ (in eV cm³ s⁻¹) in (4) by adopting the formalism in Schunk & Nagy (2009):

$$\begin{eqnarray}\begin{array}{rcl}{{ \mathcal R }}_{{en}}^{\mathrm{inelastic}} & = & {D}^{-1}\{{S}_{10}\left(1-\exp \left[98.9\left({T}_{e}^{-1}-{T}_{n}^{-1}\right)\right]\right)\\ & & +{S}_{20}(1-\exp \left[326.6\left({T}_{e}^{-1}-{T}_{n}^{-1}\right)]\right)\\ & & +{S}_{21}(1-\exp \left[227.7\left({T}_{e}^{-1}-{T}_{n}^{-1}\right)]\right)\},\end{array}\end{eqnarray} \tag{ 12 }$$

where $D=5+\exp (-326.6{T}_{n}^{-1})+3\exp (-227.7{T}_{n}^{-1})$, S₂₁ = 1.863 · 10⁻¹¹, S₂₀ = 1.191 · 10⁻¹¹, and ${S}_{10}=8.249\cdot {10}^{-16}{T}_{e}^{0.6}\exp (-227.7{T}_{n}^{-1})$.

The optical depth of the neutral atmosphere is determined by employing the numerical formula in Smith & Smith (1972) to study photoionization. Photoelectrons gain excess energy ${{ \mathcal E }}_{\mathrm{ns}}^{\mathrm{exc}}$ during the photoionization process and they lose the ionization energy of neutral oxygen, ${{ \mathcal E }}_{\mathrm{ns}}^{\mathrm{pot}}$, during electron impact ionization (Schunk & Nagy 2009), as seen from Equation (4). The number density and velocity of electrons is easy to calculate after imposing quasineutrality and expressing it in terms of the current (Tóth et al. 2012)

$$\begin{eqnarray}&&{n}_{e}=\displaystyle \frac{1}{e}\sum _{s=\mathrm{ions}}{n}_{s}{q}_{s},\end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray}&&{{\boldsymbol{u}}}_{e}={{\boldsymbol{u}}}_{+}-\displaystyle \frac{{\boldsymbol{J}}}{{{en}}_{e}}={{\boldsymbol{u}}}_{+}-\displaystyle \frac{{\rm{\nabla }}\times {\boldsymbol{B}}}{{\mu }_{0}{{en}}_{e}},\end{eqnarray} \tag{ 14 }$$

with the last equation following from Ampére's law.