Effects of an Intrinsic Magnetic Field on Ion Escape From Ancient Mars Based on MAESTRO Multifluid MHD Simulations

Ion escape has played a key role in atmospheric loss on ancient Mars due to intense solar activity. Under the existence of a strong global intrinsic magnetic field, as expected on ancient Mars, differential flows between different ion species can be important for ion escape. To assess effects of differential flows, we here developed a new global multifluid magnetohydrodynamics model with a cubed sphere grid named Multifluid Atmospheric Escape Simulations Toward Real elucidatiOn (MAESTRO). A test simulation under present‐day Mars conditions showed solar wind‐Mars interactions, for example, plasma boundaries, ionospheric profiles, and tailward outflows, consistent with observations and simulation studies. We then conducted multifluid and multispecies simulations with six different dipole field strengths under ancient Mars conditions. The multifluid cases show asymmetric outflow distributions with respect to the solar wind convection electric field, as pointed out by previous studies on present‐day Mars. Compared with multispecies cases, the multifluid representation increases the escape rates of O2+ and CO2+ by more than two orders of magnitude in the strong dipole field cases where the cusp outflow is dominant. The O+ escape rate is slightly lower in the multifluid cases with no or weak dipole field due to suppression of ion pickup in the −E hemisphere, while it is reduced by one order of magnitude in the strongest dipole field case. The existence of a dipole field can reduce the total escape rate by a factor of six. The impact on atmospheric evolution is diminished but still significant.


Introduction
Mars experienced drastic climate change due to significant loss of the atmosphere during its early period, suggested by geological and mineralogical evidence (Jakosky & Phillips, 2001).Atmospheric escape to space is a key loss process of the Martian atmosphere since sequestration of the atmosphere into the crust is thought to have played a limited role due to the absence of massive and widespread carbonate deposits (Ehlmann & Edwards, 2014).The intense activity of the young Sun supports the importance of atmospheric escape.The upper atmosphere of Mars was exposed to more severe solar X-ray and extreme ultraviolet (XUV) radiation and solar wind (Ribas et al., 2005;Wood, 2006) that enhanced atmospheric escape via heating and ionization of the neutral atmosphere.In particular, escape of the ionized atmosphere, that is, ion escape, was the dominant escape process during the early period of Mars (Dong et al., 2018;Lammer et al., 2013).On the other hand, Mars has intense and localized crustal magnetic fields, which suggests a dynamo-driven planetary magnetic field on ancient Mars (Acuña & Connerney, 1999).Studies on the magnetization and age of large impact craters indicate that the dynamo ceased at around 4.1 Ga (Lillis et al., 2013;Vervelidou et al., 2017).However, recent MAVEN observations identified magnetic signatures above a pyroclastic flow, suggesting the active dynamo at 3.7 Ga (Mittelholz et al., 2020).The existence of an intrinsic magnetic field protects the atmosphere from the solar wind by forming a large magnetosphere and inhibiting direct interactions between the solar wind and the upper atmosphere (Lundin et al., 2007).However, this idea is questioned by the fact that Earth, Venus, and Mars experience comparable ion escape despite different magnetization states (Ramstad & Barabash, 2021, and references therein).
Recent simulation studies focused on the effects of the existence and strength of an intrinsic magnetic field.Multispecies magnetohydrodynamics (MHD) simulation studies revealed that the existence of a dipole magnetic field has opposite effects on ion escape depending on the dipole field strength (Sakata et al., 2020(Sakata et al., , 2022)).The effects are remarkable on molecular ion escape caused by cold ion outflow from the ionosphere.Molecular ion escape is increased when the dipole field is weak, while it is decreased by two orders of magnitude under the existence of a strong dipole field.The threshold strength of the two opposite effects is approximately 0.1 in the ratio of the magnetic pressure of the dipole field at the equatorial surface to the solar wind dynamic pressure.On the other hand, the O + escape rate changes little in the weak dipole field cases but decreases by an order of magnitude in the strong dipole field cases.The milder effects on O + escape are due to a large contribution of ion pickup from the extended oxygen corona.The O + escape rate accounts for more than 97% of the total (O + , O 2 + , and CO 2

+
) escape rate from ancient Mars.Egan et al. (2019) also studied the effects of a dipole field on ion escape from present-day Mars based on hybrid simulations.The escape rates of O + and O 2 + increase with the dipole field strength until the magnetopause standoff distance reaches the induced magnetosphere boundary.The increased standoff distance by a dipole field reduces pickup ion precipitation into the − E hemisphere, where the solar wind convection electric field points toward the planet, enhancing ion escape.On the other hand, the stronger dipole field decreases the ion escape rates by shielding the solar wind flow.
However, there are caveats in these simulation studies due to model limitations.The multispecies MHD model is based on the single fluid approximation, that is, it assumes that all ions have the same velocity and temperature.This approximation is inappropriate in regions where ion populations with different velocities and temperatures coexist.The cusp region on Earth is a good example.The magnetosheath plasma has direct access to the ionosphere through the cusp region (Newell & Meng, 1988).Meanwhile, the cusp region is one of the major outflow channels of ionospheric ions (Nilsson et al., 2012).A more realistic representation of outflows through such regions is required to evaluate ion escape at strongly magnetized planets.Another limitation is that the fluid approximation also ignores kinetic dynamics.Multifluid MHD models partly overcome these limitations by treating ion species as separate fluids.Multifluid MHD simulations show plume-like distributions of planetary ions along the solar wind convection electric field that cannot be depicted in multispecies MHD simulations (Najib et al., 2011).Comparison studies of multispecies and multifluid MHD simulations on present-day Mars showed that planetary ion escape, particularly molecular ion (O 2 + and CO 2

+
) escape, is a few times stronger in multifluid MHD simulations (Ma et al., 2019;Regoli et al., 2018).Enhanced escape of molecular ions in multifluid MHD simulations indicates that the single fluid approximation underestimates ion outflow from the ionosphere.Hybrid simulations include ion kinetics and succeeded in reproducing solar wind-Mars interactions, including polar plumes (e.g., Jarvinen et al., 2018).However, high computational costs due to the treatment of ion particles often restrict grid resolution and consideration of ionospheric processes.A hybrid simulation study by Egan et al. (2019) focused on magnetospheric physics and used a simplified representation of the ionosphere assuming ion injection with pre-defined production rates consistent with observations at present-day Mars.The quantitative evaluation of ion escape from ancient Mars requires a self-consistent calculation of photochemistry in the ionosphere.
In this study, we introduce a new global multifluid MHD model named Multifluid Atmospheric Escape Simulations Toward Real elucidatiOn (MAESTRO) for a more precise evaluation of the effects of an intrinsic magnetic field on ion escape from ancient Mars.The model simulates the solar wind-magnetosphere-ionosphere interactions self-consistently.We describe the implementation of the model in Section 2. Section 3 shows an example simulation of present-day Mars.The results of multispecies and multifluid simulations with different dipole field strengths are shown in Section 4, and the effects of an intrinsic magnetic field are discussed in Section 5. We concluded the study in Section 6.

Multifluid MHD Equations
Here, we describe the equation system in MAESTRO.Note that the equations are written in the CGS-Gaussian units, and letters with subscript s denote variables of ion species s.
The multifluid model simulates the fluid equations of each ion species that consists of the continuity equations for the mass density ρ s , the conservation equation of the momentum density M s = ρ s u s , and the conservation equation of the energy density ε s = ρ s u 2 s /2 + P s / (γ − 1) .
where n s , u s , P s , and q s are the number density, velocity, pressure, and charge, respectively; γ is the ratio of specific heats; c is the speed of light; g is the gravitational acceleration.
S ρ s , S M s , and S ε s in the right-hand side of the equations are the source terms related to photochemical and collisional processes.
where m s is the mass, and the temperature T s is defined as P s / (n s k) with the Boltzmann constant k.
S s and L s are the production and loss rates due to photoionization, electron impact ionization, ion-neutral reactions, and dissociative recombination, respectively.Photoionization rates are calculated as described in Section 2.2.Electron impact ionization is included for H + and O + .The impact ionization rate is n e n n ν EII , where n e , n n , and ν EII are the electron number density, neutral number density, and ionization frequency, respectively.The ionization frequency depends on the electron temperature and is taken from Cravens et al. (1987).ν st is the effective momentum transfer collision frequency for ion species s due to collision with species t.The model considers ion-ion collisions in addition to ion-electron and ion-neutral collisions.u n is the velocity of the neutral atmosphere and is set to be zero in this study.The model considers five ion species: planetary H and solar wind H + .The details of reactions and collisions are described in Text S1 in Supporting Information S1.
To update the energy density of ion fluids, the model solves Equation 3 except for the photochemical and collisional terms.In calculating these terms, the model updates the pressure using the following source term for the pressure S P s instead of updating the energy density to avoid negative pressure.
The energy density is then recovered from the updated mass density, momentum density, and pressure.
For the magnetic field B, the model also solves the induction equation where the electric field E is defined by the generalized Ohm's law, including the Hall, resistive, and electron pressure gradient terms.
The resistivity η is calculated as m e ν et /n e e 2 where ν et is the total effective momentum transfer collision frequency for electrons.u + is the charge-averaged velocity of ions defined with the electric charge e: The electron number density n e and the electron velocity u e are calculated from the charge neutrality and the definition of the current. ) where the current density J is In addition, the electron pressure P e is calculated by solving the following pressure equation: K is the electron thermal conductivity coefficient used in Terada, Shinagawa, et al. (2009).
S P e includes electron heating and cooling terms in addition to photochemical and collisional terms.S e,ph is the production rate of photoelectrons by photoionization.The heating rate per photoelectron is expressed by E ph and taken to be 1.0 eV, as with Ma et al. (2019).This value is estimated based on observations of suprathermal electrons by MAVEN (Sakai et al., 2016).L e,rc kT e is the energy loss by dissociative recombination where T e is the electron temperature.L EII kT EII represents the energy loss by electron impact ionization, and the energy loss for each impact ionization is assumed to be 15 eV.L e,R(CO 2 ) and L e,V(CO 2 ) are cooling rates by inelastic collision (rotational and vibrational) with CO 2 , respectively.The formulae of the inelastic collision cooling terms are shown in Text S1 in Supporting Information S1.
The magnetic field is split into two parts B = B 0 + B 1 to improve the accuracy of simulations with a strong background field, such as an intrinsic magnetic field.B 0 is the divergence-free, curl-free, and steady background field.The model solves the fluctuating field B 1 .In this study, an intrinsic magnetic field is given as B 0 .

Numerical Implementation
We adopted the cubed sphere grid system (Ronchi et al., 1996).The cubed sphere grid consists of six identical faces, each of which is a projection of a cube's face onto the surface of the circumscribed sphere.It has a quasiuniform grid resolution over the whole sphere and can overcome numerical difficulties in the latitude-longitude grid due to the convergence of the meridians.It also provides computational advantages in parallelization because of six faces with an identical coordinate system.Figure 1 shows the six faces and the coordinate system of the cubed sphere grid.The coordinate lines in each face are arcs of great circles on which the grid points are spaced at equal angles.The grid resolution can be specified by the number of grid points along the coordinate line, represented as N. We used the grid resolution of N = 25 in this study.The angular resolution is approximately 3.6°near the center of the face and 2.5°near the vortices of the face.The formulae for the cubed sphere grid are shown in Text S2 in Supporting Information S1.
We placed 400 grid points non-uniformly in the vertical (radial) direction.
The resolution is 3 km near the inner boundary at 100 km altitude and increases with r, the distance from the center of the planet.The resolution is around 30 km at r = 1.5 R p , the typical subsolar location of the bow shock at present-day Mars, and approximately 3,000 km near the outer boundary at r = 40 R p , where R p is the radius of Mars.The total number of grid points in the simulation domain is 1.5 million.The simulation is conducted in the Mars-centered Sun Orbital (MSO) coordinate, in which the X axis points from the center of Mars toward the Sun, the Y axis points antiparallel to the orbital velocity of Mars, and the Z axis completes the right-hand system.
In MAESTRO, the equations are solved with a finite difference scheme.The numerical fluxes are calculated based on the semi-discrete second-order central scheme with a minmod limiter (Kurganov & Tadmor, 2000).The scheme requires no characteristic information other than the fastest wave speed to calculate the numerical fluxes.It is, therefore, convenient for the multifluid MHD system in which it is difficult to derive the eigenvalues and eigenvectors.The fastest wave speed is chosen among the fast magnetosonic speed for the total ion fluid and the sound speed of each ion fluid in this study.The time integration is conducted by the third-order total variation diminishing Runge-Kutta scheme (Shu & Osher, 1988), which guarantees that the extrema and total variation of the solution do not increase.The time step is calculated by the Courant-Friedrichs-Lewy (CFL) condition with a coefficient of 0.2.Note that the inclusion of the Hall term results in the whistler waves with the fastest wave speed inversely proportional to the grid resolution.This term is important in a spatial scale smaller than the ion inertia length.The typical ion inertial length in the simulations of this study was a few kilometers in the ionosphere and a few tens kilometers in the sheath region, both of which were comparable to or smaller than the grid resolution.The Lorentz force term in the momentum equations may also be numerically unstable.Since its characteristic time scale is the ion gyro period, the model limits the time step by the shortest proton gyro period.However, the time step determined by the CFL condition was typically more than one order of magnitude smaller than the shortest proton gyro period in this study.
In three-dimensional MHD simulations, some treatment of numerical divergence of the magnetic field is required.We adopted the composite diffusive/Powell method introduced by Liu et al. (2021).It is the combination of the diffusive method and the Powell method.In the diffusive method (e.g., Crockett et al., 2005), the magnetic field is modified at each time step by the following equation: which leads to diffusion of the divergence of the magnetic field.μ should be determined so as not to affect the CFL conditions and is set to be 0.4 in this study.The Powell method (Powell et al., 1999) adds the source terms proportional to the divergence of the magnetic field for the momentum and energy equations of ion fluids and the induction equation.Due to these additional terms, the divergence of the magnetic field is transported by the plasma motion and goes away from the simulation domain.The composite diffusive/Powell method can be implemented with local spatial derivatives, and it requires no additional variables.It can reduce the divergence of the magnetic field as well as other methods, such as the 9-wave method by Dedner et al. (2002).

Multispecies MHD Equations
MAESTRO is also able to perform multispecies MHD simulations.The equations for the mass density ρ, the momentum density M = ρu, the energy density U = ρu 2 /2 + P/(γ − 1) + B 2 /(8π), and the mass density of ion species sρ s are as follows.
The induction equation for the magnetic field includes the Hall and resistive term.
where the electric field is defined by The electron pressure gradient term, or the ambipolar electric field, is ignored because this term is much smaller than other terms in this study.Note that the electron pressure gradient force acts on the bulk fluid as part of the thermal pressure gradient force in Equation 18, while it acts on ion fluids as the ambipolar electric field in multifluid simulations (see Equations 2 and 9).
The electron pressure P e is also calculated by solving the electron pressure equation: The ion pressure is calculated by subtracting P e from the total pressure P.
The source terms due to photochemical and collisional processes in the right hand side of the equations (S ρ s , S M s , S ε s , S P e ) are the same as those in the multifluid MHD equations, but the ion velocities are the bulk velocity, that is, u s = u.In calculating the source term, the model updates the mass densities of ion species ρ s based on Equation 20 and then adds them up to get the mass density ρ.For the other terms, the model solves the continuity equation for ρ (Equation 17) and adjusts each ion mass density to match the total mass density.

Simulation Setup
The neutral atmosphere and photoionization rates are given as model inputs.The neutral atmosphere profiles used in this study are shown in Figure 2.For present-day Mars simulations in Section 3, we used the neutral atmosphere and photoionization rates in Shinagawa and Cravens (1989).The neutral atmosphere is based on the observations by the Viking spacecraft under low solar activity (Fox & Dalgarno, 1979;Nier & McElroy, 1977), and the photoionization rates are calculated using the solar minimum XUV flux (Cravens et al., 1980;Nagy et al., 1980).These profiles were also used in recent simulation studies (Sakai et al., 2018(Sakai et al., , 2021(Sakai et al., , 2023;;Terada, Kulikov, et al., 2009).For ancient Mars simulations in Section 4, we used a CO 2 -dominated atmosphere modeled by Kulikov et al. (2007), as with Sakata et al. (2022).Their model was a one-dimensional hydrostatic model considering the thermal balance and vibrational kinetics of molecules.The photoionization rates are calculated using the ancient solar XUV flux, which was assumed to be 50 times higher than the current flux.The current flux was based on the EUVAC model, assuming the F10.7 index of 150 (Richards et al., 1994).Both in present-day and ancient Mars simulations, the neutral atmosphere was assumed to be spherically symmetric and timeindependent.The photoionization rates were also time-independent but had a solar zenith angle (SZA) dependence in the form proportional to the cosine of SZA to take the optical depth effect into account.
In the present-day Mars simulation case, Case PM, we used typical solar wind parameters in Ma et al. (2004).The solar wind conditions were the proton number density of 4 cm − 3 , the velocity of 400 km s − 1 , the ion and electron temperatures of 1.75 × 10 5 K, and the magnetic field of 3 nT with the Parker spiral angle of 56°.Mars was assumed to have neither intrinsic nor crustal magnetic fields.
In the ancient Mars simulation cases, the solar wind proton number density, velocity, and ion and electron temperatures were 700 cm − 3 , 1,400 km s − 1 , and 6.5 × 10 5 K, respectively.The magnetic field was 20 nT with the Parker spiral angle of 60.9°.These values are taken from Sakata et al. (2022).We conducted six multifluid MHD simulation cases MF1-6.The intrinsic magnetic field was a dipole field with no tilt angle; the dipole field is purely northward at the equator.In Cases MF1-6, Mars had different dipole field strengths of 0, 100, 500, 1,000, 2,000, and 5,000 nT at the equatorial surface, respectively.For comparison, we also conducted six multispecies MHD simulation cases MS1-6 using the same neutral atmosphere model, the solar XUV, the solar wind, and the dipole field conditions.The crustal magnetic field was not included.
The input parameters are summarized in Table 1.

Test Simulation on Present-Day Mars
We performed a test multifluid MHD simulation to see the model's ability to reproduce solar wind-Mars interactions.We assumed typical conditions on present-day Mars, Case PM, described in the previous section.
Figure 3a shows the subsolar profiles of the total thermal pressure, the total dynamic pressure, the magnetic pressure, and the sum of these pressures.The bow shock is seen as a sharp transition from the dynamic pressure to the thermal pressure at around 1.45 R p .The interplanetary magnetic field (IMF) piles up in front of the planet, forming the magnetic pileup boundary at around 1.14 R p .The pileup field penetrates into the lower ionosphere because the peak thermal pressure in the ionosphere is lower than the solar wind dynamic pressure.These features agree well with previous simulation studies (Ma et al., 2004;Najib et al., 2011).The plasma boundaries are closer to the planet than the average location based on observations (Trotignon et al., 2006).As the boundaries are located further out over strong crustal field regions (Edberg et al., 2008), the closer boundaries in this test simulation can be attributed to the absence of crustal magnetic fields.
Figure 3b compares the ion number density profiles at an SZA of 43.2°on the northern and southern meridian planes.Note that the input neutral atmosphere profiles and production rates are based on the Viking observations at a SZA near 44°.Since the IMF is in the away sector of the Parker spiral, the convection electric field of the solar wind points toward the +z direction.Therefore, the Mars Solar Electric (MSE) coordinate coincides with the MSO coordinate, that is, the northern (southern) hemisphere coincides with the +E (-E) hemisphere, where the convection electric field due to deflected flows around the planet points away from (toward) the planet.The test simulation reproduced the ionospheric profiles consistent with the Viking observations and previous onedimensional simulations (Hanson et al., 1977;Shinagawa & Cravens, 1989) photochemical processes are crucial.On the other hand, the upper ionosphere above ∼350 km altitude shows asymmetric profiles with deeper penetration of the solar wind proton and a more depleted ionosphere in the +E hemisphere.The depletion of the ionosphere in the +E hemisphere is due to the polar plume escape of ionospheric ions (e.g., Sakakura et al., 2022) caused by the penetrating solar wind electric field.Such asymmetry in the upper ionosphere was also found by observations (Dubinin et al., 2018).The pileup region is between 200 and 350 km altitude, and the ionospheric flow is almost horizontal and toward the terminator.The ionospheric flow is slower in the +E hemisphere due to a stronger pileup of the IMF, which causes different trends from those seen at higher altitudes.
We also compared the cold ion outflow (tailward escape in the induced magnetotail) with a statistical study by Inui et al. (2019).They used the MAVEN measurements from July 2015 to December 2017 and analyzed O + and O 2 + outflows in the wake region.We compared the simulation results with their statistical analysis in the MSE coordinate to eliminate the effects of crustal magnetic fields.Figure 4 shows the O 2 + number density (n), the tailward velocity (− V x ), and the tailward flux (− nV x ) on the x-z plane.The color scale of each plot is similar to those used in Figure 5 of Inui et al. (2019).The ion distributions in the wake region are asymmetric between the +E and − E hemispheres.A dense O 2 + population extends to the tail in the − E hemisphere.The flow is correspondingly slower than the ambient plasma, but the tailward flux is higher in the − E hemisphere.These asymmetric features in the wake region are consistent with Inui et al. (2019).A denser and slower ion outflow with higher flux in the − E hemisphere is also pointed out by MAVEN observations and hybrid simulations in   Note.The solar wind parameters are the proton number density, velocity, interplanetary magnetic field (IMF) strength, and ion and electron temperatures.The IMF is in the away sector of the Parker spiral with the Parker spiral angle θ.B eq , the dipole field strength at the equatorial surface.XUV, X-ray and extreme ultraviolet; PM, present-day Mars case; MF, multifluid MHD case; MS, multispecies MHD case.Dubinin et al. (2019).They indicated that an asymmetric pileup of the IMF pushes the ionospheric plasma into the − E hemisphere and forms the ion trail, a dense and slow tailward flow in the − E hemisphere.The simulation results of Case PM show that the J × B force by the draped magnetic field accelerates ion outflows from the ionosphere toward the − E direction in the nightside of the +E hemisphere, which is consistent with the previous study.
The outflow fluxes of O + and O 2 + are also consistent with the observations.The outflow fluxes are calculated by averaging outward fluxes in the wake region, that is, on the part of the sphere which satisfies r = 2 R p , y 2 + z 2 < 1 R p , and x < 0. The average wake outflow fluxes of O + and O 2 + in Case PM are 1.06 × 10 6 cm − 2 s − 1 and 2.65 × 10 6 cm − 2 s − 1 , respectively.These values agree with the median fluxes in Table 1 of Inui et al. (2019): 1.04 × 10 6 cm − 2 s − 1 for O + and 2.26 × 10 6 cm − 2 s − 1 for O 2 + .The test simulation under typical conditions on present-day Mars indicates that the model can reproduce solar wind-Mars interactions, at least qualitatively, including the plasma boundaries, ionospheric profiles, and cold ion outflows, consistent with observations and simulation studies.In the next section, we applied this model to the simulations of ancient Mars under the intense solar XUV flux and solar wind to study the effects of a dipole field on ion escape.

Global Description of the Plasma
Figure 5 shows the total thermal pressure, that is, the sum of all ion pressures and electron pressure on the x-z (meridian) plane.The multifluid Cases MF1, MF3, MF5, and MF6 are in the right column, and the multispecies Cases MS1, MS3, MS5, and MS6 are in the left column.The super-fast magnetosonic solar wind flow is decelerated and compressed at the bow shock, and the dynamic pressure of the solar wind is converted to the thermal pressure in the magnetosheath.A large part of the thermal pressure is accounted for by the ion pressure in the multispecies results and by the solar wind H + pressure in the multifluid results.The thermal pressure in the magnetosheath is balanced by the magnetic pressure of the pileup IMF at the magnetic pileup boundary or that of the intrinsic magnetic field at the magnetopause.The distribution of the total thermal pressure, therefore, indicates the global configuration of the magnetosphere.
The multifluid cases show almost the same pressure distribution as the corresponding multispecies cases, indicating that the multifluid representation does not affect the global configuration of the magnetosphere.Case MF1 with no dipole field shows typical characteristics of the induced magnetic field; the magnetic pressure of the pileup IMF balances the thermal pressure in the magnetosheath.The ionosphere is strongly magnetized because the maximum total thermal pressure is ∼10 nPa, much lower than the solar wind dynamic pressure (about 2,300 nPa).In other cases with a dipole field, the existence of a dipole field changes the magnetospheric configuration.It pushed out the magnetopause, particularly at high latitudes, where the dipole field is stronger, and the solar wind dynamic pressure is reduced due to a large incident angle.The dipole field at high latitudes forms the lobe region of the magnetotail that consists of open magnetic field lines.Note that the topology of magnetic field lines will be discussed in Sections 4.2.2 and 4.2.3.In the strong dipole field cases (Figures 5e-5h), the compressed dipole field near the subsolar inhibits the penetration of the hot magnetosheath plasma.The cusps characterized by the precipitation of the hot magnetosheath plasma into the ionosphere are formed at intermediate latitudes between the polar and subsolar regions in Cases MF6 and MS6, the strongest dipole field cases.These changes in the thermal pressure distribution indicate the transformation from the induced magnetosphere to the Earth-like magnetosphere.
In contrast to the pressure distributions, planetary ion populations are quite different between the multifluid and multispecies cases.distributions in the multifluid cases are asymmetric, with a plume-like structure extending in the +E hemisphere.The population in the tail region is denser in the multifluid Cases MF1 and MF3 than in the corresponding multispecies cases.In Cases MF5 and MF6, the strong dipole field cases, there is a dense O 2 + population throughout the magnetosphere with a tenuous plume-like structure.The pattern of outflow fluxes will be investigated in the following subsections.

Cases MF1 and MS1
Figure 7 summarizes the tailward fluxes of O + , O 2 + , and CO 2 + and the topology of magnetic field lines at x = − 2 R p for Cases MF1 and MS1, no dipole field cases.The tailward flux is defined as − n s V x,s , where n s and V x,s are the number density and the x component of the velocity of ion species s, respectively.The plasma velocity V is used instead in the multispecies results.The white dashed line is a contour line of B x = 0, indicating the neutral sheet.
The topology of magnetic field lines is defined based on the connection to the planet: closed lines with both ends connected to the planet, open lines with one end connected to the planet and the other connected to the IMF, and draped lines with both ends connected to the IMF.The topology plots are also colored by the sign of B x , that is, the direction of the magnetic field: positive (negative) B x indicates the Sunward (tailward) direction.
The magnetic topology and the neutral sheet are almost identical between the two cases.It indicates that the multifluid representation does not affect the magnetospheric configuration, which is consistent with nearly the same distribution of the total thermal pressure shown in Figure 5.The magnetic field lines are draped lines in the whole region, and the neutral sheet is along the z-axis, that is, the convection electric field of the solar wind.This configuration is consistent with typical features of an induced magnetosphere under the absence of the dipole field and draping of the Parker-spiral IMF.The neutral line in the dawnside sheath region is formed due to the negative (non-zero) x component of the IMF.
In the multispecies Case MS1, the tailward fluxes of molecular ions (O 2 + and CO 2

+
) have two peaks around the neutral sheet.These fluxes are ionospheric outflows from the polar region driven by the J × B force of the draped IMF.The outflows from the polar region are also seen in Figure 6a.A ring-shaped strong flux exists in the flux distribution of O + in addition to the flux peaks around the neutral sheet.It arises from the ion pickup process in the extended oxygen corona under the strong solar XUV flux condition.An O + ion created on the sheath or solar wind region will be carried along by the solar wind flow, forming flux peaks around the planet.These results are consistent with the cases with no dipole field in the previous study with multispecies MHD simulations (Sakata et al., 2022).
On the other hand, Case MF1 shows strong asymmetric distributions with a strong flux elongating in the +E hemisphere.This flux corresponds to the plume-like structure in the number density plot on the meridian plane (see Figure 6b).The multifluid representation allows planetary ion fluid to be accelerated by the solar wind convection electric field.The plume-like flow is more tenuous but faster than the outflow in the tail region, resulting in a comparable tailward flux.Instead, the ring-shaped flux structure is indistinct due to the asymmetric flow.The plume-like flow extends to the − y direction in the sheath region (z > +3 R p ).This is due to the Lorentz force toward the − y direction acting on the outflow toward the +z direction in the IMF with negative B x .The multifluid representation also enhances ionospheric outflows around the neutral sheet.The flux distributions in the tail region have stronger peaks and a much weaker return flow, shown as the black region, compared to the multispecies case.The enhancement of ionospheric outflows is distinct in molecular ions.

Cases MF5 and MS5
Cases MF5 and MS5 were conducted with the dipole field strength of 2,000 nT at the equatorial surface.Figure 8 shows the tailward fluxes and the magnetic topology for Cases MF5 and MS5.The format of the figure is the same as Figure 7.
As in Case MS1, Case MS5 shows similar results to the strongest dipole field cases with a strength of 2,000 nT at the equatorial surface in the previous study (Sakata et al., 2022).The dipole field is dominant in the wake region, as indicated by a neutral sheet parallel to the equatorial plane.The magnetic topology also shows the formation of an Earth-like magnetosphere.The magnetosphere is characterized by a large region of closed magnetic field lines and lobe regions with open magnetic field lines.The neutral sheet extended into the northern hemisphere indicates  + , (e, f) CO 2 + , and (g, h) the magnetic field line topologies at x = − 2 R p for Cases MF5 and MS5.The format is the same as Figure 7.
the magnetopause, the boundary between the Sunward intrinsic magnetic field and the tailward IMF.It also extends in the southern-dawn side due to the kink of the IMF with positive B x just outside the magnetopause.In weaker dipole field cases, such as Cases MS3 and MS4, the neutral sheet is twisted like a "hybrid magnetosphere" recently studied in MAVEN observations (Dubinin et al., 2023;Xu et al., 2023).The twisted neutral sheet under the existence of a weak dipole field was also studied in our previous simulation studies (Sakata et al., 2020(Sakata et al., , 2022)).
The tailward fluxes of molecular ions in Case MS5 are suppressed strongly and have weak peaks in the flank region.These peaks in the flank region are due to the magnetic reconnection between the IMF and the dipole field.The open field line formed by the reconnection is carried along by the solar wind flow near the dayside ionosphere, and planetary ions escape along the field line.The detailed formation mechanism of these peaks was investigated by Sakata et al. (2020).The ring-shaped peak of O + flux is still visible but weakened compared to Case MS1 because the magnetosphere is more expanded.
While the magnetic configuration does not change, the tailward fluxes of molecular ions are much stronger in Case MF5.It suggests the enhancement of cold ion outflow through the lobe region.In addition, the fluxes have asymmetric distribution with a plume-like structure in the +E hemisphere like the no dipole field Case MF1.It inclines the +y direction because the flows through the cusp and the mantle region are accelerated by the duskward convection electric field.

Cases MF6 and MS6
Figure 9 summarizes the tailward fluxes and the magnetic topology for Cases MF6 and MS6.The two cases were conducted with the dipole field strength of 5,000 nT at the equatorial surface.The pressure ratio of the dipole field at the equatorial surface to the solar wind dynamic pressure is 4.33.This value is higher than the strongest dipole field cases in the previous studies: 1.49 in Case 5 of Sakata et al. (2020) and 0.694 in Case A5 of Sakata et al. (2022).
As in other cases, the magnetic topology and the neutral sheet are similar between the multifluid and multispecies cases.The magnetosphere is characterized by the closed field lines surrounded by the lobe region with open field lines, but its size is larger than Cases MF5 and MS5.Compared to the unmagnetized Cases MS1 and MF1, the size of the magnetotail is nearly doubled in both y and z directions.
There are only very weak tailward fluxes of O + , O 2 + , and CO 2 + in Case MS6.The flux peaks in the flank region do not exist, unlike Case MS5.This is because the magnetic reconnection occurs further away from the planet, and the escape channel responsible for the peaks is ineffective.In contrast, the multifluid Case MF6 shows enhanced fluxes with the asymmetric structure in the +E hemisphere.The size of the flux peaks also increases with the expansion of the magnetosphere.The ring-shaped structure in the tailward O + flux is absent, and the O + flux distribution is more similar to that of molecular ions than Case MF5.The expansion of the magnetosphere prevents the solar wind flow and reduces escape from the extended corona.Instead, cold ion outflow from the ionosphere remains the major escape process of O + , as well as molecular ions.

Ion Escape Rates
The escape rates of O + , O 2 + , and CO 2 + are calculated by integrating the flux on the sphere of r = 15 R p .The escape rates are summarized in Table 2 and plotted in Figure 10.The gray region shows the range where the magnetic pressure of the dipole field at the planetary surface is comparable to the solar wind dynamic pressure.The solar wind dynamic pressure was 2,300 nPa, which corresponds to the magnetic pressure with a strength of 2,400 nT.The magnetic pressure of the dipole field is comparable to the dynamic pressure at the equator when the dipole field strength at the equatorial surface (B eq ) is 2,400 nT and at the pole when B eq is 1,200 nT.
The escape rates of O + and molecular ions (O 2 + and CO 2 + ) in the multispecies cases have different dependence on the dipole field strength.Compared to Case MS1 with no dipole field, the escape rates of molecular ions are increased in Cases MS2 and MS3, the weak dipole field cases where B eq is less than 500 nT.The increase is up to 104% for O 2 + and 42% for CO 2 + .On the other hand, the escape rates are strongly reduced when the dipole field is sufficiently strong, that is, B eq is more than 1,000 nT.The escape rates in Case MS6 are approximately three orders of magnitude lower than those in Case MS1.The escape rate of O + shows similar but more moderate dependence on the dipole field strength.Case MS2 shows nearly the same O + escape rate as Case MS1, while the O + escape rate is decreased by approximately two orders of magnitude in Case MS6.Such dependence on the dipole field strength was also seen in our previous study with REPPU-Planets (Sakata et al., 2022).
The multifluid representation enhances molecular ion escape, and the enhancement is more remarkable in the strong dipole field cases.The increase in the molecular ion escape rates is a factor of about three between Case MF1 and MS1, but more than two orders of magnitude between Case MF6 and MS6 O + shows different effects of the multifluid representation.The O + escape rate is enhanced by about one order of magnitude in the strongest dipole field case.However, the O + escape rate is slightly lower in Case MF1-4 than in the corresponding multispecies cases.As a result, the escape rate dependence on the dipole field strength is weaker in the multifluid cases.The weak dipole field cases show 32% and 10% increases in the O 2 + and CO 2 + escape rates, respectively.The decreases are 35% for O 2 + and 80% for CO 2 + in the strong dipole field cases.The O + escape rate decreases with the dipole field strength, but the decrease is not as strong as that in the multispecies cases.  .The gray area is the range of the dipole field strength at the equatorial surface where its magnetic pressure becomes comparable to the solar wind dynamic pressure at different latitudes from the poles (lower limit: 1,200 nT) to the equator (upper limit: 2,400 nT).Note that the x-axis is scaled by the cubic root to show results clearly in the weak dipole field cases.
As with the density and flux distributions, the escape rates also exhibit asymmetry about the convection electric field of the solar wind in the multifluid cases.We divided the escape rate into two parts: one in the +E hemisphere (z > 0) and one in the − E hemisphere (z < 0).In all cases and all ion species (O + , O 2 + , and CO 2 + ), the escape rates in the +E hemisphere are a few times higher than those in the − E hemisphere.We also divided the escape rates in the multispecies cases, but the escape rates in both hemispheres are almost the same.It is consistent with the symmetric distributions of the densities and fluxes in the multispecies cases.
In Case MF1-4, no or weak dipole field cases, the O + escape rate in the +E hemisphere is nearly the same as the escape rate in the +E hemisphere of the corresponding multispecies case.In Case MF1, for example, the O + escape rate in the +E hemisphere is 2.36 × 10 27 s − 1 , accounting for 62% of the total O + escape rate.It is close to that in Case MS1, 2.34 × 10 27 s − 1 .Since the magnetosphere is not so expanded when the dipole field is absent or weak, ion pickup from the extended oxygen corona contributes to O + escape.Pickup ions float together with the solar wind flow in the multispecies simulation due to the single-fluid approximation.In the multifluid model, they are accelerated by the convection electric field and form a plume-like structure in the +E hemisphere.The velocity of such ion flows eventually reaches the solar wind velocity.Therefore, the escape rate in the +E hemisphere due to these processes is expected to be the same in both simulations.On the other hand, the O + escape rate in the − E hemisphere is lower in the multifluid cases than in the multispecies cases.A large part of pickup ions that originated in the − E hemisphere precipitates into the planet or escapes through the +E hemisphere, resulting in weaker ion pickup escape in the − E hemisphere.The results in Cases MF1-4 are consistent with this interpretation.This interpretation is invalid in the stronger dipole field cases because the contribution of ion pickup becomes small due to the expansion of the magnetosphere.Higher escape rates in the +E hemisphere are also seen in molecular ions.This is partly because the convection electric field can affect outflow from the ionosphere by penetration into the ionosphere, which enhances outflow in the +E hemisphere.In addition, the outflow can be accelerated by the convection electric field in a distant tail region.The escape rates in the +E (− E) hemisphere are increased (decreased) as they are calculated on a more distant sphere.The difference between the escape rate from each hemisphere is the smallest in Case MF6, the strongest dipole field case.The expansion of the magnetosphere reduces the penetration of the convection electric field and thus its effects on ion escape.

Discussions
The multispecies cases show consistent features with the previous study based on another multispecies MHD model (Sakata et al., 2022).However, the escape rates are different despite the same solar conditions and the neutral atmosphere.For example, the O + , O 2 + , and CO 2 + escape rates in the no dipole field case (Case MS1) are 4.68 × 10 27 , 3.61 × 10 25 , and 2.33 × 10 25 s − 1 , respectively, while those in the previous study (Case A1) are 5.17 × 10 27 , 1.84 × 10 25 , and 5.38 × 10 24 s − 1 , respectively.One reason for this discrepancy is the different treatments of the ion and electron temperatures.As coefficients of chemical reactions and collisions depend on the ion and electron temperatures, differences in the temperatures affect the photochemistry and dynamics of the ionosphere.Both models calculate the total thermal pressure P, which is the sum of the ion pressure n i kT i and the electron pressure n e kT e .In REPPU-Planets used in the previous study, the ion temperature T i and the electron temperature T e are assumed to be the same.In addition, the ion and electron temperatures have lower limits to reproduce temperature profiles in the lower ionosphere, which is consistent with the Viking observations.The detailed treatment of temperatures in REPPU-Planets is described in Terada, Kulikov, et al. (2009).On the other hand, MAESTRO solves the electron pressure equation and thus calculates the electron pressure and electron temperature independently.The ion pressure can be calculated by subtracting the electron pressure from the total thermal pressure in the multispecies cases.In addition, it includes photoelectron heating and inelastic collisional cooling, which are not included in REPPU-Planets.Another possible reason is the difference in the number of ion species.REPPU-Planets considers 11 ion species, while MAESTRO considers only five species and treats planetary and solar wind H + separately.The number of considered reactions is also different between the two models.These differences may affect the ion composition in the ionosphere.
The multifluid cases show significantly different features in the plasma distributions and the tailward fluxes from the multispecies cases.They are asymmetric with respect to the solar wind convection electric field.The multifluid model has separate momentum equations of ion fluids and allows planetary ion flows to have different velocities from the solar wind flow.In contrast, all ion species have the same bulk velocity, often dominated by the solar wind protons carrying large momentum, in the multispecies model.The separation of the momentum equations also allows the Lorentz force and the electric field to act on each ion fluid.Therefore, the distributions of the number density and the tailward fluxes show a plume-like structure in the +E hemisphere due to the acceleration by the convection electric field of the solar wind.It is a typical difference between multispecies and multifluid simulations indicated by previous studies (Ma et al., 2019;Najib et al., 2011;Regoli et al., 2018).
The multifluid representation also affects ion escape rates.In Case MF1 with no dipole field, the escape rates of O 2 + and CO 2 + are a factor of 3.9 and 2.6 higher than those in Case MS1, respectively, while the escape rate of O + is decreased by 19%.These effects on ion escape rates are also consistent with the previous studies.12 shows 1D profiles of the number density, the thermal pressure, and the vertical velocity of each ion species, as well as the magnetic field at 57.6°north latitude on the meridian plane, located in the region where the O 2 + downward flow occurs in Case MS6.The results of Cases MF6 (MS6) are plotted as solid (dashed) lines.For Case MS6, the total thermal pressure and the mass-averaged vertical velocity are plotted instead.The magnetic field of the two cases is nearly the same and almost radial, indicating that the profiles are located in the cusp region in both cases.The total thermal pressure profiles are also similar.The multifluid representation does not affect the global configuration of the pressure balance between the dipole field and the solar wind, as indicated in the total thermal pressure plots in Figure 5.Note that the thermal pressure is susceptible to a slight difference in the magnetic field because of low plasma beta.On the other hand, the density and velocity profiles are quite different.The solar wind proton shows strong downflow in Case MF6.Precipitation of solar wind protons is the dominant dynamics in the cusp.The precipitation of hot solar wind particles is also indicated by the dominance of the solar wind proton pressure in the thermal pressure profile in Figure 12b.The flows of planetary ions are weak and downward in the ionosphere and upward above ∼700 km altitude, except planetary H + with a strong downflow due to the momentum transfer from the solar wind H + .Case MS6, in which all ions have the same bulk velocity, shows the downward flow above an altitude of approximately 300 km.As the precipitation of solar wind proton is dominant, planetary ions cannot flow out through the cusp region.The different ion dynamics affect the structure of the ionosphere.The ionosphere is compressed above 300 km altitude due to the downflow in Case MS6.In contrast, the ionosphere inflates to higher altitudes due to a much weaker downflow in Case MF6.The number densities of molecular ions are more than two orders of magnitude higher at 700 km altitude, where the outflows emerge, resulting in strong upward fluxes of molecular ions shown in Figure 11.
This study assumed the steady solar XUV flux and solar wind conditions.However, they are constantly fluctuating, as represented by solar events such as flares and coronal mass ejection (CME) events.Disturbances of the upstream conditions and magnetosphere impact ion escape from planets.MAVEN observations indicated enhanced ion escape from Mars during CME events (Jakosky et al., 2015).At Earth, the O + outflow from the dayside cusp increases with the Kp index, an index of the geomagnetic field disturbance (Slapak et al., 2017).The IMF direction is also an important factor.A simulation study assuming weakly magnetized Mars pointed out that IMF rotation changes ion escape rates by more than one order of magnitude (Sakai et al., 2023).As observations of solar-like stars suggest that CME events occurred more frequently on the young Sun (Kay et al., 2019;Odert et al., 2017), it is expected that ion escape during CME events has played a significant role in the evolution of planetary atmospheres.
This study investigated the effects of a dipole field with a strength of up to 5,000 nT at the equatorial surface.This strength is lower than the estimated strength of the intrinsic magnetic field on ancient Mars, which can be comparable to the geomagnetic field (30,000-60,000 nT) (Weiss et al., 2008).The study of ion escape under the existence of such a strong intrinsic magnetic field would be valuable for evaluating the contribution of ion escape from ancient Mars.However, MHD simulation studies overlook some important processes at strongly magnetized planets.At Earth, ion outflow drivers are roughly categorized into suprathermal electrons and electromagnetic energy input from the magnetosphere (Strangeway, 2005).Suprathermal electrons enhance the ambipolar electric field by heating thermal electrons, which induces ion outflows.Electromagnetic processes are also important drivers.Heating by wave-particle interactions is a typical example.Ions perpendicularly heated by resonant interactions between particles and waves are accelerated upward by the magnetic mirror force.Such outflows through the cusps contribute to a large part of ion outflows at Earth (Yau & André, 1997).MAESTRO includes electron heating by photoelectrons and the ambipolar electric field by the electron pressure gradient to represent partly the effects of suprathermal electrons, but a large part of these processes is not included in the model.Recent global MHD models on Earth used ionospheric outflows derived from other models as the inner boundary condition.Brambles et al. (2010) used an empirical correlation between the O + cusp outflow and the Poynting flux by Strangeway (2005).However, it is questionable whether the empirical model on Earth can be applied to other planets, such as ancient Mars.Coupling an MHD model and a physical-based outflow model is more desirable for versatility.Glocer et al. (2009) coupled a multifluid MHD model with the Polar Wind Outflow Model that solves field-aligned transport equations from the ionosphere.Our model should be improved for simulations with a stronger dipole field to represent more realistic ion outflows.
The effects of an intrinsic magnetic field revealed in this study can be applied to other planets including exoplanets, to some extent.Ion escape plays a more important role on Earth-sized exoplanets because photochemical escape is inhibited by higher escape energy.In addition, potentially habitable exoplanets around K-and M-type stars are expected to be exposed to intense stellar XUV radiation and stellar wind due to the proximity to the central star.The existence of an intrinsic magnetic field and its effects on ion escape should be important from the perspective of atmospheric evolution and habitability on exoplanets.Although factors controlling ion escape, such as the XUV spectrum, stellar wind, planetary neutral atmosphere, and planetary magnetic field, are poorly constrained on most exoplanets, simulation studies on the effects of an intrinsic magnetic field on exoplanets may have implications for the habitability of exoplanets.The thermal pressure of five ion species, electron.The total pressure is the sum of ion and electron pressures.For Case MS6, only the total pressure is plotted.(c) The vertical velocity (V r ) of five ion species.For Case MS6, the bulk velocity is plotted.(d) The magnetic field in the spherical coordinate (B r , B θ , B ϕ ).

Conclusion
This study investigated the effects of an intrinsic magnetic field on ion escape from ancient Mars to understand the possible contribution of ion escape to atmospheric loss.As multispecies MHD simulations used in our previous studies (Sakata et al., 2020(Sakata et al., , 2022) ) may underestimate ion outflows from the ionosphere, we developed a new multifluid MHD model named MAESTRO.
MAESTRO is a global multifluid MHD model implemented with the cubed sphere grid, which consists of six identical faces with a quasi-uniform grid resolution.The model solves the continuity, momentum, and energy equations of each ion fluid, as well as the induction equation and the electron pressure equation.It also includes major photochemical processes and collisions for seamless plasma simulations from the ionosphere (100 km altitude) to the magnetosphere and the solar wind region (40 planetary radii).In this study, the model considers five ion species: planetary H + , O + , O 2 + , CO 2 + , and solar wind H + .The test simulation under the present-day Mars conditions reproduced the plasma boundaries, ionospheric profiles, and wake outflows consistent with previous observations and simulation studies, indicating that the new model has a good ability to simulate solar wind-Mars interactions.
We conducted multifluid MHD simulations with MAESTRO under the intense solar XUV and solar wind conditions at ancient Mars.Six simulation cases with different dipole field strengths from 0 to 5,000 nT at the equatorial surface were conducted to see the dependence on the dipole field strength.We also conducted multispecies MHD simulations under the same solar and dipole field conditions with the same model for comparison.
The magnetic configuration and the flux distributions change with the dipole field strength.The effects in the multispecies cases are consistent with the previous study (Sakata et al., 2020(Sakata et al., , 2022)).The multifluid cases show nearly the same magnetic configurations as the multispecies cases, but the plasma dynamics differ significantly.The multifluid representation enhances cold ion outflow, as indicated by stronger tailward fluxes and higher escape rates of molecular ions (O 2 + and CO 2

+
).The separation of momentum equations allows ionospheric ions to flow out from the ionosphere independently of precipitating solar wind protons.It also affects the plasma transport in the upper ionosphere, leading to different ion compositions and outflows.These effects are most remarkable in the multifluid case with the strongest dipole field, in which the escape rates of molecular ions (O 2 + and CO 2 + ) are two orders of magnitude higher than the corresponding multispecies case.The multifluid representation also allows the acceleration of planetary ions by the solar wind convection electric field.The multifluid cases show highly asymmetric plasma distributions characterized by a plume-like structure in the +E hemisphere.Ion pickup escape is suppressed in the − E hemisphere, which results in a slightly lower O + escape rate in the multifluid cases with no or a weak dipole field.In contrast, the O + escape rate is enhanced by one order of magnitude in the multifluid case with the strongest dipole field compared to the multispecies case because cold ion outflow becomes the dominant O + escape process instead of ion pickup.
When the dipole field is weak or absent, the total ion escape rates are comparable in the multifluid and multispecies cases and reach the order of 10 27 s − 1 , indicating that the multifluid representation does not influence the impact of an intrinsic magnetic field on atmospheric evolution.In contrast, the total escape rate is significantly different between the multifluid and multispecies cases when the dipole field is strong enough to sustain the solar wind dynamic pressure and form an Earth-like magnetosphere.Due to the multifluid representation, the total escape rate in the strongest dipole field case is enhanced by one order of magnitude compared with the multispecies cases.The existence of a strong dipole field can reduce the total escape rate by up to a factor of six, even in the multifluid simulations, suggesting that it could have influenced atmospheric evolution.

•
A new global multifluid magnetohydrodynamics (MHD) model of solar wind-Mars interactions with a cubed sphere grid system named Multifluid Atmospheric Escape Simulations Toward Real elucidatiOn is developed • The multispecies MHD cases underestimate outflows of ionospheric ions in strong dipole field cases where cusp outflow is important • The effects of the multifluid representation on the O + escape rate are small in no and weak dipole field cases where ion pickup is dominant Supporting Information: Supporting Information may be found in the online version of this article.

Figure 1 .
Figure 1.A cubed sphere grid used in the model with a grid resolution of N = 25.The coordinate lines of each face are also plotted.

Figure 2 .
Figure 2. The neutral atmosphere profiles used in this study.(a) Present-day Mars atmosphere based on the Viking measurements.(b) Ancient Mars atmosphere modeled by Kulikov et al. (2007).

Figure 3 .
Figure 3. 1D profiles in Case PM.(a) Pressure profiles along the subsolar line.The black line is the total thermal pressure, and the blue line is the total dynamic pressure.The green line is the sum of thermal, dynamic, and magnetic pressures.(b) Number density profiles at a solar zenith angle of 43.2°.Solid (dashed) lines are the profiles in the northern (southern) hemisphere on the x-z (meridian) plane.
Figure 6 shows the O 2 + number density in the same format as Figure 5.Note that the MSE coordinate coincides with the MSO coordinate in all simulation cases because the IMF is in the away sector of the Parker spiral.In the multispecies cases, the O 2 + number density is symmetrical between the +E and − E hemispheres.The dense population of O 2 + arises around the polar region and extends to the tail region in no or weak dipole field cases (Cases MS1 and MS3).The strong dipole field cases (Cases MS5 and MS6) do not show such dense populations, indicating suppression of ion outflow from the ionosphere.On the other hand, the plasma

Figure 4 .
Figure 4.The contour plots of the present Mars case on the x-z (meridian) plane.(a) O 2 + number density (n i ).(b) O 2+ tailward velocity (− V i,x ).(c) the O 2 + tailward flux (− n i V i,x).The value range of the colormap in each plot is the same as that in Figure5ofInui et al. (2019) to see variations in the tail region.Note that the maximum values of the contour plots are 1.08 × 10 5 cm − 3 , 174 km s − 1 , and 2.62 × 10 12 m − 2 s − 1 .

Figure 5 .
Figure 5. Contour plots of the total thermal pressure on the x-z (meridian) plane for Cases MF1/MS1, MF3/MS3, MF5/MS5, and MF6/MS6.The left column (a, c, e, g) is the multispecies magnetohydrodynamics (MHD) results, and the right column (b, d, f, h) is the multifluid MHD results.

Figure 6 .
Figure 6.Contour plots of the O 2 + number density on the x-z (meridian) plane for Cases MF1/MS1, MF3/MS3, MF5/MS5, and MF6/MS6.The format of the figure is the same as Figure 5.

Figure 7 .
Figure 7.The tailward fluxes of (a, b) O + , (c, d) O 2 + , (e, f) CO 2 + , and (g, h) the magnetic field line topologies at x = − 2 R p for Cases MF1 and MS1.Case MF1 (MS1) results are in the right (left) column.The white dashed lines in the tailward flux plots show the neutral sheet.The magnetic field line topologies are categorized by the topology: closed, open, and draped.They are also categorized by the sign of B x , that is, the direction of the magnetic field: positive (Sunward) and negative (Tailward).

Figure 8 .
Figure 8.The tailward fluxes of (a, b) O + , (c, d) O 2+ , (e, f) CO 2 + , and (g, h) the magnetic field line topologies at x = − 2 R p for Cases MF5 and MS5.The format is the same as Figure7.

Figure 9 .
Figure 9.The tailward fluxes of (a, b) O + , (c, d) O 2+ , (e, f) CO 2 + , and (g, h) the magnetic field line topologies at x = − 2 R p for Cases MF6 and MS6.The format is the same as Figure7.

Figure 11 .
Figure 11.Contour plots of the radial flux of O 2 + at 1,000 km altitude in Cases MF6 and MS6.Upward (downward) fluxes are plotted in red (blue).The dotted lines are the open-closed boundary, that is, the boundary between open and closed magnetic field lines.The subsolar point is located at (0, 0) in the longitude-latitude coordinate.

Figure 12 .
Figure12.1D profiles at 57.6°north latitude on the meridian plane of Case MF6 (solid) and Case MS6 (dashed).(a) The number density of five ion species.(b) The thermal pressure of five ion species, electron.The total pressure is the sum of ion and electron pressures.For Case MS6, only the total pressure is plotted.(c) The vertical velocity (V r ) of five ion species.For Case MS6, the bulk velocity is plotted.(d) The magnetic field in the spherical coordinate (B r , B θ , B ϕ ).

Table 1
Input Parameters for the Simulation Cases

Table 2
Escape Rates of O + , O 2 + , and CO 2 Regoli et al. (2018)compared multispecies and multifluid MHD simulations of present Mars with crustal magnetic fields.They showed higher ion escape rates in the multifluid simulation.The escape rates of O 2 + and CO 2 + are about four and five times higher, respectively, while the O + escape rate increased by only a factor of ∼1.5.Ma et al. (2019) conducted similar comparisons on present Mars.Their Case 3 (with no crustal field) showed that the escape rates of O 2 + and CO 2 + are two and three times higher, respectively, but the O + escape rate decreased by 10%.Large enhancement in molecular ion escape indicates that the multifluid representation significantly affects cold ion outflow from the ionosphere.In addition, Ma et al. (2019) studied the effects of an independent electron pressure equation that calculates the electron temperature and the electron pressure gradient force selfconsistently.The electron pressure gradient force affects ion dynamics via the ambipolar electric field, and the electron temperature affects some chemical reactions, such as the dissociative recombination of CO 2 + and O 2 + with electrons and electron impact ionization.They found that including the electron pressure equation increases the escape rate by up to a factor of two.Enhancement of molecular ion escape is more remarkable in stronger dipole field cases.Case MF6 with the strongest dipole field shows more than two orders of magnitude higher escape rates of molecular ions than the corresponding multispecies cases.The contour plots of the O 2 + radial flux in Case MS6 and MF6 in Figure 11 indicate how cold ion outflow is enhanced in the multifluid case.The dotted lines are the open-closed boundary, the boundary between open and closed magnetic field lines.In both cases, the upward O 2 + flow occurs at high latitudes where the magnetic field lines are open, indicating cold ion outflow through open field lines.However, the O 2 + flow is very weakly upward or rather downward around − 50°to 50°in latitude and mid-latitudes in the multispecies cases.Figure