Planetary Magnetic Field Control of Ion Escape from Weakly Magnetized Planets

TL;DR

This study uses 3D hybrid simulations to explore how weak planetary magnetic fields influence ion escape, revealing a balance where magnetic strength initially increases escape but then inhibits it, depending on stellar wind conditions.

Contribution

It demonstrates the complex role of weak magnetic fields in ion escape, providing new insights into planetary atmospheric loss mechanisms and their dependence on magnetic and stellar wind parameters.

Findings

Ion escape rate peaks at an optimal magnetic field strength.

Stronger magnetic fields can both enhance and inhibit ion escape.

Power-law relationships describe how escape parameters vary with magnetic field.

Abstract

Intrinsic magnetic fields have long been thought to shield planets from atmospheric erosion via stellar winds; however, the influence of the plasma environment on atmospheric escape is complex. Here we study the influence of a weak intrinsic dipolar planetary magnetic field on the plasma environment and subsequent ion escape from a Mars sized planet in a global three-dimensional hybrid simulation. We find that increasing the strength of a planet's magnetic field enhances ion escape until the magnetic dipole's standoff distance reaches the induced magnetosphere boundary. After this point increasing the planetary magnetic field begins to inhibit ion escape. This reflects a balance between shielding of the southern hemisphere from ``misaligned" ion pickup forces and trapping of escaping ions by an equatorial plasmasphere. Thus, the planetary magnetic field associated with the peak ion escape rate is critically dependent on the stellar wind pressure. Where possible we have fit power laws for the variation of fundamental parameters (escape rate, escape power, polar cap opening angle and effective interaction area) with magnetic field, and assessed upper and lower limits for the relationships.

Tables

Table 1. Table 1: Upstream conditions used to drive the models. Here u 𝑢 u is the stellar wind velocity, n p subscript 𝑛 𝑝 n_{p} is the stellar wind H + superscript 𝐻 H^{+} number density, T p subscript 𝑇 𝑝 T_{p} is the stellar wind H + superscript 𝐻 H^{+} temperature, B s ​ w subscript 𝐵 𝑠 𝑤 B_{sw} is the interplanetary magnetic field (IMF), P s ​ w subscript 𝑃 𝑠 𝑤 P_{sw} is the dynamical pressure, and v A subscript 𝑣 𝐴 v_{A} is the Alfven speed.

$u$	$[-350,0,0]\mathrm{km}/s$
$n_p$	4.9 ${\mathrm{cm}}^{-3}$
$T_p$	59200 K
$B_{sw}$	$[-0.74,5.46,-0.97]\mathrm{nT}$
$P_{sw}$	1.0 nPa
$v_A$	55 km/s

Table 2. Table 2: Summary of derived scaling relationships for polar cap solid angle ( Ω Ω \Omega ), mass outflow ( M ˙ ˙ 𝑀 \dot{M} ), escape power ( P e ​ s ​ c subscript 𝑃 𝑒 𝑠 𝑐 P_{esc} ), and bow shock interaction area ( A 𝐴 A ) with magnetic field along with upper and lower limits of the scaling law. Lower limits are calculated empirically from the simulations while upper limits are postulated as in the text.

Parameter	Scaling	Lower Limit	Upper Limit
$\Omega$	$\propto B^{-0.33}$		N/A
$\dot{M}$	$\propto B^{-0.67}$	$R_s\sim R_{IMB}$, $B\sim B_{max}$	a) may flatten beyond $P_B\sim 35P_{sw}$ b) cusp dynamics may dominate past $P_B\sim {\left(5.26\right)}^6{P}_{sw}$
$P$	$\propto B^{-0.60}$	$R_s\sim R_{IMB}$, $B\sim B_{max}$	a), b)
$A$	$\propto B^{0.41}$	$R_s\sim R_{IMB}$, $B\sim B_{max}$	Bow shock models needed

Equations

$$m_i\frac{d\overrightarrow{v_i}}{dt}=q_i(\vec{E}+\overrightarrow{v_i}\times \vec{B}),$$

$$\vec{J}=\nabla \times \vec{B}/{\mu }_0,$$

$$\vec{E}=-\overrightarrow{U_e}\times \vec{B}+\eta \vec{J},$$

$$\frac{\partial \vec{B}}{\partial t}=-\nabla \times \vec{E}.$$

$$q(\chi ,z^{\prime })=Q_0\exp [1-z^{\prime }-\sec (\chi )\ast e^{-z^{\prime }}],$$

$$R_s=R_P{\left(\frac{P_B}{P_{sw}}\right)}^{1/6}=R_P{\left(\frac{2B_P^2}{{\mu }_0{\rho }_{sw}{u}_{sw}^2}\right)}^{1/6}$$

$${\theta }_c=\arcsin \left(\sqrt{\frac{R_p}{R_s}}\right)$$

$${\Omega }_c=4\pi (1-\cos ({\theta }_c))=4\pi (1-\sqrt{1-R_p/R_s}).$$

$$\Phi =\int n_i(\overrightarrow{u_i}\cdot \vec{\hat{r}})r^{2}d\Omega$$

$$P_{esc}=\int m_i{n}_i{u_i}^2(\overrightarrow{u_i}\cdot \vec{\hat{r}})r^{2}d\Omega$$

$$\bar{E}=P_{esc}/(\mu \Phi )$$

$${\dot{M}}_{in}=A({\rho }_{sw}{u}_{sw})$$

$${\dot{P}}_{in}=A({\rho }_{sw}{u}_{sw}^3)$$

Full text

1 Introduction

Atmospheric escape is capable of shaping a planet’s atmospheric composition and total mass, and thus the planet’s long-term habitability. Loss to space of atmospheric particles is thought to have played a key role in the atmospheric evolution of both Mars (e.g. Anders and Owen, 1977; Jakosky and Phillips, 2001; Jakosky et al., 2018) and Venus (e.g. Dayhoff et al., 1967; Kulikov et al., 2006). This evolution may have prevented these planets from being currently habitable like the Earth, despite their likely similarities in initial formation. As we begin to discover an abundance of exoplanets, it is important to understand the processes that shape their evolution, including atmospheric escape. This will be particularly interesting for potentially habitable terrestrial planets, such as Proxima-b (Anglada-Escudé et al., 2016; Meadows et al., 2018) and those in the Trappist-1 system (Gillon et al., 2017; Dong et al., 2018).

Hydrogen atmospheres accreted early in a planet’s life are thought to escape due to the thermal energy of the particles, either slowly through Jean’s escape or quickly due to a thermal wind (Hunten, 1973). In both of these processes a significant portion of the particle velocity distribution exceeds escape velocity due to their thermal motion. This thermal energy can be gained from photo-heating (e.g. Yelle, 2004), giant impacts (Liu et al., 2015), or core-accretion (Ginzburg et al., 2016). Thermally-driven atmospheric escape is thought to lead to the gap between super-earths and mini-neptunes (Owen and Wu, 2017; Ginzburg et al., 2018; Biersteker and Schlichting, 2018).

Although thermal loss is an important process for dictating the evolution of the primary atmosphere or hydrogen nebula of a planet, secondary atmospheres created by planetary out-gassing are composed of gravitationally bound heavier species and have more difficulty escaping thermally after the primary atmosphere has been lost. Additional energy sources such as photo-chemical production (McElroy et al., 1977; Hunten, 1982), sputtering (Haff and Watson, 1979; Hunten, 1982), charge exchange (Chamberlain, 1977), or ion escape are necessary to drive significant escape of secondary atmospheres over the course of a planet’s lifetime (Hunten, 1982). Ion escape encompasses the collection of escape channels where escape occurs through loss of ionized species, including ion pickup (e.g. Luhmann and Kozyra, 1991), ion bulk escape (e.g. Brain et al., 2010), and the polar wind (e.g. Banks and Holzer, 1968; Yau et al., 2007). Both the energization of the particles through these processes and their eventual escape trajectories (or lack thereof) are critically affected by the dynamic plasma environment created by the interaction of the planet and the stellar wind.

The strength of a planet’s magnetic field is a critical factor in shaping a planetary plasma environment. In the solar system the only clearly habitable planet, Earth, has a strong intrinsic magnetic field driven by an internal dynamo, while Venus and Mars do not. This has led many to speculate that an intrinsic magnetic field is a critical component required to shield the planet from atmospheric erosion due to the solar wind (Driscoll and Bercovici, 2013; Lundin et al., 2007; do Nascimento et al., 2016; Driscoll, 2018), we call this the magnetic umbrella hypothesis. However, present measurements of escape rates suggest that all three planets are losing their atmosphere at similar rates (Strangeway et al., 2005), and the magnetic umbrella hypothesis has been called into question (Kasting, 2012; Brain et al., 2013). Understanding the complex interplay between the planetary plasma environment, solar wind driving, and ion loss is thus of critical importance.

Some investigation into the effects of an intrinsic magnetic field on ion loss has been made using global plasma models. Dong et al. (2017) compared escape rates from Proxima-b using a multi-species MHD model with a magnetized and an unmagnetized case and found larger escape rates in the unmagnetized case, while Garcia-Sage et al. (2017) found using a polar wind outflow model that the higher ionizing radiation levels of EUV radiation fluxed of M dwarfs compared to Earth may increase the effectiveness of the polar wind by orders of magnitude. Global hybrid modeling work from Kallio and Barabash (2012) looked at ion escape from a Mars type planet with a global dipole of 0, 10, 30, and 60 nT, and found maximum ion escape from the 10 nT case. Gunell et al. (2018) combined empirical measurements at Earth, Mars, and Venus with semi-analytic models to analyze the influence of planetary magnetic field over a large range of magnetic field values, finding multiple peaks in escape rate over the range that varied in strength and value by planet. Sakai et al. (2018) found an increased escape rate for a weakly magnetized Mars (100 nT) compared to an unmagnetized Mars in a multi-species MHD model. In combination these results suggest that there is much more to learn about about ion escape as it relates to planetary magnetic fields.

Weakly magnetized planets are an especially interesting regime as they represent the transition from unmagnetized planets with induced magnetospheres to magnetized planets with intrinsic magnetospheres. In the weak-field limit the scale size of the stellar wind interaction region is not dramatically changed. However, a global planetary magnetic field is typically approximated to the first order as a magnetic dipole centered at the planet, this introduces a new axis of symmetry. Furthermore, some studies indicate that planets around M-dwarfs may be weakly magnetized (Gaidos et al., 2010). Analysis of the Kepler data has shown that planetary systems are common around M-dwarfs (Kopparapu, 2013), and these systems also show the best promise of observing exoplanet atmospheres (Shields et al., 2016). Thus it is important to understand what effect weak magnetic fields will have on the subsequent ion escape.

Here we study the influence of weak planetary magnetic fields on the plasma environment of a planet and subsequent ion loss through a series of hybrid plasma simulations. Section 2 describes our methods, Section 4 details influence of the magnetic field on the planetary plasma environment, Section 5 describes the corresponding influence on ion escape, Section 6 contains a discussion of the results, and Section 7 summarizes the conclusions.

2 Methods

To analyze the effects of planetary magnetic field strength on ion escape we have run seven simulations of the solar wind interaction with increasingly magnetized planets. These simulations were performed using the hybrid plasma model RHybrid (Jarvinen et al., 2018). In a hybrid model the ions are treated as macroparticle clouds while the electrons are treated as a charge-neutralizing fluid. By treating ions as macroparticles ion kinetic effects associated with finite ion gyroradii are included, while treating the electrons as a fluid allows simulating a much larger volume than a corresponding fully kinetic approach.

Each ion macroparticle represents a group of ions that have the same velocity ($v_{i}$), central position ($x_{i}$), charge ($q_{i}$), and mass ($m_{i}$). Each cloud has the same shape and size as a grid cell and clouds move inside the grid cells obeying the Lorentz force such that

[TABLE]

where $\vec{E}$ and $\vec{B}$ are the electric and magnetic fields respectively. The electron charge density then follows from the quasi-neutral assumption with all ion species accumulated in the grid and total new charge density assumed to be zero in each cell.

Ampere’s law is used to calculate the current density from the magnetic field

[TABLE]

and then the electric field is found using Ohm’s law

[TABLE]

where the $\eta$ is the explicit resistivity profile, and $U_{e}$ is the electron fluid velocity. The value of $\eta$ was chosen to be $\eta=0.02\mu_{0}\Delta x^{2}/\Delta t$, such that the magnetic diffusion time scale becomes $\tau_{D}=\mu_{0}L_{B}^{2}/\eta=50\Delta t$, for the magnetic length scale $L_{B}\simeq\Delta x$. This choice adds diffusion in the propagation of the magnetic field (Ledvina et al., 2008), stabilizing the numerical integration, while ensuring that the magnetic field diffuses on timescales longer than the timestep $\Delta t$ . Note that the resistivity is not explicitly included in the Lorentz force. Finally, the magnetic field is advanced using Faraday’s Law

[TABLE]

Jarvinen et al. (2018) and references therein contain details of the numerical scheme. Highly parallelized RHybrid and its sequential predecessor code HYB have been used to study ion escape from Mars (Kallio and Janhunen, 2002; Kallio et al., 2010; Diéval et al., 2012), Venus (Jarvinen et al., 2009, 2013), Mecury (Kallio and Janhunen, 2003), Titan (Kallio et al., 2007), and exoplanets around M-dwarfs (Egan et al., 2019).

The results presented in this paper are in the Planet Stellar Orbital (PSO) coordinate system, such that the planet is at the origin, the x direction points from the planet against the undisturbed upstream stellar wind velocity vector, the z direction is perpendicular to the orbital plane of the planet, and y is the completion of a right-handed coordinate system. The physical distances are shown in planetary radii ($R_{p}$) that correspond to 3390 km ($=1$ Martian radius).

The externally forced stellar wind and interplanetary magnetic field (IMF) was chosen to be the same as in Egan et al. (2018) and is a case of a nominal solar wind interaction with Mars. The corresponding parameters are listed in Table 1. Note that the magnetic field is largely in the +y direction and thus the corresponding motional electric field is largely in the +z direction.

The ion production representative of the planetary ionosphere is implemented via Chapman profiles such that

[TABLE]

where $z^{\prime}$ is the normalized height parameter given by $z^{\prime}=(z-z_{0})/H$, $\chi$ is the solar zenith angle, $z_{0}$ is the reference height, and $H=16$ km is the scale height. Such a profile assumes an isothermal atmosphere that is ionized by plane-parallel, monochromatic radiation in the EUV. This is a simplified method to inject planetary ionospheric ions in the simulation using a pre-defined production rate and spatial emission profile near the exobase and the inner boundary of the model; limitations arising from such a choice are discussed in Section 6. The peak location of ion production occurs at an altitude of $z_{0}=300$ km, while the ion absorption lower boundary occurs at 150 km. The global emission rates $Q_{0}$ for each ion are picked to correspond to nominal values at Mars today ($Q_{0}(O_{2}^{+})=2\times 10^{25}$ s*-1*, $Q_{0}(O^{+})=1.4\times 10^{25}$ s*-1*) (Jarvinen et al., 2016, 2018; Egan et al., 2018, 2019).

The planetary magnetic field is enforced at the planet surface, and in the non-magnetized case is implemented as a super-conducting sphere. The range of planetary magnetic fields runs from unmagnetized to an equatorial surface dipole field of $B_{p}=150$ nT, which corresponds to approximately 1/400th of the Earth’s equatorial surface magnetic field. The magnetic dipole is aligned with the z axis such that equatorial dipole field point along $+z$.

Each simulation was run on a $240^{3}$ grid, with boundaries at $\pm 5R_{P}$ in the Y and Z directions and $[-6,+4]$ in the X direction. This leads to a corresponding resolution of $\Delta x=141$ km, which is comparable to an estimate of the ion gyroradius for $O^{+}$ ($\gamma=mv_{\perp}/qB\sim 390$ km). Each simulation was run for 100,000 time steps with $\Delta t=0.01s$ (satisfying the Courant condition), corresponding to roughly 8 stellar wind crossing times through the simulation domain.

3 Unmagnetized and Magnetized Planetary Paradigms

For unmagnetized planets with some atmosphere, the conductivity of the ionosphere induces a magnetosphere which acts to stand off the stellar wind in the absence of a planetary magnetic field. This creates a transition region that includes the Magnetic Pileup Boundary (MPB) due to the pileup of the IMF as it drapes around the planet and the Induced Magnetosphere Boundary (IMB) marked by a transition from stellar wind to planetary plasma (Nagy et al., 2004). The dynamics of ion escape are thus dominated by the plasma environment associated with the induced magnetosphere.

As the stellar wind directly interacts with the ionosphere at unmagnetized planet, the stellar wind motional electric field accelerates pickup ions into a plume (Cloutier et al., 1974; Luhmann and Kozyra, 1991; Dong et al., 2015). Additionally, cold ions can flow out on the night-side of the planet in tail escape (Fränz et al., 2015; Engwall et al., 2009).

For magnetized planets the intrinsic magnetic field acts to stand off the incoming stellar wind through a balance of the magnetic pressure and incoming stellar wind dynamic pressure. For a planetary magnetic field given by a dipole this pressure balance can be used to define the magnetic standoff distance Rs as follows

[TABLE]

where the additional factor of 2 in the numerator comes from the compression of the magnetic field at the magnetopause (Cravens, 2004). Because the stellar wind has less direct access to the ionosphere, the dynamics of the escaping ions are dictated by the plasma environment associated with the intrinsic magnetosphere and the topology of the field lines.

Magnetic topology is defined by the structure of magnetic field lines and their relationship to the surrounding environment. Open magnetic field lines have one end connected to the planet with the other connected directly to the stellar wind, while closed field lines have both foot-points connected to the planet. As a planet’s magnetic field increases, the area of the ionosphere that is magnetically connected to the stellar wind via open field lines can be expected to decrease.

The locations of closed and open field lines can be estimated analytically for magnetized planets using the magnetic stand-off distance. Because the magnetic standoff distance is associated with the last closed field line, the latitude of the polar cap, or region of open field lines is thus given by

[TABLE]

and the corresponding polar cap solid angle is then

[TABLE]

When the standoff distance becomes less than the planet radius the polar cap angle becomes saturated at $4\pi$, or twice the angle of an entire hemisphere.

As ions are constrained to gyrate around magnetic field lines, magnetic topology has key implications for determining ion escape rates. At strongly magnetized planets ions escape along open field lines at the poles in a process known as cusp escape (Lockwood et al., 1985). If the ions are accelerated due to the electron pressure gradients along the poles it is called the polar wind (Banks and Holzer, 1968; Yau et al., 2007).

The transition between magnetized and unmagnetized systems is a transition between these two paradigms; unmagnetized planets are dominated by induced magnetosphere dynamics and direct acceleration of planetary ions via the stellar wind motional electric field, while magnetized planets are dominated by their intrinsic magnetospheres and global magnetic topology. The following results sections explore this transition in plasma environments and the corresponding changes ion escape.

4 Plasma Environment

In this section we analyze the structure of the plasma environment and magnetic topology near the planet in our set of runs with an increasing planetary magnetic field strength. Even a moderate change in environment and topology associated with a weakly magnetized planet can affect the ion dynamics near the planet and thus the ion outflow.

4.1 Tail Twisting

The tail magnetic morphology is critically affected by both the IMF and the dipole structure. In an unmagnetized planet the stellar wind magnetic field drapes around the planet in a roughly symmetric fashion, creating a central current sheet within the induced magnetosphere. Observations at Mars have shown that the current sheet is very dynamic with local variations ascribed to both twisting and flapping motions (DiBraccio et al., 2017, 2018; Chai et al., 2019), which may be associated with changes in response to upstream IMF or stellar wind pressure, kink-like flapping due to waves propagating along the current sheet, or reconnection of crustal fields (Luhmann et al., 2015). For magnetized planets the intrinsic dipole field also creates a current sheet in the tail oriented perpendicularly to the dipole axis. This current sheet also experiences flapping motion (Speiser and Ness, 1967) and is strongly influenced by reconnection events (Nakamura et al., 2006).

In this set of simulations the dipole axis and the IMF are roughly perpendicular, making the topology of the current sheet sensitive to the relative strengths of these fields. Figure 1 shows slices through the tail at $x=-2$ for each simulation for both $O_{2}^{+}$ number density and the x-component of the magnetic field, which illustrate the location of the current sheet. As the field strength increases the tail becomes more twisted such that the central area is dominated by the dipole contribution while farther out is dominated by the IMF draping. This structure shows the transition in relative influences of the induced and intrinsic magnetospheres, as well as the effect on ion escape morphology.

4.2 Magnetic Topology Mapping

Figure 2 shows magnetic field line tracing through our models with orange lines indicating closed field lines and blue indicating open field lines. The lines were traced in both directions from a spherical shell at 400 km altitude. These field line tracings reflect the above intuition where stronger magnetic fields create larger regions with more closed magnetic field lines.

Because the magnetization of the planet is weak, the interplanetary magnetic field (IMF) carried by the stellar wind has a large impact on the structure of the magnetosphere. This is similar to what has been found at Mercury (Slavin, 2004; Luhmann et al., 1998; Kabin et al., 2000) and Ganymede (Kivelson et al., 1997) where the external magnetic field induces polar alignment. Figure 2 shows this influence where the $+y$ oriented IMF causes the open field lines at the poles to twist towards the dawn ($+y$) and dusk ($-y$) sides in the northern and southern hemispheres respectively, such that at larger distances the field lines align with the IMF.

The closed field lines attached to the planet on both sides also show twisting to align with the IMF. This causes the closed field lines to attach asymmetrically in longitude; a field line attached at a dawnward longitude in the northern hemisphere will close on the duskward longitude in the southern hemisphere. As the magnetic field strength increases, this effect decreases due to the larger strength of the planetary field in comparison to the IMF.

The lower portion of Figure 2 shows a latitude-longitude map of the field line topology at the same altitude where the field line tracings began, with dashed lines indicating the critical latitude $\theta_{c}$ given by Equation 7. While field lines outside $\pm\theta_{c}$ are indicated as open within the analytic formulation described earlier, the field line mapping shows significant area of closed field lines outside this region, particularly for the lowest magnetic field strength (10 nT) and on the night side (longitude $\sim 180$). This occurs because of the aforementioned asymmetries, as well as the neglect of the induced magnetospheric pressure in the analytic formulation.

4.3 Magnetic Standoff

In the upper panel of Figure 3 the magnetic topology is compared across magnetic field strength by assessing the area of the polar cap open to the stellar wind, an estimate of the region from which ions may escape. The “full field line trace” is calculated by summing the solid angle of all open field lines identified in the tracing associated with Figure 2. The “axially symmetric field trace” is calculated finding the region of closed field lines along the noon longitude and assuming a polar axis symmetry. The “analytic cap area” is calculated assuming a polar axially symmetric system with a critical latitude given by Equation 7. The discrepancy between the two calculated values is due to the lack of symmetry in the system because of day-night variation, and symmetry breaking induced by the IMF as discussed in the following subsection. Both values are useful to show due the physical accuracy (topology area) and directness of the comparison with the analytic value (symmetric area).

In general the analytic formulation overestimates the polar cap solid angle, especially at small magnetic field values. This occurs because the analytic estimate must be completely open for all field values where the magnetic field pressure is less than the stellar wind pressure. However, this once again neglects the influence of the induced magnetosphere.

We also show a best fit scaling relationship between the polar cap angle and magnetic field for the full field line trace, where $\Omega_{c}\propto{B_{p}}^{-0.34}$. A scaling law exponent of $-1/3$ can be recovered by expanding Equation 8 around small values of $R_{p}/R_{s}$, and then substituting Equation 6 for $R_{s}$. This scaling matches well even to small magnetic field values ($B_{P}\sim 10$ nT).

The bottom panel of Figure 3 shows the transition from induced to intrinsic magnetospheres through a variety of metrics designed to assess standoff of the stellar wind by the planet. The analytic stand-off altitude (Equation 6) is calculated assuming pressure balance between the incoming stellar wind and the intrinsic magnetosphere of the planet. This is most clearly related to the pressure equilibrium point (the altitude along the sub-stellar line at which the magnetic pressure is equivalent to the stellar wind pressure); however, the pressure equilibrium point occurs at much higher altitude in the magnetosheath due to the expected thickness of the transition region, corresponding to a few proton gyroradii ($\gamma_{H^{+}}\sim 730$ km) (Cravens, 2004).

The analytic stand-off altitude is also closely related to the magnetic topology, as shown in Equation 7. Here, rather than calculating the latitude associated with the last closed field line, we show the altitude along the sub-stellar line of the last closed field line. This matches fairly well with the analytic stand-off distance but occurs slightly higher due to the the same asymmetry and induced field pressure support as described in Section 4.2.

Finally, we show a measure of the magnetospheric boundary through identification of the steepest gradient in electron number density, a method used to identify the induced magnetospheric boundary at unmagnetized planets (Vogt et al., 2015). At unmagnetized planets this can be thought of as the “top” of the ionosphere, where there is a transition from planetary to stellar wind plasma. At magnetized planets this becomes the outer region of the plasmasphere (see Section 5.3), where closed magnetic field lines can confine plasma. This measurement thus has a clear relationship to the closed field line boundary, while also having a well defined meaning for unmagnetized planets.

In Figure 3 we see that the last closed field line and electron density gradient altitudes agree fairly well with each other and the analytic standoff distance at higher magnetic field strengths. At lower magnetic field strengths, the electron density gradient flattens out due to the intrinsic gradient associated with the ionosphere. This transition occurs where the magnetic stand-off crosses the unmagnetized induced magnetosphere boundary (847 km) (Trotignon et al., 2006). Critically, this point is also identical to the magnetic field associated with maximum ion outflow as shown in Section 5.

This transition between induced and intrinsic magnetospheres is shown schematically in Figure 4. The induced magnetosphere boundary remains constant with magnetic field as it is determined by the ionospheric conductivity and pressure support against the stellar wind. For the weakest magnetic fields ($B_{P}<50$ nT) the induced magnetosphere is the dominant source of pressure support against the stellar wind, the last closed field lines occur below the boundary, and the electron density gradient is then associated with the induced boundary. At stronger magnetic fields ($B_{P}>50$ nT) the intrinsic magnetosphere is the dominant source of pressure support against the stellar wind, the last closed field lines occur above the induced magnetosphere boundary, and the electron pressure gradient is associated with trapped plasma within the closed magnetic field line region. The transition point occurs when the intrinsic magnetosphere standoff distance reaches the induced magnetosphere models; for this set of simulations this occurs at $B_{P}=50$ nT. In the following section we will discuss the implications this pattern has for ion escape.

5 Ion Escape

As discussed in Section 4, the transition between unmagnetized and magnetized planets affects the stand-off of the stellar wind and the opening angle of the polar cap, which in turn affect the pickup ion and cusp escape processes. Here we analyze the trends in ion escape rates and energy and connect them to the influence of the surrounding plasma environment.

5.1 Escape Rates and Energies

Figure 5 shows the escape rate , escape power, and average ion escape energy (escape rate divided by escape power) for each simulation run. These parameters are calculated as

[TABLE]

[TABLE]

[TABLE]

respectively, where $n_{i}$ is the ion number density, $u_{i}$ is the ion velocity, $m_{i}$ is the ion mass, and $\mu$ is the ion mass expressed in units of the proton mass. We calculate these rates by integrating over a sphere far from the planet such that the sphere is large enough ($r=4.0R_{P}$) that the flux value is insensitive to changes in radius, and all outbound flux escapes and does not return to the planet, even for our most strongly magnetized simulation cases.

We find that the escape rate increases for magnetic fields up to $B_{P}=50$ nT and then begins to decrease again, with a maximum variation of a factor of $\sim 3$. The magnetic field associated with peak escape rate ($B_{max}$) is equivalent to the magnetic field where the induced magnetosphere altitude equals the intrinsic magnetosphere altitude. This divides the ion escape trends into different regimes illustrated in Figure 6; weaker escape associated with predominately induced magnetospheres, stronger escape associated with magnetic shielding and equivalent induced and intrinsic magnetospheres, and weaker escape associated with strong magnetospheres and plasmasphere trapping. Each of these regimes will be explored more fully in the following subsections. Additionally, we show a power law fit for simulations $B_{P}\geq 50$ nT where we find that the escape rate scales as ${B_{P}}^{-0.67}$.

Escape power shows a steady decrease with magnetic field, and a power law scaling of ${B_{P}}^{-0.60}$ for $B_{P}\geq 50$ nT. Average escape energy is calculated by dividing integrated escape power over the spherical shell by integrated escape flux. For both $O_{2}^{+}$ and $O^{+}$ Figure 5 shows the reduced escape energy is highest for small magnetic field values before reaching an asymptote at around 50 nT.

Both power laws calculated for escape rate and escape power have clear lower limits at $B_{p}=50~{}nT$, a point we define as $B_{max}$. These power laws also likely have upper limits that may prevent their extension to very strongly magnetized systems like Earth, as indicated by the various shaded regions in Figure 5. This is discussed further in Section 6.4.

Figure 7 shows a histogram of the full velocity distribution of escaping particles through the $R=4.0\;R_{P}$ sphere, weighted by $\textrm{area}\times n_{i}\times v_{i}$ and normalized such that the integral over the distribution is 1. Dashed lines indicate the average escape energies. This figure shows that the decrease in escape energy with increasing magnetic field is caused by a corresponding decrease in high energy particles in the tail of the velocity distributions. As the magnetic field begins to stand off the stellar wind, fewer ions are exposed to the strong motional electric field, and thus fewer ions are accelerated to high velocities.

5.2 Southern Hemisphere Shielding

The increase in escape rate with increasing magnetic field is caused by stellar wind standoff in the southern hemisphere decreasing the flux of pickup ions accelerated towards the planet and preciptating in the inner boundary. While the motional electric field in the northern hemisphere causes ion pickup and the plume feature, in the southern hemisphere the motional electric field accelerates the ions into the planet, preventing their escape. This effect is illustrated in Figure 6, with the direction of the motional electric field indicated by the green arrow. While particles injected in the southern hemisphere of the 10 nT simulation are accelerated towards the planet and do not escape, particles injected in the same location in the 50 nT simulation do escape due to the increased stand-off of the stellar wind.

Figure 8 shows slices of the total motional electric field, total electron velocity, and ratio of the fluid velocity and the ion velocity for the 10, 50, and 100 nT simulations. In the 10 nT case strong electric fields due to high electron velocities permeate the region of planetary ions. In the 50 and 100 nT cases the increasing magnetic pressure slows the stellar wind down much more before reaching the planet, decreasing the magnitude of the motional electric field. The ions in these more strongly magnetized cases can then flow farther from the planet along open field lines before being exposed to the strong electric field, and are thus accelerated down-tail rather than into the planet. This effect will saturate roughly when the standoff distance reaches the altitude of the IMB at the terminator, which is illustrated occurring at 50 nT in Figure 3.

5.3 Plasmasphere Trapping

The decrease in escape rate at higher magnetic field strengths is due to plasmasphere trapping, where ions become trapped within the closed field line regions. This is depicted in Figure 6 where particles escape in the 50 nT simulation, but are confined within the closed magnetic field line in the 100 nT simulation. This effect is related to the decrease in polar cap solid angle illustrated in the lower panel of Figure 3. As more area of the planet is wrapped in stronger closed magnetic field lines, ions will have greater difficulty escaping these closed field regions.

As the magnetic field of the planet increases creating a region of closed magnetic field lines, the planet develops a region of plasma protected from stellar wind. This toroidal region of dense cold plasma connects to the top of the ionosphere and is known at Earth as the plasmasphere. The development of the plasmasphere is shown in Figure 9, where higher densities of planetary ions occur within the planetary dipole field where the motional electric field becomes small. Because the magnetic field stands off the stellar wind, the motional electric field becomes very small within the the plasmasphere. The ion motions are then dominated by gyromotion along the magnetic field lines, and thus the ions are trapped within the plasmasphere.

Figure 5 shows a best fit scaling relation for the escape rate of the plasmasphere trapping dominated simulations ($\geq 50$ nT), where $\Phi\propto{B_{P}}^{-2/3}$ or $\Phi\propto\Omega_{c}^{2}$. This is a steeper relationship than typically expected where $\Phi\propto\Omega_{c}$ (e.g. Gunell et al., 2018). This is likely due to the production rate variation with stellar zenith angle. Because more ions are produced in the equatorial region, the relationship of escape rate with polar cap angle will be steeper than $\Phi\propto\Omega_{c}$. Thus, the escape rate relationship with magnetic field will be sensitive to the ion production in relation to the magnetic field topology. This is further discussed in Section 6.

5.4 Mass Flux and Power Coupling

Though intrinsic magnetic fields were long thought to shield atmospheres from being stripped by interaction with a stellar wind, it has also been pointed out recently that intrinsic fields increase the cross-sectional area through which a planet interacts with the wind, and therefore the amount of energy or momentum it could ‘collect’ from the wind to drive escape (Brain et al., 2013). A common technique used to discuss the relationship between input and output from a given reservoir is through the idea of a coupling efficiency. This idea has also been used in the context of planetary atmospheres where the mass loss rate ($\dot{M}_{out}$) or the escaping power ($P_{out}$) are coupled to the mass flux ($\dot{M}_{in}$) or power input ($P_{in}$) from the stellar wind with some efficiency $\epsilon$ (e.g. Strangeway et al., 2005; Dubinin et al., 2006; Ramstad et al., 2017b; Egan et al., 2019). The key questions then become how to define and parameterize the inflow flux/power and identifying the typical efficiency of various processes.

We put this framework in a more general form by including a correlation constant $k$ such that $\dot{M}_{out}=\epsilon(\dot{M}_{in})^{k}$ (or $P_{out}=\epsilon(P_{in})^{k}$), where $k=1$ indicates a traditional linear coupling with a constant efficiency.

We define the inflow properties by scaling the stellar wind parameters with an interaction area defined by the cross-sectional area of the bow shock at the terminator ($x=0$) as

[TABLE]

and

[TABLE]

where $\rho_{sw}$ and $u_{sw}$ are the density and velocity of the undisturbed stellar wind. The area A is calculated by identifying the shock boundary in the $x=0$ plane and summing the corresponding area in the $x=0$ plane within the boundary.. We identify the boundary in this way rather than an area calculated using a symmetric shock radius $R_{shock}$ such that $A=\pi R_{shock}^{2}$ to account for bow shock asymmetry. As the magnetic field increases the effective obstacle size presented by the planet also increases; this then increases the cross sectional area of the bow shock. The outflow properties are determined in the same method described in earlier in Section 5.1.

The resulting inflow and outflow fluxes and power for each simulation are shown in Figure 10, as well as the cross sectional areas calculated for each simulation. We find that the cross sectional area increases with magnetic field as roughly $A\propto{B_{P}}^{0.41}$ or $A\propto R_{S}^{1.23}$ (by Equation 6).

Neither the mass flux nor the power are coupled with a linear efficiency, nor does either relationship show a constant power law correlation over the entire range. This is indicative of the same change in ion escape processes as described in the preceding sections on escape trends, despite a constant scaling in cross-sectional area. For the higher magnetic field values where plasmasphere trapping is the dominant process affecting escape, the mass flux and power correlations have indexes of $k=-1.60$ ($\epsilon=-0.48$) and $k=-1.43$ ($\epsilon=25.8$) respectively. As the inflow mass (power) scales with area ($A\propto{B_{P}}^{0.41}$) and the outflow mass (power) scales as $\propto{B_{P}}^{-0.66}$ ($\propto{B_{P}}^{-0.6}$), both correlation constants follow straightforwardly.

6 Discussion

6.1 Stellar Driving Dependence

All the simulations shown in this paper were run with the same nominal Martian solar wind driving; however, stellar winds are dynamic and can dramatically change ion loss (Ramstad et al., 2017b; Brain et al., 2017). Here we discuss the impact of our choice in stellar wind on our results.

As shown in Equation 6, the stellar wind pressure controls the standoff distance associated with an intrinsic magnetosphere. The induced magnetosphere boundary has also been shown to depend weakly on dynamic pressure and extreme ultra-violet input, but much less so than the intrinsic magnetosphere (Ramstad et al., 2017a). As the magnetic field associated with the peak escape rate, $B_{max}$, occurs when the intrinsic magnetosphere reaches the induced magnetosphere, it will also scale with stellar wind pressure as $B_{max}\propto P_{sw}^{1/2}$. Therefore a planet experiencing a stellar wind pressure increase by orders of magnitude, as may be likely for a planet in the habitable zone around an M-dwarf (Vidotto et al., 2015; Garraffo et al., 2017), will have a much larger $B_{max}$.

As shown in Section 4, due to the relative strengths of the planetary dipole and the IMF, the IMF critically affects the magnetic topology of the system. The IMF can be described by a magnitude, the clock angle $\phi$ (angle in the YZ plane, perpendicular to the sub-stellar line), and the cone angle $\theta$ (angle between the IMF and the sub-stellar line). Our study has assumed an IMF that is perpendicular to the dipole axis, thus a motional electric field that is aligned with the dipole axis. In the opposite case where the IMF is aligned with the dipole axis, the southern hemisphere shielding effect would be mitigated as the motional electric field will stand-off or accelerate into the plasmasphere rather than the poles. Thus, there is a critical dependence on the IMF clock angle relative to the dipole axis, necessitating further study for different configurations.

Similarly, the angle between the dipole magnetic field axis and the sub-stellar line also impacts the escape rate relations we have shown here. The ion production method we have chosen is dependent on the stellar zenith angle, such that the highest production occurs at the sub-stellar point. A dramatically different configuration would occur for a planetary pole aligned with this peak in ion production; however, this is a less likely scenario due to spin-orbit alignment. Further, planets orbiting close to their stars have fast orbital velocities. In the case when the orbital velocity is the same as the stellar wind speed, the upstream stellar wind vector and the sub-stellar line have an angle of 45∘.

6.2 Model Limitations

In a global hybrid model like RHybrid, specific escape rates depend on the inner boundary ionospheric emission rates. The emission rate of inospheric ions is a free parameter, which is fixed by comparing a model to in situ observations from solar system planets (e.g. Jarvinen et al., 2009, 2018). While we show escape rates here that are comparable to current escape rates at Mars (Lundin et al., 1990; Jakosky, 2017; Brain et al., 2015), we caution against the over-interpretation of these rates due to their dependence on the lower boundary, and instead concentrate on the relative differences.

Furthermore, our implementation of ionospheric emission is driven by a predefined Chapman ion production profile. Currently, accurately resolving ionospheric dynamics on the same domain as the global magnetosphere is computationally infeasible. Although some simulations suites include a one way coupling of an ionosphere implementation to a global magnetosphere (e.g. Glocer et al., 2009; Brecht et al., 2016; Modolo et al., 2016), we chose to restrict our ionospheric production to a simple approximation. This allows us to focus on the magnetospheric physics associated with an increasing magnetic field, without introducing an additional layer of uncertainty from coupling an additional model.

We have also neglected to include an ionospheric electrodynamics model, and instead use a constant resistivity value above the lower boundary and zero resistivity at the lower boundary. Global thermospheric structure (Roble et al., 1982) and ionospheric-magnetospheric coupling (Ridley et al., 2004) are both affected by ionospheric electrodynamics through mechanisms such as current closure, atmospheric Joule heating, and Alfven wave dissipation.

Impact ionization, neutral structure, and resistivity are all affected by auroral precipitation (Cravens et al., 1995; Deehr et al., 1980; Tinsley, 1979; Robinson et al., 1987; Fuller-Rowell and Evans, 1987), which will change with magnetic field strength. In the weakly magnetized cases we have chosen here the polar cap solid angle (and the corresponding precipitation per area) varies by a factor of a few across the simulations. While this is likely to become important at large magnetic field strengths (see Section 6.2), it is likely a secondary effect behind the basic ion production method we use here.

6.3 Comparison With Existing Results

In general our conclusions agree with similar results from previous analytic and numerical studies, but add a deeper understanding of the physics controlling the transition from unmagnetized to weakly magnetized planets.

An analytic study by Gunell et al. (2018) showed a similar peak in ion escape rates at the point where the standoff radius reaches the IMB; however, their model enforces a sharp transition in escape processes at this point rather than a gradual build up. This study also uses the typical assumption of linear dependence of ion escape with the polar cap angle, whereas we find a stronger dependence.

Kallio and Barabash (2012) used a hybrid model to study weakly magnetized planets and found a peak ion escape rate for a Mars-like planet with a 10 nT magnetic field. Our peak escape rate occurs at a higher point due to the stronger stellar wind conditions. Furthermore, the same group studied the magnetic morphology of the same system and found similar trends in tail twisting and IMF control; however, the effects were weaker due to the relative weakness of the IMF in comparison with the dipole strengths for a similar stellar wind pressure (Kallio et al., 2008).

Blackman and Tarduno (2018) studied the impact of magnetic fields on stellar wind erosion of planetary atmospheres using an analytic estimation of the maximum escape rate for a given system. This model assumes there is a direct relationship between the inbound and outbound mass flux, an assumption our results contradict in the weakly magnetized regime; however, they largely concentrate on magnetic fields estimated for ancient Earth, which are larger than those studied here.

Other studies examining the escape rates of exoplanets with and without magnetic fields have been performed (e.g. Dong et al., 2017); however, these results are not directly comparable as the magnetic fields studied are much stronger than the range considered here.

6.4 Limits of Derived Scaling Laws

In Table 2 we summarize the scaling laws we have fit for various parameters from the simulations with $B_{P}\geq 50$ nT, along with estimations for the lower and upper limits of the power law applicability. Although we can calculate error bars on the scaling relations from the model results, these errors do not reflect the sensitivity of the calculated power law to the model assumptions. Instead, here we discuss some of the main parameters and assumptions these laws are most sensitive to.

As shown in Section 4.3, the polar cap solid angle should scale as ${B_{P}}^{-1/3}$ for large magnetic field values. This matches our calculated scaling law very well, even to small magnetic field values leading us to place an approximate lower limit of $R_{S}\sim R_{P}$. While magnetic topology is a fairly straightforward feature to interpret from our models, it is also sensitive to the magnetic dipole orientation with respect to the IMF and sub-stellar line. This relationship does not have a clear upper limit for reasonable values of planetary magnetic field; however, as the polar cap region becomes smaller care must be taken to understand the asymetries and dependence on IMF (e.g. Kabin et al., 2004; Merkin and Goodrich, 2007)

We find a mass loss rate that scales more strongly with magnetic field than predicted, likely due to the variation in ion production with SZA. This result is quite sensitive to the ion production mechanism, as well as the magnetic topology. From the simulations here we have shown the lower limit will be $R_{S}\simeq R_{IMB}$ or equivalently $B_{P}\simeq B_{max}$.

The upper limit is more difficult to estimate. First, the power law will become more shallow with decreasing polar cap angle, likely approaching the predicted $\dot{M}\propto\Omega\propto B^{-1/3}$ relationship. Additionally, as the polar cap angle becomes small, cusp dynamics associated with ionospheric electrodynamics and energization from particle precipitation will become more dominant. Strangeway et al. (2005) showed that cusp ionospheric outflow at the Earth is associated with polar cap Poynting flux and electron precipitation leading to Joule dissipation and electron heating, causing increased ion upwelling and subsequent outflow. Both of these effects are associated with a large power input to a small polar cap area coupling to the ionosphere.

A polar cap solid angle an order of magnitude smaller than our base case will occur at $\Omega_{c}\sim 0.1\times 4\pi$, or $P_{B}\sim(5.26)^{6}P_{SW}$. For this set of simulations this would occur for magnetic fields of $B_{P}\sim 3500$ nT, which is still much smaller than the Earth’s surface field, where cusp escape processes are known to be relevant (Li et al., 2012). While the escaping ion flux may in fact continue to decrease past this point, the dominant physical mechanisms are likely to be quite different due to the intensity of precipitating particle flux in the small auroral region and the scale of the electrodynamic systems.

As shown in Figure 5, the average escape energy per particle becomes constant with increasing magnetic field and the corresponding total power of ion outflow decreases as $P_{esc}\propto B^{0.6}$. Thus, if the average energy of escaping particles remains constant, the $P_{esc}$ relationship is determined by the escape flux relationship, and will have the same upper and lower limits. As the escape energy is determined by the stellar wind velocity and corresponding motional electric field, this is likely to be true as long as the wind maintains its momentum through an encounter with the planetary obstacle. This will be true in the absence of significant mass-loading.

The interaction area is found to scale as $A\propto B^{0.41}$, which is shallower than one might expect for an interaction area that scales with a length scale ($R_{S}$) squared, or $A\propto B^{2/3}$. While we recover the standoff radius relationship as expected (see Figure 3), the interaction area is more sensitive to the stellar wind parameters and corresponding shock properties. Furthermore, the shape of the shock will be affected by the shape of the obstacle; the assumption of a spherically symmetric obstacle is particularly poor for small magnetic field values (as shown in Section 4).

7 Conclusions

Understanding the influence of planetary magnetic field on ion escape is a critical area for assessing the potential habitability of exoplanets. Here we have presented a systematic study of the influence of a weak magnetic dipole on the plasma environment and corresponding ion escape from a Mars-like exoplanet using a hybrid plasma model.

We have found that an intrinsic magnetic field does not always act to decrease the ion escape rate, and instead leads to an increase in escape rate up until the point where the magnetic standoff reaches the terminator IMB before then beginning to decrease again. This reflects a balance between southern hemisphere shielding of ions from the motional electric field and increased plasmasphere trapping due to smaller polar cap area. The value of the magnetic field associated with maximum ion escape is therefore dependent on the stellar wind pressure.

Where possible we have calculated scaling laws for the relationship of various fundamental quantities and magnetic field. The scaling laws are interpreted in the global context of planetary atmospheric escape theory. Lower limits for these scaling laws have been derived from the simulation results, while upper limits are postulated via predicted limits on the observed processes.

Going forward, more study is needed to assess effects associated with the dynamics of the stellar wind for planets with weak magnetospheres due to the extreme influence of the stellar wind on these systems, as well as a larger array of IMF configurations relative to the planetary dipole. Furthermore, the effects of stronger planetary magnetic fields must be assessed using a simulation platform that fully couples a polar wind model.

8 Acknowledgments

Hilary Egan was supported by the Department of Energy Computational Science Graduate Fellowship, Maven , and a generous donation from Marcia Grand. This research used resources of the National Energy Research Scientific Computing Center (NERSC), a U.S. Department of Energy Office of Science User Facility operated under Contract No. DE-AC02-05CH11231. Global hybrid simulations were performed using the RHybrid simulation platform, which is available under an open source license by the Finnish Meteorological Institute (https://github.com/fmihpc/rhybrid).