IMF-induced escape of molecular ions from the Martian ionosphere

Since Mars does not possess a significant global intrinsic magnetic field, the solar wind interacts directly with the Martian ionosphere and can induce ion escapes from it. Phobos-2 and recent Mars Express (MEX) observations have shown that the escaping ions are O + as well as molecular O+2 and CO + 2 . While O + escape can be understood by the ion pick-up of non-thermal O corona extended around the planet, regarding the heavy molecular O +2 and CO + 2 , which are buried in the lower ionosphere, a novel escape mechanism needs to considered. Here we attack this problem by global magnetohydrodynamic (MHD) simulations. First, we clarify the global structure of the streamlines that result from the interaction with the solar wind. Then, by focusing on the streamlines that dip into the low-altitude part of the dayside ionosphere, we investigate the escape path of the molecular ions. The effects of the interplanetary magnetic field (IMF) on the molecular ion escape process are investigated by comparing the results with and without IMF. IMF has little effect on O escape via ion pick-up mediated by solar wind electron impact ionization of the O corona. O +2 and CO + 2 are shoveled from the low-altitude regions of the dayside ionosphere by magnetic tension in the presence of IMF. These ions are pulled by the U-shaped field lines to the north an south poles, and at the terminator, they are concentrated in the noon–midnight meridian plane. These ions remain confined to the noon–midnight plane as they are transported to the nightside to form the tail ray. Then they escape along the streamlines open to the interplanetary space. Under a typical solar wind and IMF condition expected at Mars, O , O+2 and CO + 2 escape fluxes are 8 .0× 10 23, 3.5× 1023 and 5.0×1022 ion s−1, respectively, which are in good agreement with the MEX observations.


Introduction
The first detailed measurements of the magnetic field at Mars by the Mars Global Surveyor (MGS) (Acuña et al., 1998(Acuña et al., , 1999) ) has established that Mars, like Venus, does not have a large-scale magnetic field.Although localized strong magnetic fields of crustal origin have been discovered, the absence of large-scale magnetic field implies that the solar wind interaction with the Martian ionosphere would be Venus-like.Indeed, magnetic structures, such as magnetic field rotations at the ionopause and magnetic field, suggesting Venus-like interaction have been reported (Cloutier et al., 1999).In 1989 the Phobos-2 mission discovered O + outflow from Mars.The ion loss routes are either ion pick-up by the solar wind or accelerated beam (up to several keV) from the ionosphere.The estimated O + outflow rate is ∼ 3.0×10 25 s −1 (Lundin et al., 1989).The resent observations by Mars Express (MEX), which was launched in 2003 and carried an ion mass analyzer (IMA) onboard, has shown that the molecular ions of O + 2 and CO + 2 are also escaping from the ionosphere.The loss rates are 1.6×10 23 s −1 , 1.5×10 23 s −1 and 8.0×10 22 s −1 , for O + , O + 2 and CO + 2 , respectively (Barabash et al., 2007a).The purpose of this study is to understand the escape processes of molecular O + 2 and CO + 2 ions from the Martian ionosphere.For this, we need to understand (1) convection in the dayside low-altitude ionosphere, (2) transport of ionospheric plasma from dayside to nightside and (3) the structure of the nightside ionosphere.We believe that a global three-dimensional MHD simulation including a realistic ionosphere model, and that with sufficient spatial resolution in the low-altitude part of the ionosphere, is one of the best approaches to the problem.Below, crucial elements that show up in our study are reviewed.
Because simultaneous observations of magnetic field and charged particles have not been made in the Martian Published by Copernicus Publications on behalf of the European Geosciences Union.
ionosphere, a number of theoretical models for the solarwind-ionosphere interaction of an unmagnetized body have been based on the observations at Venus.One of the prominent discoveries is the detection of large-scale magnetic field of solar wind origin in the dayside ionosphere well below the ionopause (Russell and Vaisberg, 1983).Ionopause is the plasma boundary between the solar wind and the ionospheric plasma.The interplanetary magnetic field (IMF) does not penetrate through this boundary in a "frozen-in" situation.The models developed to account for the formation mechanism of the large-scale magnetic field followed the "convection-diffusion" scenario.In this scenario, the key elements are the downward transport of magnetic flux from the ionopause region and the enhanced magnetic diffusion due to smaller conductivities in the lower ionosphere (Luhmann et al., 1984;Cravens et al., 1984;Phillips et al., 1984;Shinagawa et al., 1987;Shinagawa and Cravens, 1988;Shinagawa, 1996a,b).The model successfully reproduced the observed altitude profiles of the magnetic field that has a long enough decay time in order to explain the observations.However, they do not explain the reason why the downward transport is generated.Jin et al. (2008) discussed the relationship between the vertical convection and the solar wind dynamic pressure.Observationally, IMF penetration occurs when the solar wind dynamic pressure is relatively high.They suggested that the vertical convection is directed downward when an ion deficit occurs in the lower ionosphere.As the solar wind total pressure increases, the ionospheric pressure also increases in order to balance with the solar wind total pressure.The increase of the ionospheric pressure leads to an increased ionospheric electron density because the temperature at lower ionosphere does not almost change.The increase of the electron density leads to enhanced chemical loss of O + 2 by recombination.Thus an ion deficit occurs in the lower ionosphere, and a downward flow from the ionopause occurs to compensate for the ion deficit in the lower ionosphere.
As for the observations of the Martian ionosphere, the Viking 1 and 2 landers are the only spacecraft that observed ion profiles in the ionosphere (Hanson et al., 1977), but they did not measure the magnetic field.MGS explored magnetic field structures globally in the low-altitude region.However, there has not been a detailed observation of IMF penetration into the ionosphere, because simultaneous charged particle observations are absent.According to the Viking observation, the ionospheric pressure at Mars is found to be insufficient to balance the solar wind dynamic pressure (Hanson and Mantas, 1988).Therefore, the case of Mars resembles that of Venus to some extent.Theoretical models show that the magnetic field of the solar wind origin penetrates into the Martian ionosphere and sustains the solar wind dynamic pressure (Shinagawa and Cravens, 1989;Jin, 2004).
Regarding the transport of the ionospheric plasma from dayside to nightside, Shinagawa (1996a,b) performed twodimensional simulations for a Venus case.It was shown that the horizontal pressure gradient force is sufficient to produce the observed acceleration of ionospheric plasma from dayside to nightside, and the nightside ionosphere can be maintained by the plasma supply of this flow pattern.Because the simulated spatial range is 120-500 km, however, the upper part of the nightside region was out of the scope of the study.
The part of a field line embedded in the dayside stagnation region moves slowly, while those in the solar wind move at a rather fast speed.As a consequence, the magnetic field lines are stretched and draped around the planet.A threedimensional simulation does take into account the magnetic tension force arising from this field-line deformation that is crucial for understanding the nightside ionosphere structure.Tanaka and Murawski (1997) first calculated the nightside structure of the Venus ionosphere by a three-dimensional MHD simulation that covered the altitude range from 140 km to 10 R V (R V : Venus radius) including a realistic ionosphere model.This simulation suggests that the O + ions from the ionosphere are transported to a few R V altitude in the nightside, and then pulled anti-planetward by the magnetic tension.
Similarly, a Martian case has been simulated by Ma et al. (2004).The three-dimensional MHD simulation included a realistic ionosphere and the crustal field.The nightside structure of the Martian ionosphere, however, was not the focus of the study.
Various hybrid simulation studies have been made of the solar-wind-ionosphere interaction of unmagnetized planets (Modolo et al., 2006(Modolo et al., , 2005;;Bößwetter et al., 2004;Kallio andJanhunen, 2002, 2001;Shimazu, 2001Shimazu, , 1999;;Brecht, 1997;Brecht et al., 1993).The hybrid simulations treat ions as particles while electrons are considered as a massless charge-neutralizing fluid.The hybrid simulations do have the advantage of including ion kinetic effects.These simulations, however, did not incorporate realistic ionosphere models, and the grid resolution in the ionosphere has been relatively low; that is, the radial grid size is typically ∼ 100 km in the ionosphere.Using a hybrid model, it has been difficult to study the convection in the lower ionosphere, which is the crucial element for understanding the molecular ion loss from Mars, and where the importance of the ion kinetic effects is minimal because of the collisional nature of the local plasma.
These lead us to think that a global three-dimensional MHD simulation including a realistic ionosphere model, and that with sufficient spatial resolution in the lower ionosphere, is one of the best approaches to understand the molecular ion escape from the Martian ionosphere.Indeed Ma and Nagy (2007) calculated ion escape fluxes by using a multispecies three-dimensional MHD simulation code including a realistic ionosphere.This code included ionospheric processes such as ion production and loss due to photochemical reactions, gravity and collisions with neutral atmosphere, with a sufficiently small (10 km) radial grid spacing in the ionosphere.The estimated escape fluxes of O + , O + 2 and CO + 2 are 7.2×10 23 , 1.9×10 23 and 1.3×10 23 ion s −1 , respectively, under a solar minimum condition.While the estimated escape fluxes are consistent with the MEX observations, Ma and Nagy (2007) did not elucidate the ion escape mechanism.Najib et al. (2011) calculated ion escape fluxes by using a multifluid three-dimensional MHD simulation code that treats asymmetry due to the convection electric field included in the ionosphere process with a small (10 km) radial grid spacing in the same way as Ma and Nagy (2007).The estimated O + 2 escape flux is higher than that of Ma and Nagy (2007) compared to the escape rates of the other species.They discuss that this higher escape flux might be due to the new dynamics observed through their model as asymmetric plumes.However, Najib et al. (2011) did not elucidate the ion escape mechanism.
In this paper, we discuss the ion escape mechanism of the molecular ions such as O + 2 and CO + 2 .Ions of O + 2 and CO + 2 are mainly produced in the lower ionosphere.Therefore, in order to have these ionospheric molecular ions to escape by a large amount, the ions produced at such low altitudes need to be transported upward to reach the altitudes where ions find themselves on streamlines that are open to the interplanetary space.We will show that the presence of finite IMF and its magnetic tension in the dayside part is crucial for the escape of these two ion species from Mars.

Governing equations
A set of three-dimensional multispecies MHD equation is developed, which includes the continuity equations separately for H + , O + , O + 2 and CO + 2 ions, a common momentum equation, Maxwell's equations and the energy equation.These equations are summarized below: where ρ s is the mass density for the s ion component (s = H + , O + , O + 2 and CO + 2 ), ρ is the total mass density, u is the common velocity, B the magnetic field and p the total thermal pressure.γ is the ratio of specific heats (and taken to be 7/5).µ 0 is the magnetic permeability.
The source terms are S a = m s α s − β s ρ s , (5) where α s is the photochemical production rate per volume, m s is the mass per particle, β s is the chemical loss rate per particle for ion species s, G is the gravity constant, M is the mass of the planet, ν in is the ion-neutral collision frequency, D B is the magnetic diffusion coefficient, k B is the Boltzmann constant, and T 0 (h) is the assumed temperature of newly produced plasma as a function of altitude based on the Viking lander observations (Hanson and Mantas, 1988).This temperature model is shown in Fig. 1.N s is the number density of ion species s (N s = ρ s /m s ) and N e is the electron density that we assumed to be equal to the sum of N s .The ion-neutral collision frequency was set at where m e is the electron mass, ν en electron-neutral collision frequency, ν ei electronion collision frequency and e the electron charge.ν ei and ν en are given as ν ei = 54.5 × N e /T 3/2 e s −1 (T e , the electron temperature) and ν en = 3.68×10 −8 [CO 2 ] s −1 , respectively (both taken from Schunk and Nagy, 2000).We replaced T e with the plasma temperature T p (T p = p/N e /k B /2) since we did not separate ion and electron temperatures in the present model.

The neutral environment and photochemical reactions
Neutral atmospheric profiles are necessary for calculating ion production and loss rates due to photochemical reactions, and ion-neutral and electron-neutral collision frequencies.The neutral atmosphere has the thermal components and the nonthermal component.The non-thermal component is represented by oxygen atoms that are produced by the dissociative recombination of molecular ions O + 2 .The radial variation of the number density in the gravity field of the planet can be expressed as where n i (h) is the number density at the height h for the i neutral component (i = H, O, CO 2 ), T i is the neutral temperature.The equivalent temperature of non-thermal O is assumed to be 2900 K.The density of the non-thermal oxygen atom at the altitude (h 0 = 400 km) was set at 1.0×10 4 cm −3 .This radial profile is consistent with the profile of Kim et al. (1998) that was obtained by a Monte Carlo simulation.The thermal components are composed of H, O and CO 2 .The same constant temperature of 200 K is assumed for the altitude higher than 200 km.The neutral densities at the altitude h 0 = 200 km were set to 1.5 × 10 7 cm −3 for O, 3.0 × 10 7 cm −3 for CO 2 and 9.0 × 10 5 cm −3 for H.At altitudes lower than 200 km, the density profiles of H, O and CO 2 were taken from Krasnopolsky (1993) for a solar minimum condition.The altitude profiles of the temperatures model are from Krasnopolsky and Gladstone (1996).The neutral atmosphere model used in this study is shown in Fig. 1.
In the present model, photochemical reactions for four ion species (H + , O + , CO + 2 and O + 2 ) are taken into account.The considered photochemical reactions and their reaction coefficients are listed in Table 1.The photoionization frequencies at Mars were set at 1.0×10 −7 s −1 for O + and 2.6×10 −7 s −1 for CO + 2 under a solar minimum condition, which were taken from Torr and Torr (1985) for Earth's case and adjusted to the heliocentric distance of Mars.The absorption of EUV flux by the neutral atmosphere is included in the model.The photoionization frequencies are assumed to depend on solar zenith angle as ∝ cos(SZA).Reaction rates with solar wind protons and electrons need special care because they depend strongly on the temperature of incident ions and electrons.Reaction rates and average energy of ejected particles for the proton-neutral charge exchange (CE) and for the electron impact ionization (EII) reactions are obtained by Monte Carlo simulations, whose details are given in Jin et al. (2006).

Numerical method
A triangular grid has been generated on spherical surfaces.In the construction process of the grid system, it is desirable that two-dimensional spherical surfaces are covered by control volumes of similar size, because the integration time step is restricted by the smallest control volume.The number of grids on the spherical surfaces is 8000, which corresponds to 93 km × 93 km at the inner boundary.The altitude range Table 1.List of chemical reaction and rates considered in the model.Units are cm 3 s −1 .The reaction coefficients for the ion-neutral charge exchange (Reactions R1, R2 and R4) were taken from Anicich (1993), those for the recombination (Reactions R5 and R6) were taken from Schunk and Nagy (2000) and those for the protonneutral charge exchange (Reactions CE1 and CE2) and for the electron impact ionization (Reactions EII1 and EII2) were taken from Jin et al. (2006).

No
Reaction Rate coefficient of the simulated region is from 120 km (the bottom of the ionosphere) to 12 R M (R M : Mars radius).By changing bin size along the radial direction, Tanaka and Murawski (1997) numerically simulated the solar-wind-Venus interaction with sufficient resolutions for the entire region of interaction including the formation of the bow shock and detailed interaction processes occurring in the ionosphere.A similar idea was adopted in the present model in treating both the ionosphere and the boundary region with sufficient spatial resolutions ( r = 10-4967 km).To solve the MHD equations, the Lax-Wendroff scheme was used.This scheme has a secondorder accuracy both in time and space, but generates numerical oscillations near discontinuities and shocks.Therefore numerical diffusion is required to suppress the oscillations.
In the present calculation, numerical oscillations were kept at Ann. Geophys., 31, 1343-1356, 2013 www.ann-geophys.net/31/1343/2013/minimum by controlling the coefficients of numerical diffusion, which are of the third order as compared to the secondorder numerical accuracy of the Lax-Wendroff code.For the outer boundary, the variables were fixed at the solar wind values (as given in Sect.2.4) at the upstream boundary.Free boundary conditions are used for the downstream boundary.
For the inner boundary at the bottom of the ionosphere, the vertical plasma velocity was assumed to vanish, and the magnetic field magnitude was also assumed to vanish.

Solar wind parameters
For simulations with finite IMF, the solar wind parameters are as follows: plasma density of n sw = 3.0 cm −3 , velocity of V sw = 400 km s −1 , magnetic field of B sw = 4.0 nT and temperature of T sw = 1.0 × 10 5 K.The direction of IMF was assumed to be perpendicular to the solar wind flow.To clarify the role of IMF in the molecular ion escape process, comparison between the cases with and without IMF is made.In the case of no IMF, the solar wind parameters are chosen such that the total pressure, the velocity and the temperature in the post-shock region are the same as the corresponding value from the finite IMF case.Jin (2004) suggests that IMF penetration into the ionosphere occurs when the solar wind total pressure is relatively high.The present solar wind condition enables IMF to penetrate into the ionosphere.

Simulation results
The system was allowed to develop in time, and eventually a steady state was obtained.Figure 2 shows how IMF lines deform as they interact with the Martian ionosphere.The figure also shows the coordinate system: The solar wind flow is toward −x direction.The IMF is directed in the +y direction and the z axis completes the right-handed coordinate system.The x-y plane (the z = 0 plane) is called the equatorial plane, and the x-z plane (the y = 0 plane) is called the meridian plane hereafter.The field line in panel (a) is at z = 1 R M initially.With this shallow impact, the field line is only weakly bent as it is convected anti-sunward.On the other hand, the field line shown in panel (b), which initially has z = 0.1 R M , impacts deep into the ionosphere.Then the part of the field line embedded inside the ionosphere is slowed down so greatly that this part is delayed substantially from its ends in the solar wind.This results in significant bending of the field line -that is, the U-shaped field line shown in the panel.The field line exerts magnetic tension on the ionospheric plasma.This is one of the key elements of the present study.

The Martian nightside density structure
In our simulation, we obtain the densities of four ion species (H + , O + , O + 2 and CO + 2 ) by integrating the equation of continuity for each species.Each density distribution shows dif- ferent structures as shown in Fig. 3.The top half of each panel shows the meridian plane, while the bottom half shows the equatorial plane.We used a multispecies MHD model.Therefore our model results are plane symmetry; that is, the z > 0 result is a mirror image of the z < 0 result, and the y > 0 result is a mirror image of the y < 0 result.A ray structure can be recognized at the center of the tail.The ions in the tail ray are composed of, in orders of decreasing density, H + , O + 2 , O + and CO + 2 .Figure 4 shows the altitude profiles of the ion densities at the subsolar point.We compare our density distribution with Viking 1 lander observation in Sect.cm cm -3 -3 to high altitude on the nightside, O + 2 and CO + 2 need to be shoveled from the lower ionosphere in the dayside and accelerated into the tail ray.It is this mechanism that will be studied in this paper.
In order to understand the role of the solar wind magnetic field in the molecular ion transport, we show the density distribution of O + 2 and O + obtained in the case without IMF in Fig. 5. Comparison with Fig. 3 shows the following: (1) along the tail axis at 1 R M altitude, O + 2 density is less than 0.1 cm −3 for the case without IMF.The O + 2 density at the same location for the case with IMF is 10 cm −3 , implying 2 orders of magnitude reduction in the molecular ion density at the high-altitude position in the tail.That is, the ray structure of O + 2 forms only in the presence of IMF.
(2) O + 2 starts to be elevated to higher altitude already at the terminator for the case with IMF, while such a signature is absent in the case without IMF.This implies that only in the presence of IMF is O + 2 shoveled from the dayside lower ionosphere and transported to higher altitude in the meridian plane as they are convected tailward with the bent field line.(3) In contrast, regarding O + , the density pattern from the case without IMF shows a similar pattern to that on the equatorial plane from the case with IMF (Fig. 3c with the flow pattern.The flow pattern in the meridian plane is heavily affected by magnetic tension.O + 2 elevated to higher altitude upon crossing the terminator is shifted towards the equatorial plane in the nightside.The pair of flows from north and south collide on the tail axis at 0.5 R M above the surface.Then the converging flows are redirected to the anti-sunward direction.
One can also see that, in the equatorial plane, there is little connection between the low-altitude O + 2 in the nightside ionosphere and the O + 2 at higher altitude forming the tailward jetting ray.We conclude that the transport of the heavy ions in the meridian plane from dayside to nightside is crucial for forming the anti-sunward-flowing ray in the tail.Since the ions in the ray are eventually lost into the interplanetary space, the ray formation is the process that causes the heavy molecular ions to be lost from the planet.
What is missing in the above argument is the flow pattern in the dayside, which is now shown in Fig. 7. Color contours for V X , V Y and V Z at the ionospheric altitude of 400 km are shown.The tailward flow (V X < 0 and V Z > 0 in the Northern Hemisphere) converges towards the meridian plane (V Y > 0 on the Y < 0 side, and vice versa) at SZA > 75 • .This flow pattern is created as the ions are accelerated by the magnetic tension of a U-shaped field line whose apex is embedded in the ionosphere.In a sense, one may say that the escaping ionospheric ions are dragged by the solar wind ions.However, because the momentum exchange is mediated by the deformed field lines, the process involves the formation of the tail ray. Figure 8 shows the resultant O + 2 density concentration on the meridian plane at the terminator.The flow pattern of these concentrated O + 2 beyond the terminator is what has been discussed in the previous paragraph (Fig. 6).
In order to see where the O + 2 in the tail ray originates, we shown in Fig. 9  2 on the open stream lines are accelerated by the magnetic tension and escape to the interplanetary space.This O + 2 flux creates the molecular tail-ray structure.

Quantifying the ion escape flux
In MHD simulations, the ion escape rate is usually estimated by numerically integrating the escaping ion flux across a plane placed perpendicular to the solar wind at a reasonably large distance behind the planet (Ma and Nagy, 2007).This method is hereby called Method 1. We, however, found that this method tends to give a value larger than the one estimated independently from the chemical balance between the ion production and loss rates.The reason is the following: in MHD simulations, an artificial viscosity should be included to avoid numerical oscillations.The artificial relaxation tends to be large at a sharp boundary, one of which is the one between the stagnant ionospheric plasma and the flowing solar wind plasma.This region is the key to determining the molecular ion process.Indeed, the escape fluxes in the cases without IMF (Table 2) are ∼ 10 22 s −1 in the case where 10 km is the minimum grid size, and ∼ 10 23 s −1 in the case of 20 km.Therefore we decided to estimate the ion escape flux using another method as well.
Here we introduce the idea of calculating the ion escape flux that is less vulnerable to the artificial diffusion.This method, which is called Method 2 hereafter, is as follows.
In the steady state that we are dealing with, the ion escape flux F j can be also calculated by where V is volume, Q j is the ion production rate and L j is the chemical loss rate for ion species j .The crucial point here is that the integration should be made only over the volume through which the streamlines that are open to the interplanetary space pass.The ionospheric ions can escape only if they are situated on these streamlines.If the integration is made over the whole simulation box, we cannot distinguish the ion escape flux from the artificial diffusion flux.Method 2 is less damaged by the worst effect of the artificial viscosity that brings fluxes on a closed streamline onto an open streamline via diffusion in the direction perpendicular to the stream-lines.The interface between a closed and an open streamline tends to be associated with enhanced velocity shear and thus more vulnerable to the artifact.
Table 2 compares ion escape fluxes estimated by using Methods 1 and 2 and those obtained by MEX observations.In the case of no IMF, the escape fluxes of O + 2 and CO + 2 are negligible when calculated by Method 2, while they are estimated to be significant by Method 1.The values obtained by Method 2 have very small absolute values of less than 10 18 s −1 and a negative sign.We note that the escape flux calculated by Method 2 can become negative.There are ions that come onto the open streamlines by artificial diffusion and that are lost via chemical reactions.If this artificial loss is more than the integration of the source term along the open streamlines (which is ∼10 14 s −1 in this case), the flux is calculated to be negative.That is, the numerical diffusion tends to reduce the estimated escape flux by Method 2, while it tends to increase the estimated flux by Method 1.The true escape flux should be on the order of, or less than ∼ 10 14 s −1 , which is negligible as compared to the estimate by Method 1 (Table 2).On the other hand, in the case of simulation with a finite IMF, O + 2 and CO + 2 escape fluxes obtained independently by Methods 1 and 2 are in reasonably good agreement (the largest difference is the factor of 2 for the O + 2 flux).We consider that errors due to the artificial diffusion coming into these escape flux estimates are within the acceptable range.Now we discuss the ion loss rate obtained by Method 2. In the presence of IMF, it is found that O + 2 and CO + 2 escape flux are 3.5×10 23 s −1 and 5.0×10 22 s −1 , respectively.These values are in relatively good agreement with the MEX observation.The molecular ion loss rates are negligible in the absence of IMF, which indicates that the shoveling from the dayside lower ionosphere and acceleration into the tail ray require the magnetic tension.The escape flux of O + is rather insensitive to IMF.This fact is consistent with the interpretation that most of the escaping O + are produced by ionizations of the non-thermal O corona that is extended to very high altitudes.Note that the ionized O + is not accelerated by the convection electric field in the case without IMF.However, 62 % of non-thermal O made by the dissociative recombination of O + 2 has more than 2 eV, which is the energy of escape velocity for oxygen (Nagy and Cravens, 1988).Therefore we consider that O + has more than 2 eV after non-thermal O is ionized.O + can escape because O + has more than 2 eV.The magnetic tension effects are crucial for the dynamics in the dayside lower ionosphere and in lifting the streamlines into the tail ray.While most of the molecular ions escape along this path, O + are already on the solar-wind-like streamlines when they are ionized.Indeed, a step in Method 2 allows us to confirm that the O + escape rate is sufficiently provided by

Molecular ion production on open streamlines
From the arguments deployed above, it is now clear that, for the molecular O + 2 to escape, the crucial issue is how to situate the molecular ion production region (Q O + streamlines.In order to clarify this better, we plot in Fig. 10 the color contour of the O + 2 production/loss rates.The horizontal axis is SZA, while the vertical axis is the altitude.The three panels show no IMF case (top), the equatorial plane for the case with IMF (middle) and the meridian plane for the case with IMF (bottom), respectively.
As shown in the top panel, the O + 2 production region is located at the altitude of 200 km in the dayside.The produced O + 2 are transported downward and are subject to the loss process at the altitude of 150 km, where the recombination reaction of Reaction (R6) in Table 1 occurs (Jin et al., 2008).In the nightside, the O + 2 production region is located in a thin region at the altitude 200 km.The production is due to the precipitation of O + from higher altitude.O + impact to neutral CO 2 at the altitude of 200 km and O + 2 is generated by the Reaction (R2).Also shown in Fig. 10 are the flow vector (white) and the open/closed demarcation line.In the absence of IMF, the open/closed boundary is located above 400 km altitude -that is, far above the altitude of the O + 2 production region (top).We have seen that O + 2 escapes little in such cases because the production region is not embedded in the region of "open" streamlines.
When IMF is present, the field lines are deformed to exert tension force on the ionospheric plasma.In particular, magnetic tension changes the flow direction on the nightside from planetward to anti-planetward.This shifts the open/closed boundary downward in the equatorial plane as shown in the middle panel.This, however, is not enough to make it touch the O + 2 production region.Note that the magnetic tension effects are smallest on the equatorial plane.
In • is on open streamlines.This indicates that the CO + 2 escape mechanism is the same as that for O + 2 .The production rate of CO + 2 , however, is less, which leads to the smaller escape flux of CO + 2 shown in Table 2.

Discussion and summary
We have shown that O + 2 and CO + 2 can escape from the Martian ionosphere in the presence of finite IMF because O + 2 and CO + 2 are shoveled from low-altitude regions of the dayside ionosphere by the magnetic tension.Without IMF, O + 2 and CO + 2 transported from the dayside form a vortex in the nightside and return to the ionosphere.With a finite IMF, O + 2 and CO + 2 ions transported to the nightside form a ray structure at the center of the tail.They are pulled anti-planetward by the U-shaped field lines, and the ions do not return to the Martian ionosphere.For this process to be effective, IMF-lines need to penetrate into the low-altitude ionosphere where the molecular ions are produced.It is noted that, for the typical solar wind parameter adopted in our simulation, IMF penetrates into the low-altitude ionosphere to shovel the molecular ions.
We have shown that the spatial distribution of plasma in the tail puts on a ray structure at the center.The distribution of ion fluxes observed by MEX at the altitude of 1 R M on the nightside is shown by Barabash et al. (2007a).They show that the ions including the molecular ions are concentrated at the center of the tail and create a sheet-like structure corresponding to the current sheet in the induced magnetosphere tail.We can find this sheet-like structure in the meridian plane in our simulation.The same tail structure in Venus was observed by Venus Express (VEX) (Fedorov et al., 2011).The observed distribution of ion fluxes is consistent with our simulation result and would signify the effects of the magnetic tension.However, we should point out differences between our simulation result and the observations.The observations suggest that the observed heavy ion escape fluxes are higher in the positive-z hemisphere, which is the direction of the solar wind convection electric field, than in the negative-z hemisphere (Barabash et al., 2007a).Because our simulation is a multispecies MHD code, it does not include the effect of ion-finite gyro radius and asymmetry due to the convection electric field.In our simulation, the distribution of heavy ions is symmetrical to the z direction.The observed asymmetry would be understood by taking into account the finite-radius effect under the z-directed solar wind convection electric field by using multifluid MHD simulation or hybrid simulation, which is out of the scope of the present study.Another difference between our simulation result and the observations is low-energy ions distributed as a ring (like a smile) in the Southern Hemisphere observed by MEX and VEX (Barabash et al., 2007a,b).These ions are observed on the field lines instantly connected to the dayside ionosphere.This escape process does not occur in this paper's results.However, in the case of low solar wind dynamic pressure in our simulation, heavy ions at low energy escape on the field lines in the vicinity of equatorial plane because the convection in the dayside ionosphere is upward and the heavy ions that flow at high-altitude are transported from dayside to nightside.This escape process may indicate the observations of low-energy ions distributed as a ring.Fränz et al. (2010) suggests that a ratio of O + /O + 2 is almost 1 at an altitude of 290-500 km at the terminator.In our evaluated escape fluxes, the ratio is ∼ 2.3 derived from Table 2.The difference of the ratio is due to difference of the evaluated region.The non-thermal O is distributed at an altitude of a few R M and escapes by pick-up process.Our escape flux includes the ion escape at high altitude because our evaluated region is all simulation boxes.In our simulation, the O + /O + 2 ratio at altitude 290-500 km at the terminator is almost 1.This is consistent with their result.
We compared heavy ions escape fluxes in our simulation with the observation of Barabash et al. (2007a) as shown in Table 2.The escape fluxes observed at tail region are consistent with our simulation results.However, the escape fluxes are estimated in the energy range to be more than 30 eV.Because MHD escape rate contains the whole energy range, a realistic MHD escape rate should be larger than the measured values.It is important that we mention this inconsistency.Our MHD model cannot treat the non-MHD escape process such as that of sputtering.We need to consider the non-MHD escape process in order to compare the observed escape flux more carefully.In addition, our simulation model does not include crustal field.The crustal magnetic field can protect the low-altitude ionosphere from IMF.This should reduce the total escape rate provided by the presented mechanism.Lundin et al. (2009) estimated the escape flux including the low-energy escape fluxes to be less than 50 eV at the terminator region.Those escape fluxes are O + = 2.1 × 10 24 , O + 2 = 1.4×10 24 , CO + 2 = 3.5×10 23 s −1 .These escape fluxes are larger than our simulation results.Reasons for these differences could be that (1) these escape fluxes were evaluated at the terminator region supposing that these escape fluxes at the terminator were x-axial symmetry and that (2) a part of observed flux could not escape but could return to the nightside ionosphere.Consequently, these values are larger than our simulation results.
We should discuss the limitations of the used multispecies MHD model.In our simulation, all ion species have the same bulk velocity.This supposition is not useful in the case of a high-temperature region such as the solar wind, but it is useful in the case of a low-temperature region such as the ionosphere because bulk velocity of pick-up ions become the background fluid velocity and the gyro radius is proportional to the square root of temperature.Molecular ions are produced at low altitude in the ionosphere.Therefore we could estimate these escape fluxes by integrating production and loss on the open streamline in the ionosphere.However, the molecular ions produced at low altitude are transported to high altitude.At high altitude, the gyro radius becomes longer, and the escape path of these ions would cause asymmetry as mentioned above.Moreover, a part of an ion picked up at high altitude such as O + may return to Mars because gyro radius at the solar wind region is a few times the Martian radius.We do not consider the effect of this return.
We have shown the ionosphere density distribution at subsolar point.We should compare our density distribution with the Viking 1 lander observation (Hanson et al., 1977).Den-sity distribution below 250 km is the same as photochemical equivalence, while density distribution above 250 km depends on the solar wind dynamic pressure.Therefore it is difficult to compare our result with the Viking 1 lander observation at SZA = 45 • directly.Shinagawa and Cravens (1989) and Fox and Hać (2009) simulate the altitude profiles of the ion densities using one-dimensional simulation, and compare with the Viking observation.Shinagawa and Cravens (1989) point out that the ionosphere should have downward velocity and horizontal ion loss due to large-scale horizontal plasma convection in order to explain the observed ion profiles such as O + 2 and O + .However, they suggest that the horizontal loss rates are artificially reduced in their model because the horizontal loss is strong when the ionosphere has the downward velocity.Fox and Hać (2009) show that their density distribution is consistent with the Viking observation when vertical velocity in the ionosphere is upward.These simulation results indicate that it needs the ion loss due to horizontal or upward convection in the ionosphere to explain the observation profile.In our simulation under the used solar wind parameters, the density profiles at SZA = 60 • in the meridian plane is in good agreement with the Viking 1 lander observation because ions such as O + 2 and O + are lost due to the horizontal and upward convection produced by magnetic tension.The density profiles at SZA = 45 • , where the vertical velocity is downward, in our simulation are smaller than the densities of the Viking observation because the horizontal ion loss is strong.In order to reproduce the Viking density profiles at SZA = 45 • in our simulation, we need to calculate a case of low solar wind dynamic pressure in which the vertical convection is upward in the ionosphere.
We have also studied how the escape fluxes of different ion species depend on the solar wind dynamic pressure.Two additional cases with IMF, in which the solar wind density and thus the dynamic pressure are switched to 0.1 and 10 times the value used in the run shown above, have been performed.The escape flux of O + is found to be in direct proportion to the solar wind density.This is understandable if we note that the rate of electron impact ionization of O corona is proportional to the solar wind electron density.What is interesting is that the escape flux of the molecular ions (O + 2 and CO + 2 ) has a weaker dependence on the solar wind pressure than O + .The escape flux of O + 2 varies from 2.0×10 23 to 9.0×10 23 s −1 when the solar wind density is changed by 2 orders of magnitude.The molecular ions are produced in the low-altitude region where CO 2 is abundant.When the solar wind density is elevated by a factor of 10, since the open/closed boundary is already located at below 200 km altitude, where strong collisions with neutrals exist, it is not possible to lower the boundary substantially.The escape flux is only slightly increased.When the solar wind density is reduced to 1/10, under the low dynamic pressure, the convection in the dayside ionosphere becomes upward (Jin, 2004) ionosphere along the open streamlines.However, since the ions are produced deep in the ionosphere and the upward speed is low (less than 100 m s −1 ) at low altitudes, the ions are quickly lost by recombination before reaching higher altitudes.As a result, the escape flux is slightly reduced.We compare the obtained solar wind dependence of ion escape fluxes with observations.Lundin et al. (2008) suggested that there is a strong correlation between the solar wind dynamic pressure and the planetary ion loss rate.They focused on the heavy ions in the energy range of 30-800 eV, and reported that the total escape flux increases from 6.0 × 10 23 to 6.0 × 10 24 s −1 as the solar wind pressure increases from 0.1 to 10 times the nominal value.However, our simulation result indicates that the total (O + , O + 2 and CO + 2 ) escape flux increases from 2.0×10 23 to 9.0×10 24 s −1 for the same range of solar wind pressure variation.The reason why the observed variation is smaller may be that the upper limit of the energy range dealt in the observation is not high enough to include those O + that are picked up at relatively high altitudes where the flow speed is high.When the O + flux is not fully included, the contribution of O + 2 and CO + 2 , which have a weaker dependence on the solar wind pressure, to the total ion escape flux becomes more substantial.
Here let us discuss the molecular ion (O + 2 ) escape for the Venus case.When the solar wind dynamic pressure is high, IMF penetrates into the dayside ionosphere and the molecular ions' escape can be expected.In order to simulate the Venus case, we change the neutral atmospheric model from Mars to Venus.The solar wind parameter and photoionization frequencies are arranged appropriately as well.In the results, we have found that the escape flux of O + 2 is estimated to be ∼ 10 23 s −1 when the solar wind dynamic pressure is high; that is, when IMF penetrates into the dayside ionosphere.Meanwhile, the O + escape flux for this case is ∼ 10 25 s −1 , which is larger by an order of magnitude than the Martian case.The reason for this enhanced O + flux is simply that the solar wind electron density at Venus is larger than at Mars.Barabash et al. (2007b) show that the escape of O + from Venus is observed, while molecular ion (such as O + 2 and CO + 2 ) escape is not observed.Our simulation results show that the O + 2 and CO + 2 do escape at the same rate as for the Martian case when the solar wind dynamic pressure is high.The molecular ion escape, however, would not have been detected by VEX observation because it would have been masked by the far larger flux of O + .

Fig. 1 .
Fig. 1.The neutral atmosphere profiles (left) and the temperature profile (right) used in the present model.

Fig. 2 .
Fig. 2. The distribution of O + 2 density is shown by color contours, and the three-dimensional configuration of magnetic field lines is shown by white lines.The direction of the solar wind flow is toward −x direction.The IMF is directed in the +y direction and the z axis completes the right-handed coordinate system.The field line in panel (a) is at z = 1 R M initially.The field line in panel (b) is at z = 0.1 R M initially.
4. There are peaks of H + and O + at ∼ 240 km altitude.The peak of H + is created by the chemical reaction in the ionosphere.The heavy ions such as O + 2 and CO + 2 are produced at the low-altitude region, and the densities are peaked at ∼ 120 km altitude.While the www.ann-geophys.net/31/1343

Fig. 3 .Fig. 4 .
Fig. 3.The distribution of H + (top left), O + (top right), O + 2 (lower left) and CO + 2 (lower right) densities shown by color contours.The solar wind flows from the right.The top half of each panel represents the meridian plane, while the bottom half represents the equatorial plane.

Fig. 5 .
Fig. 5.The distribution of O + 2 and O + density for the case without IMF.The vertical axis is y 2 + z 2 because of x-axial symmetry.
the distribution of O + 2 flux (color) at the ionospheric altitudes.The horizontal axis is SZA and the vertical axis is the altitude.Here we focus on the altitude range from 100 to 850 km.The top panel shows the meridian plane.The bottom panel shows the equatorial plane.The black lines in the top and the bottom panels are the demarcation lines between the open and the closed streamlines.Here we define the topology of stream lines as follows.Stream lines that are connected to infinity of downstream are defined to be "open".Thus particles forming O + 2 flux on the open region escape to the interplanetary space.The lines that have ends in the ionosphere are defined to be "closed".Particles forming O + 2 flux on the close region do not escape, because the stream lines are connected to the ionosphere.Of course, we assume here that the simulation has already reached a quasi-steady state.Comparison of two panels shows that the O + 2 flux in the x-z (meridian) plane is much larger than the x-y (equatorial) plane.O + 2 is produced at a low altitude of 200 km.In the x-y plane, O + 2 production region is not on the open region but is connected to the lower ionosphere.In the x-z plane, the O + 2 production region at 200 km is in the open region.Produced O +

Fig. 6 .Fig. 7 .
Fig. 6.The distribution of O + 2 density (color) superposed by the velocity vectors (white arrows).The solar wind flows from the right.The left panel represents the meridian plane, while the right panel represents the equatorial plane.

Fig. 9 .
Fig. 9. Contour plots of the O + 2 flux distribution.The horizontal axis is in the solar zenith angle and the vertical axis is the altitude.Here we are focusing into the ionospheric altitude.The flow patterns are superposed by the white arrows.The top panel shows the meridian plane, and the bottom panel shows the equatorial plane.The label "open" ("closed") indicates the region that is filled by open (closed) streamlines.The black lines in the top and the bottom panels are the demarcation lines between these two regions.

Fig. 10 .Fig. 11 .
Fig. 10.Contour plots of the net O + 2 production/loss obtained by our MHD simulation.The arrows on the panels indicate velocity vectors.The top panel corresponds to the values in the case of no IMF.The middle panel corresponds to the values in the equatorial plane containing IMF.The bottom panel corresponds to the values in the meridian plane affected by magnetic tension.An enlargement of the white rectangular area in the bottom panel is shown.These are shown as functions of the solar zenith angle (horizontal axis) and the altitude (vertical axis).The label "open" indicates the region in which the open stream lines exist.The label "closed" indicates the region in which the closed stream lines exist.The red arrows indicate the characteristic flow direction.The flow direction on the nightside is planetward in the case without IMF, while the flow direction on the nightside is anti-planetward in the case with IMF.The flow direction on the dayside is downward in both cases; that is, the solar wind penetrates into the ionosphere.