Effects of a dipole‐like crustal field on solar wind interaction with Mars

A three‐dimensional four species multi‐fluid magnetohydrodynamic (MHD) model was constructed to simulate the solar wind global interaction with Mars. The model was augmented to consider production and loss of the significant ion species in the Martian ionosphere, i.e., H+, O2+, O+, CO2+, associated with chemical reactions among all species. An ideal dipole‐like local crustal field model was used to simplify the empirically measured Martian crustal field. Results of this simulation suggest that the magnetic pile‐up region (MPR) and the velocity profile in the meridian plane are asymmetric, which is due to the nature of the multi‐fluid model to decouple individual ion velocity resulting in occurrence of plume flow in the northern Martian magnetotail. In the presence of dipole magnetic field model, boundary layers, such as bow shock (BS) and magnetic pile‐up boundary (MPB), become protuberant. Moreover, the crustal field has an inhibiting effect on the flux of ions escaping from Mars, an effect that occurs primarily in the region between the terminator (SZA 90°) and the Sun–Mars line of the magnetotail (SZA 180°), partially around the terminator region. In contrast, near the tailward central line the crustal field has no significant impact on the escaping flux.


Introduction
The interaction between the solar wind and the Martian atmosphere has attracted wide attentions because, unlike Earth, Mars lacks an intrinsic global magnetic field. But crustal magnetic fields are present on Mars and play an important role in the planet's interaction with the solar wind. The strongest crustal fields of Mars are located primarily in the southern hemisphere; the strongest sources are at latitudes of 30°S and longitudes between 120°W and 210°W (e.g. Acuña et al., 1999).
As the supersonic solar wind approaches the vicinity of the Mars, the plasma flow is thermalized and deflected around the obstacle. A bow shock (BS) comes into being as the outermost boundary ; a magnetic pile-up boundary (MPB) forms due to compression and deflection of the interplanetary magnetic field (IMF) around Mars . Several influences affect the position and the shape of the Martian BS: the planet's crustal field (Edberg et al., 2008), variations associated with the Martian seasons, solar wind dynamic pressure, and the incident solar extreme ultraviolet (EUV) irradiance (Hall et al., 2016(Hall et al., , 2019. Results from global magnetohydrodynamic (MHD) simulation in conjunction with observations from MAVEN have demonstrated the existence of a Martian hemispheric asymmetry: the BS in the southern Martian hemisphere is displaced at a higher altitude due to the presence of a significantly stronger crustal field . The planet's MPB, which is defined as the altitude at which the solar wind thermal pressure balances with the total magnetic pressure (Holmberg et al., 2019), appears to be the transition region in which the planetary exosphere plays a major role in determining the plasma properties (Ma YJ et al., 2008). In the southern hemisphere, where crustal fields provide extra magnetic pressure (Dubinin et al., 2008), the MPB altitude can be as high as 400-500 km. Based on data from every Mars Global Sur-veyor (MGS) MPB crossing, Crider et al. (2002) reported dependence of the Martian MPB altitude on planetary latitude, revealing that the average value for dayside southern hemisphere crossings is 200 km higher than for dayside northern hemisphere crossings; they found also that the terminator distance of the MPB is around 5180 km near the crustal field and 5040 km outside of the region.
Many researchers have constructed and improved different threedimensional models to simulate solar wind interaction with Mars in terms of physical assumptions, boundary condition treatments, and implementation schemes. Among the various numerical models that have been developed are multi-species single-fluid MHD models, multi-species multi-fluid MHD models, single-fluid MHD-with-test-particle models, hybrid models, and particle-in-cell models. They have been used to study the rate and mechanism of ion escape on Mars under specified conditions or during a long history (e.g. Ma YJ et al., 2015;Terada et al., 2009;Dong CF et al., 2018), the planet's magnetic field topology and structure (e.g. Liemohn et al., 2007;Kallio et al., 2006;Fang XH et al., 2018), the influence of crustal field rotation (e.g. Ma YJ et al., 2014a;Fang XH et al., 2015), plasma boundaries (e.g. Simon et al., 2007;Harnett and Winglee, 2006;Xu SS et al., 2016), solar wind impact on the ionosphere (e.g. Harnett and Winglee, 2007;Ma YJ et al., 2014b;Dubinin et al., 2017) as well as phenomena in the magnetotail, such as twist structure and reconnection (e.g. DiBraccio et al., 2018;Ma YJ et al., 2018).
Ancient Mars was a warm and wet planet (Lasue et al., 2013). The reason it became cold and dry has been a topic of interest for many years. Since Mars has no global intrinsic magnetic field, its atmosphere is exposed to, and interacts directly with, the solar wind. Several mechanisms have been proposed to interpret the ion escape from Mars, including Jeans escape (Chassefière and Leblanc, 2004), photochemical escape Lillis et al., 2015Lillis et al., , 2018Fox and Hác, 2014), ion sputtering (Johnson and Leblanc, 2001), ionospheric outflow (Fränz et al., 2015) and ion pick-up (Rahmati et al., 2017(Rahmati et al., , 2018 The presence of crustal fields on Mars complicates understanding of the planet's interaction with the solar wind and has an important impact on explanations of ion escape from Mars (e.g. Ramstad et al., 2016;Dong CF et al., 2015;Ma YJ et al., 2014a;Ledvina, 2012, 2014). It has been suggested that the crustal field on Mars may have two opposite effects on ion loss. One is a shielding effect, affecting solar wind penetration and atmospheric/ionospheric stripping (Brecht and Ledvina, 2014;Nilsson et al., 2011). The other is an escape-fostering effect acting on the transterminator ion flux area and day-night connection (Fang XH et al., 2017).
In this paper, a three-dimensional multi-fluid magnetohydrodynamic (MHD) simulation was performed to clarify the role played by the Martian crustal field on the boundary layer location and shielding effect region during the solar wind-Mars interaction. Because it is not easy to distinguish many physical phenomena from the crustal magnetic field, owing to its complicated nature and distribution on Mars, our simulation employs a simple dipole magnetic field model with magnitude comparable to that of the crustal field. This simplification has the advantage, compared to models based on empirical measurements of the crustal magnetic field, that different processes can be isolated by refining specific Martian remnant magnetic field models, revealing more easily the physics underlying the role of each element.

Model Description
The multi-fluid MHD model used in this paper is composed of Navier-Stokes (NS) equations for the significant ion species in the Martian ionosphere, which include conservative equations for the plasma flow with respect to continuity, momentum, and energy. The NS equations for individual ion species are augmented by the interaction of the electromagnetic effects.
Here we consider four main ion species in the Martian ionosphere, i.e., H + , O 2 + , O + , CO 2 + . The governing equations for individual species can be expressed as follows (Najib et al., 2011): where . Here and q s are the individual mass density, velocity, pressure of the ions, number density, and charge, respectively. is the magnetic field, J is the identity matrix. p e , n e are the pressure and number density of electrons, respectively. The source terms , , and in Equation (1) represent the variations due to the chemical reactions among all the species for mass, momentum, and energy respectively. The inelastic collisions considered in this paper include charge exchange, photoionization, and recombination. The source terms can be written as follows: u u u n where S s and L s denote production and loss rate for species s, respectively; ν st is the collision frequency between species s and t; and T n represent velocity and temperature for neutrals, respectively; T s and m s are temperature and mass for ion species, respectively; κ is the Boltzmann constant, and γ is the ratio of specific heats (taken to be 5/3).
The neutrals in the atmosphere considered in this paper include CO 2 , O 2 , H; their densities are taken from the model calculations of Bougher et al., (2000). Spherical symmetric profiles of neutrals are used for simplification. Since the model takes account of chemical reactions, including photoionization, charge exchange, and dissociative recombination, a self-consistent three-dimensional ionosphere could be generated. The associated chemical reactions and the corresponding reaction rates that we consider in the model are adopted from Ma YJ et al. (2004) and Schunk and Nagy (2009).
In the multiple fluid MHD model proposed by Najib et al. (2011), the cosine function of the solar zenith angle was used to obtain the optical depth effect in calculating photoionization rates of the neutrals, which would make the solar flux be zero on the night side. In fact, the Chapman function was proved to describe more precisely the M2 layer of the Martian ionosphere (Withers, 2009), so we use it to calculate photoionization rates in our model.
Our calculations are performed in the Mars-centered Solar Orbital (MSO) coordinate system, with x axis pointing from Mars towards the Sun, z axis perpendicular to the x axis and positive toward the north celestial pole, and y axis completing the right-handed coordinate system. The computational domain is set to be , where R M is the radius of Mars (R M = 3396 km).
The second order entropy conditioned scheme (Dong HT et al., 2002) has been proven to be a highly accurate and robust procedure for solving 3-D Euler equations, because the scheme requires no adjustable parameters to guarantee satisfaction of the entropy condition. The finite difference discretization of Equation (1) is performed under a general curvilinear coordinate system (ξ, η, ζ). Accordingly, Equation (1) can be written symbolically as where covers the conservative variables, F, G, H are the convective fluxes, and F v , G v , H v the viscous fluxes.
As in conventional fluid dynamics, we temporally disregard the source terms in Equation (2), and consider only the case of ideal MHD flow. The finite difference discretization of Equation (2) for ideal MHD flows can be deduced by a splitting procedure as (3) The numerical flux , , and were proposed by the entropy conditioned procedure (Dong HT et al., 2002).
The procedure was implemented to solve all the non-homogeneous N-S equations for individual ion species, one by one, followed by calculation of the un-coupled magnetic induction equation using a second order MacCormack scheme.
Since the general curvilinear coordinate system was adopted and transferred from the physical MSO coordinate system, we obtain high resolution for the region with intense variations of physical (B x , B y , B z ) = (−1.6, 2.5, 0) nT parameters by refining the physical grid. The lowest boundary of the computational domain was set to be 100 km above the Martian surface; the smallest grid size was 30 km, and the O 2 + , O + , CO 2 + densities were taken to be values in photochemical equilibrium. The density of H + was set to be 0.3 of the solar wind density. The solar wind density and velocity were chosen to be 4 cm -3 and 500 km/s, respectively. The interplanetary magnetic field (IMF) follows a Parker spiral orientation of 56 degree in the x-y plane with a magnitude of 3 nT; that is, in the MSO coordinate system.
In order to simulate the interaction of solar wind with Mars in the presence of a local dipole magnetic field, we used a three-dimensional function of an ideal dipole-like magnetic field to mimic the interaction under the impact of the crustal magnetic field. The local dipole-like magnetic field can be expressed as follows: where M z and r, respectively, denote the magnetic moment along the z axis and the altitude above the Martian surface. The point in the MSO coordinate represents the position of the dipole field.
and (x 0 , y 0 , z 0 ) = (2500 km, 0, -1443 km) were chosen for this simulation, which means that the dipole field is located at 30°S and its magnitude at 200 km altitude is roughly 600 nT, which is consistent with observed values of the local crustal magnetic field.

Asymmetry Feature of Structures and Plasma Flow
In our model scenario the magnetic field configuration and plasma distribution during the solar wind interaction with Mars evolve self-consistently. A steady solution is obtained through converging in time with specific initial and boundary conditions. Figure 1 demonstrates an overview of the magnetic field lines in the x-y and x-z plane, superposed on the distribution of magnetic intensity in the no-dipole case. On the dayside, the typical largescale plasma configurations, such as detached bow shock and induced magnetosphere, are obviously apparent. The IMF pile-up and draping structure can be observed inside the magnetic pileup region (MPR) in the x-y plane of Figure 1b.
A south-north MPR asymmetry can be seen clearly in the meridian plane (Figure 1a), especially where hosts the strongest magnetic field. However, this asymmetric feature of the MPR in the meridian plane was not reported in results from the single-fluid multi-species MHD model (Ma YJ et al., 2004). Thus, it is reasonable to interpret the asymmetry as possibly caused by the movement of individual species, since the multi-species MHD model assumes that all ion species share the same velocity and temperature . Nevertheless, the MPR region in the equatorial plane is in a symmetric pattern. Figure 1 also indicates that the Martian MPB is the outer boundary of the place where the draping effect becomes prominent. Such a signature can be applied to identify the MPB (Bertucci et al., 2003a), particularly when the data of the magnetic field intensity are ambiguous (Nagy et , 2004). This approach was demonstrated in Bertucci et al. (2003b) and also applied to Venus.

(u u u + × B B B)
The convection electric field vectors are shown in Figure 2a, superposed on the magnitude of the average ion velocity in the meridian plane for the no-dipole case. It can be seen that the convection electric field points primarily northward, offering an explanation for the formation of a plume structure in the Martian magnetotail (Regoli et al., 2018;. In the southern hemisphere, ions are picked up and accelerated toward Mars by the convection electric field, resulting in precipitation in the atmosphere. However, in the northern hemisphere, the convection electric field sweeps ions away from Mar, leading to the formation of the plume flow. In the equatorial plane, as shown in Figure 2b, the velocity profile is symmetric about the x-axis, whereas Figure 2a shows it to be asymmetric about the x-axis in the meridian plane. Individual ion velocity decoupling could be responsible for the velocity asymmetry (Najib et al., 2011), which has not been detected in the single-fluid MHD model.
Plume flow of individual ion species is a potential channel for ion escape Dong CF et al., 2014;Fang XH et al., 2008). The contours of densities for all four individual ion species are plotted in Figure 3. The three main ion species with large fractions in the ionosphere, ionosphere, O 2 + , O + and CO 2 + possess highly asymmetric density distributions. The contours of O 2 + and CO 2 + are quite similar, whereas the contour of O + exhibits a slightly different shape. Najib et al., (2011) suggested that this difference may result from the ionization of the neutral oxygen corona, which pushes O + outwards. The Lorentz force exerted on the individual ion species is expressed as , where the electric field is calculated by using the generalized Ohm's law, which can be expressed as follows: It can be easily proven that the term can lead to an asymmetric flow in the meridian plane, as long as the magnetic  field is in the X-Y plane (Dong CF et al., 2014). On the other hand, from the particle simulation point of view the asymmetry can also be explained by the convection electric field (Fang XH et al., 2010). In fact, there is no essential difference between these two explanations except that the ions of a certain species move together in the multi-fluid MHD code, while they move individually in a particle simulation code. Therefore, the multi-fluid MHD model can capture the asymmetries of the ion flow arising from heavy ion acceleration by the convection electric field. By using MAVEN observation data,  suggested that the plume-like channel is a constant structure of O + escape, which suggests that our multi-fluid MHD simulation is a reliable method to study the configuration of the escaping flux of heavy ions.

Effects on the Locations of BS and MPB
The temperature of solar wind protons would be increased in the magneto-sheath and deceased upon the magnetic pile up region . In the magnetic pile up region, the magnetic pressure increases due to the compression of magnetic field. Here we assume, as a matter of simple gas-dynamics, that the bow shock is defined as the place where the dynamic pressure and thermal pressure are in balance, whereas the magnetic pileup boundary layer is the place where the thermal pressure is equal to the magnetic pressure. Figure 4a shows pressure distributions along the Sun-Mars line at subsolar points in the no-dipole case, which also includes magnetic pressure , thermal pressure of the solar wind P Tsw , and solar wind dynamic pressure P sw .
Closer inspection of Figure 4a shows that, in the absence of the di-pole field, and given the gas-dynamic consideration mentioned above, the transition of the dominated region from thermal pressure to magnetic pressure takes place at 1.28R M and the bow shock is located at 1.52R M .
In comparison, Figure 4b demonstrates the profiles, in the dipole case, of the pressure ingredients along the Sun-Mars line at subsolar points. Since the magnitude and distribution of the ideal dipole magnetic field model used in this paper is comparable to the practical crustal magnetic field descriptions suggested by several previous investigations (Purucker et al., 2000;Cain et al., 2003;Morschhauser et al., 2014), the influence of the crustal field in shaping the Martian plasma environment can be studied qualitatively. It can be seen in Figure 4b that our model locates BS and MPB at 1.56R M and 1.32R M , respectively, which is consistent with observations (Vignes et al., 2000. When the model includes the presence of a dipole magnetic field, the BS and MPB are pushed towards the sun. These behaviors are understandable because the dipole field increases the magnetic pressure in the ionosphere. Since the magnitude of the dipole field in our model was chosen to represent realistically the scale of the crustal field, we assert that our model with this dipole field simulation can be used to reveal crustal field influences on the interaction process between the solar wind and the Martian upper atmosphere.

Effect on Ion Escape Flux
The multi-fluid MHD model offers a practical, feasible, approach to calculating the number densities and velocities of the escape flux of individual ion species. Figure 5 presents the distribution of escape fluxes of O + , O 2 + and CO 2 + in the meridian plane for both nu- The lower panels in Figure 5 present the distributions of ion escape fluxes for O + , O 2 + and CO 2 + in the case of a dipole magnetic field. The most obvious discrepancy between the two cases is that the escape flux intensity in the region between the terminator and the tailwards Sun-Mars line decreases remarkably when the dipole field is present, resulting in two distinct escape channels. In order to check differences of escape flux in the plume region in detail, in Figure 6 we plot the variations of escape fluxes for heavy ions along a part of circles with radius of 7R M (white solid curve in the lower panel of Figure 5) in the plume region, where θ is defined as the rotation angle from 45° to 225° going counterclockwise from the solar zenith angle (SZA) of 45°. It is readily seen that in the no-dipole case the escape flux in the outflow region is much higher than that of the dipole case. However, a notable difference between the two cases occurs in the region between the terminator (θ = 90°) and the tailward Sun-Mars line (θ=180°), particularly around the terminator, whereas the difference in the tailward region is subtle. Our numerical results are in good agreement with observations .   Observations and simulations demonstrate that the crustal field has an inhibiting effect on the ion escape from Mars, and the terminator region and a channel around the tailward Sun-Mars line are two main escape channels Fang XH et al., 2010). Our simulation results indicate that the inhibiting effect of crustal field on the ion escape mainly occurs in the region between the terminator and the Sun-Mars line of the Martian magnetotail.

Summary
In this paper, a three-dimensional multi-fluid magnetohydrodynamic (MHD) model was established to study the solar wind global interaction with Mars under the influence of an ideal dipole-like crustal magnetic field. The model was augmented by consideration of production and loss of the significant ion species in the Martian ionosphere, i.e., H + , O 2 + , O + , CO 2 + , associated with chemical reactions among all species. The Navier-Stokes equations for individual ion species and the magnetic induction equation were calculated by entropy conditioned and MacCormack schemes, respectively, both in the second order. An ideal dipole magnetic field model with the same intensity and location as the strongest empirically measured crustal magnetic field was used.
Our simulations effectively reproduced the known large-scale phenomena, such as the bow shock (BS), the magnetic pile-up boundary (MPB), and the ionosphere et al. Since velocities and densities for all ion species were calculated separately, we are able to conclude that individual ion velocity decoupling causes the asymmetries of average ion velocity and magnetic field as well as ion density plume structure. Comparison between the cases of with/without dipole magnetic field indicate that the dipole magnetic field located in the southern hemisphere pushes the BS and MPB to higher altitudes towards the sun. Locations of the BS and MPB predicted in the dipole magnetic field case agree reasonably with observations and with results from previous simulations.
Our results confirm that the crustal field has an inhibiting effect on ion escape. In the absence of a dipole magnetic field, ions escape through a plume region of the one-fourth quadrant between the terminator and the Sun-Mars line. However, the dipole crustal magnetic field causes ion escape principally concentrated in the region around the terminator and the Sun-Mars line. That is to say, under the impact of a local dipole magnetic field, ion escape is severely suppressed in the transition region between the terminator and the Sun-Mars line, whereas the escape flux is only slightly decreased around the Sun-Mars line.