Estimating ion escape from unmagnetized planets

We propose a new method to estimate ion escape from unmagnetized planets that combines observations and models. Assuming that upstream solar wind conditions are known, a computer model of the interaction between the solar wind and the planet is executed for different ionospheric ion production rates. This results in different amounts of mass loading of the solar wind. We then obtain the ion escape rate from the model run that best fits observations of the bow shock location. As an example of the method, we estimate the heavy-ion escape from Mars on 1 March 2015 to be 2×10²⁴ ions s⁻¹, using a hybrid plasma model and observations by the Mars Atmosphere and Volatile Evolution (MAVEN) and Mars Express (MEX) missions. This method enables studies on how escape depends on different parameters as well as studies on escape rates during extreme solar wind conditions; moreover, the technique is applicable to studies of escape in the early solar system and at exoplanets.

1 Introduction

Ion escape to space is important for the evolution of planetary atmospheres. Neutral atoms or molecules in the upper parts of the atmosphere can be ionized by factors such as ultraviolet (UV) photons, charge exchange, and electron impacts. The newly created ion can then be energized by electric fields and transported away by the stellar wind, overcoming gravity, resulting in atmospheric loss.

For planets in our solar system, we can observe the present-day escape of planetary ions by directly observing the ion flux near a planet. This is done using an ion detector on a spacecraft and gives us the flux of ions along the trajectory of the spacecraft. As the flux of escaping ions is highly variable both temporally and spatially, accurately estimating the escape of ions can require observations over many years to get an average escape rate. Investigating how the escape rate of ions depends on different parameters (e.g., upstream solar wind conditions) is even more difficult due to the large amounts of observations needed to get sufficient statistics.

Another way of estimating ion escape rates is to use computer models of the solar wind interaction with a planet. An advantage of models compared with observations is that the full three-dimensional escape is obtained at every instance. In addition, there are no limitations with respect to sensitivity, energy range, or field of view, in contrast to the limitations on observations. However, accurately estimating ion escape requires that the model contains all of the important physical processes in sufficient detail. It is questionable if this is the case at present.

Here, we propose an alternative method of estimating ion escape that uses both observations and computer models. Instead of directly observing the escaping ions, we use observations of other plasma quantities near the planet. We then employ a parameterized model and find the set of model parameters that gives the best fit between the model and the observations. The model escape rate for the best fit parameters gives us an estimate of the ion escape rate.

This approach allows us to use data sets traditionally not used for ion escape estimates, such as magnetic field and electron observations. We can also estimate the escape rate from a very small set of observations, during one orbit of a spacecraft around a planet or during one flyby of a planet.

To illustrate this general method, we estimate the escape rate of ions from Mars during one bow shock crossing of the Mars Atmosphere and Volatile Evolution (MAVEN) spacecraft (Jakosky et al., 2015) using magnetic field observations of the bow shock location, observations of the upstream solar wind conditions, and a hybrid plasma model. We also use Mars Express (MEX) (Barabash et al., 2007) observations of electrons to verify our findings.

The location and shape of the Martian bow shock have been the topic of many observational studies that have aimed to establish where the bow shock is located and its controlling parameters. Recently Hall et al. (2016) found a seasonal dependence of the bow shock location and, later, also a solar cycle dependence (Hall et al., 2019). The seasonal dependence should be due to the changing distance of Mars from the Sun; however, it is not straight forward to deduce if it is due to changing solar wind pressure or changing UV insolation, as both scale in the same way with distance from the Sun. Moreover, both the solar wind and the UV insolation change over a solar cycle. Regarding the effect of crustal magnetic fields on the bow shock location, Gruesbeck et al. (2018) found such a dependence; furthermore, a modeling study by Fang et al. (2017) found that a large part of the variability in escape may be due to the crustal fields, reducing and enhancing escape, depending on the location. Regarding the shape of the bow shock, Vignes et al. (2002) found that it is furthest away from the planet in the hemisphere in the direction of the solar wind convective electric field.

However, it has been noted that, apart from the upstream solar wind conditions, the factor controlling the location of the shock is the amount of mass loading of the solar wind by ionospheric ions. This was noted for Venus by Alexander and Russell (1985) and later for Mars by Vignes et al. (2002) and Mazelle et al. (2004). The UV flux, atmospheric state, ionospheric chemistry, magnetic anomalies, and similar factors will all affect the location of the bow shock, although only indirectly through the amount of mass loading. Therefore, given upstream solar wind conditions, we should be able to use the amount of mass loading as a free parameter when modeling the location of the bow shock.

2 Method

The algorithm presented in this paper is general and can be applied to any model of the interaction between Mars and the solar wind. Here, to illustrate the method, we use a very simple hybrid model.

We now describe the hybrid plasma solver used; the adaptation for Mars; the parameters used; the observations of magnetic field, ions, and electrons used; and, finally, the algorithm to estimate ionospheric ion escape.

2.1 Hybrid model

In the hybrid approximation, ions are treated as particles, and electrons are treated as a massless fluid. The trajectories of the ions are computed from the Lorentz force, given the electric and the magnetic fields. The electric field is

$$\begin{array}{}\text{(1)}& E={\displaystyle \frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{I}}}}\left(-J_{\mathrm{I}}\times B+J\times B-\mathrm{\nabla }{p}_{\mathrm{e}}\right)+\mathit{\eta }J,\end{array}$$

where ρ_I is the ion charge density, J_I is the ion current density, p_e is the electron pressure, and η is the resistivity. The current is computed from $J={\mathit{\mu }}_{\mathrm{0}}^{-\mathrm{1}}\mathrm{\nabla }\times B$, where ${\mathit{\mu }}_{\mathrm{0}}=\mathrm{4}\mathit{\pi }\cdot {\mathrm{10}}^{-\mathrm{7}}$ is the magnetic constant.

Faraday's law is then used to advance the magnetic field in time:

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

Further details on the hybrid model used here, the discretization, and the handling of vacuum regions can be found in Holmstrom et al. (2012).

2.2 Mars model

In the hybrid simulation domain, Mars is modeled as a resistive sphere, of radius R, centered at the origin, where all ions that hit the obstacle are removed from the simulation.

The ionosphere is represented by the production of a single species of ions according to an analytical Chapman ionospheric profile (Holmstrom and Wang, 2015), where the production rate of ions is given by

$$\begin{array}{}\text{(2)}& \begin{array}{rl}P(h,\mathit{\chi })& =p{e}^{\mathrm{1}-y-\sec \mathit{\chi }{e}^{-y}},y=\left(h-h_{\mathrm{0}}\right)/H,\\ h& \ge \mathrm{0},\mathrm{0}\le \mathit{\chi }<\mathit{\pi }/\mathrm{2}.\end{array}\end{array}$$

Here, h [m] is the height above the planet surface; χ [rad] is the solar zenith angle (SZA); p [m⁻³ s⁻¹] is the maximum production along the sub-solar line, at height h₀ [m]; and H [m] is the scale height. Ions are then randomly placed according to this production function and are also given a random thermal velocity corresponding to a temperature of 200 K.

We note that the ion production is a free parameter in such a model. This is in contrast to models that self-consistently include ionospheric chemistry and a neutral atmosphere (Brecht et al., 2016). Usually this free parameter is seen as a limitation of the model; however, here, we use this as an advantage in order to find a best fit to observations. This means that the exact processes in the ionosphere that produce the ions or how they are transported to the top of the ionosphere is not important.

Regarding the composition of escaping ionospheric ions, a study by Carlsson et al. (2006) based on MEX observations found flux ratios of O${}_{\mathrm{2}}^{+}$ $/$ O${}^{+}=\mathrm{0.9}$ and CO${}_{\mathrm{2}}^{+}$ $/$ O${}^{+}=\mathrm{0.2}$. Later, Rojas-Castillo et al. (2019) estimated that O${}_{\mathrm{2}}^{+}$ $/$ O${}^{+}=\mathrm{0.76}$. Using MAVEN data, Inui et al. (2019) found a flux ratio of O${}_{\mathrm{2}}^{+}$ $/$ O${}^{+}=\mathrm{1.2}$ and reported that CO${}_{\mathrm{2}}^{+}$ contributed less than 10 % to the total heavy-ion flux. In summary, observations indicate that the flux of escaping O⁺ and O${}_{\mathrm{2}}^{+}$ ions are of similar magnitude and that the CO${}_{\mathrm{2}}^{+}$ flux is not significant.

As the code used here only handles one single-charged ionospheric ion species, we carry out separate simulation runs using O⁺ and O${}_{\mathrm{2}}^{+}$ to investigate the effects of composition on the ion escape rate estimates. Note also that there is no neutral H or O corona in the model. Thus, those populations of exospheric pick-up ions are missing. This means that the mass loading due to the corona is missing in the simplified model that we use. The total mass loading can be compensated for by more ions from the ionospheric source; however, in reality, the spatial distribution would be different. We also do not include any alpha particles in the solar wind, which should have some effect on the solar wind interaction.

2.3 Model parameters

The coordinate system used is MSO (Mars Solar Orbital) coordinates with the solar wind flowing along the −x axis, with density n_sw, speed v_sw, and temperature T_sw. The upstream interplanetary magnetic field is B_sw. The computational grid has cubic cells of size Δx, and the time step is Δt. The computational domain is $-\mathrm{11}\mathrm{000}\le x\le \mathrm{10}\mathrm{000}$, $-\mathrm{33}\mathrm{600}\le y\le \mathrm{33}\mathrm{600}$, and $-\mathrm{33}\mathrm{600}\le z\le \mathrm{33}\mathrm{600}$ km. On the upstream boundary, after each time step, we insert solar wind protons so that the number of particles per cell there is constant. In the y and z directions, we have periodic boundary conditions. The produced ionospheric ions have a weight (how many real ions they represent) that is chosen such that the weight is similar to that of the solar wind protons. The model parameters and their values are listed in Table 1.

Table 1 The parameters used for all simulation runs.

2.4 Observations

We use MAVEN magnetic field (Connerney et al., 2015; Dunn, 2021) and ion (Halekas et al., 2017) observations to determine upstream conditions and bow shock location. The orbit is chosen such that the solar wind conditions are steady (so that the conditions should be unchanged while MAVEN is inside the bow shock) when we do not have observations of the solar wind. This also allows us to use a simulation that does not have time-dependent upstream solar wind conditions. To verify our results, we also use MEX observations of electrons (Lundin et al., 2005; Frahm et al., 2006) to locate bow shock crossings. The upstream solar wind parameters are listed in Table 1.

2.5 Algorithm

We now describe the algorithm (illustrated in Fig. 1) for estimating the escape of ionospheric ions from observations of the upstream solar wind and the location of the bow shock:

1. We start with an observed state of the upstream solar wind – the magnetic field and the solar wind density, velocity, and temperature.

2. We then carry out several runs of the hybrid model for these upstream solar wind conditions, using different ionospheric ion production rates.

3. Next, we find the simulation run that has a bow shock location that best corresponds to the observed location. This could be done quantitatively (e.g., by a least squares fit), but here we visually compare the simulations and the observations.

4. The escape for the best fit simulation run is then computed, and this will be our estimate of the escape rate of ionospheric ions at the time of the bow shock observation.

Figure 1 An illustration of the algorithm to estimate the ion escape rate. The Sun is to the left, and Mars is the red and black disk. Using fixed upstream solar wind conditions from observations in the hybrid model, we vary the ionospheric heavy-ion (blue dots) production rate for different simulation runs. The bow shock location (green lines) in each simulation run is then compared to the observed bow shock location. The estimated escape rate corresponds to that of the simulation run that best fits the bow shock location.

3 Results

As an example of the proposed algorithm, we performed 10 simulations with the production rates $p_i=\mathrm{0.1},\mathrm{0.2},\mathrm{\dots },\mathrm{1.0}$ [cm⁻³ s⁻¹]. In Fig. 2, we present a comparison of MAVEN observations and two hybrid runs (R₁ and R₂) with an O⁺ ionosphere and different ion production rates that were judged to best fit the observations. We see that there is a fairly good agreement between the models and the observations at the bow shock and in the magnetosheath. However, closer to the planet, in the induced magnetosphere, the agreement is not as good. The magnetic field is much larger and more variable than in the model. This is not surprising because we have a model with a very simplified ionosphere, no magnetic anomalies, and low spatial resolution. Moreover, in the magnetosphere, the proton velocities in the model are much higher than observed. However, the proton density is very small, close to zero, in much of this region, as observed and as shown in the simulations. Near the exit from the induced magnetosphere, the observed density is larger than in the simulations, resulting in a similar dynamic pressure. Furthermore, the variability of B_z in the magnetosheath is smaller in the models than in the observations.

As expected and as seen when comparing the two model runs, the location of the bow shock moves outward when the ionospheric ion production is increased, resulting in a larger mass loading.

The escape rate never reaches a steady state due to the intrinsic variability of the induced magnetosphere. Therefore, we determine the escape rate by averaging the flux of ionospheric ions along −x in the simulation domain, from $x=-\mathrm{5000}$ km to the outflow boundary. We then average these computed escape values over time, between 200 and 590 s, with a time step of 10 s.

Figure 2 A comparison of MAVEN observations (blue) with two O⁺ model runs (R₁ in orange and R₂ in green) at 490 s. Panels (a), (c), and (e) show the magnetic field magnitude, proton velocity, and proton number density, respectively. Panels (b), (d), and (f) show the three respective magnetic field components (B_x, B_y, and B_z) in MSO coordinates. The location of the magnetosheath as seen from the observations is also indicated in gray. Thus, the induced magnetosphere is between the two gray regions.

Table 2 The different best fit simulation runs, the ionospheric ion species, the maximum ionospheric production rates, and the resulting escape rates.

In Table 2, we show the results of simulation runs R₁–R₄, performed for O⁺ and O${}_{\mathrm{2}}^{+}$ ionospheres with different ion production rates. For the best fit runs using O⁺, the approximate escape rate is 2.0×10²⁴, whereas it is 1.5×10²⁴ s⁻¹ for O${}_{\mathrm{2}}^{+}$. As we have a mixture of the two ion species in reality, the escape rate should be between these values, and we can estimate the actual escape as 2×10²⁴ s⁻¹.

We can note how the escape rate for the best fit simulation runs depends on the species of the escaping ionospheric ions. The best fit has 25 % less escaping O${}_{\mathrm{2}}^{+}$ ions compared with O⁺. It is not, however, directly proportional to the total mass of the escaping ions; in that case, we would expect a 50 % reduction. Thus, it is not only the amount of mass loading that characterizes the best fit, the dynamics of the escaping ions is also important.

Looking at the escape rate in Table 2 for the same species but for different ion production rates, we see that it is weakly dependent on the production rate. The two best fit runs in Fig. 2 have less than a 1 % difference in the escape rate.

To verify the location of the bow shock in the two best fit model runs, we also use MEX observations of the bow shock in electron data. In Fig. 3, we plot the proton number density from the same two hybrid runs as in Fig. 2, although now along the MEX orbit, along with the location of the bow shock crossings observed by MEX. The agreement is fairly good, even if the observed bow shock is a few minutes earlier than that seen in the model runs. This corresponds to a distance of a few hundred kilometers, which is comparable to the simulation grid cell size. One reason for this could be that we have not used an aberrated solar wind velocity (it flows along the −x axis). This should result in a tilt of the whole magnetosphere and bow shock.

Figure 3 The proton number density along the MEX orbit for the two best fit simulation runs (R₁ in orange and R₂ in green), at 490 s of simulation time. The blue vertical lines show the two bow shock crossings seen in MEX electron data.

4 Discussion

For our example case, we found an escape rate of 2×10²⁴ ions s⁻¹. This is in the range of recently published estimates for the escape rate at Mars (Ramstad et al., 2015; Brain et al., 2015; Dong et al., 2017).

An assumption is that, given upstream conditions, the mass loading determines the bow shock location. Although we find an escape rate for a single orbit that matches observed escape rates, this assumption needs to be tested in more detail. An ongoing investigation is to apply the method to a large number of orbits and verify that the model-estimated escape rates are consistent with observed escape rates.

The proposed method is directly applicable to unmagnetized planets. For magnetized planets, the bow shock location is mainly determined by the upstream solar wind and the strength of the dipole field. Escape at magnetized planets occurs in the cusp regions; how this affects the bow shock location as well as if the presented method could be adapted to magnetized planets would need further investigation.

The bow shock location has been found to depend on the location of the magnetic anomalies relative to the solar wind flow (Fang et al., 2017). It is unclear if this is because the fields “push out” the boundaries or because the fields increase ion escape. The latter may not require crustal fields in the model used in our algorithm, as the parameter that we vary is the amount of ions near Mars available to escape. If the crustal fields in a specific geometry enhance escape, this will be captured in the algorithm because the best fit bow shock will be further out; in contrast, if the crustal fields in a specific geometry depress escape, the bow shock will be closer to the planet.

We used a hybrid plasma model in the algorithm. It would, however, be possible to use another type of plasma model, such as a magnetohydrodynamic (MHD) model, that can predict the location of the bow shock for different amounts of mass loading by ionospheric ions, given upstream solar wind conditions.

5 Conclusions

In the past, ion escape has been estimated either by computer models or from observations. Models have the problem that every physical process has to be present in the model. In contrast, observations suffer from variability, requiring the averaging of data over years or even decades. The method proposed here uses a model and observations together. In this way, we overcome the difficulties of each approach. We then get an estimate of the escape using just one observation. A model is used to estimate a global property (ion escape) from a local observation (bow shock location). As we use a model, there are not the physical limitations – in terms of energy coverage and field of view – that are present in observations. In particular, low-energy escaping ions are difficult to observe.

This opens up the possibility of estimating escape during flybys of unmagnetized planets, in the past and in the future. It also allows for the estimation of escape twice per orbit given only a magnetometer and an ion detector. This enables detailed studies on how the escape depends on different parameters, which has been difficult in the past due to the years of observations needed to collect enough statistics. It also makes the study of escape during transient events, like extreme solar wind conditions, possible, which is something that is important for studies on escape in the past as well as escape at exoplanets.