Simulated solar wind plasma interaction with the Martian exosphere: influence of the solar EUV flux on the bow shock and the magnetic pile-up boundary

The solar wind plasma interaction with the Martian exosphere is investigated by means of 3-D multi-species hybrid simulations. The influence of the solar EUV flux on the bow shock and the magnetic pile-up boundary is examined by comparing two simulations describing the two extreme states of the solar cycle. The hybrid formalism allows a kinetic description of each ions species and a fluid description of electrons. The ionization processes (photoionization, electron impact and charge exchange) are included self-consistently in the model where the production rate is computed locally, separately for each ionization act and for each neutral species. The results of simulations are in a reasonable agreement with the observations made by Phobos 2 and Mars Global Surveyor spacecraft. The position of the bow shock and the magnetic pile-up boundary is weakly dependent of the solar EUV flux. The motional electric field creates strong asymmetries for the two plasma boundaries.


Introduction
Due to the absence of a global magnetic field of the internal origin (Acuña et al., 1998), the upper atmosphere and the dayside ionosphere interact directly with the solar wind plasma.The Magnetometer/Electron reflectometer (MAG/ER) experiment on board Mars Global Surveyor (MGS) spacecraft has revealed the presence of crustal magnetic fields (Connerney et al., 2001).
Correspondence to: R. Modolo (ronan.modolo@cetp.ipsl.fr)The total magnetic field of the crustal sources exceeds, in few regions, 1500 nT at ∼100 km altitude (Acuña et al., 1999) but is also highly localized (Acuña et al., 2001;Connerney et al., 2001) and thus only weakly influences the global interaction between the solar wind and the planet.Therefore this interaction is similar in many ways to that of Venus or comets (Acuña et al., 1998).The MGS measurements have provided us a clear identification of the boundaries such as the bow shock (BS) and the magnetic pile-up boundary (MPB), located downstream the bow shock and upstream the ionopause.A new signature of this plasma boundary, namely, the enhancement of the magnetic field draping, have been recently emphasized both for Mars (Bertucci et al., 2003a) and Venus (Bertucci et al., 2003b).An extensive study of these two main boundaries have been reviewed by Mazelle et al. (2004) and Nagy et al. (2004).
Numerical simulations, with a different approach, have been performed during the last ten years to characterize the solar wind interaction with the exosphere of the planet (Brecht and Ferrante, 1993;Kallio and Janhunen, 2001;Ma et al., 2004).This interaction is strong enough to modify the Martian atmospheric escape and might has influenced the chemical evolution of the planet.
Single-fluid MHD models reproduce successfully some globals properties of the Martian environment such as the BS or the ionospheric plasma flow (Liu et al., 1999) and can provide meaningful results in the ionosphere where the spatial resolution can reach about 50 km (Ma et al., 2004) but failed to reproduce some important observed signatures.Two ions populations (light and heavy ions), at least, appear to be essential in the formation of the MPB and thus bi-fluid MHD model (Sauer et al., 1994) or hybrid models are mandatory to describe fully this boundary.
Published by Copernicus GmbH on behalf of the European Geosciences Union.R. Modolo et al.: The Martian plasma boundaries This study is based on the development of an hybrid model which includes a kinetic description for ions, that allows a full description of the dynamics of each ion species.In this model a self consistent treatment of the ionization processes is involved and the ionization rates are not fixed but calculated locally and for each planetary ion species using the values of a neutral density, the ionization frequencies of the electron impact, photoionization and cross sections for charge exchange reactions.The planetary plasma is greatly modified by the solar EUV flux which changes the photoionization frequencies but also governs the extension of the neutral corona.We focus this study on the influence of the solar EUV flux on the location of two main plasma boundaries: the bow shock and the magnetic pile-up boundary.This work continues the previous paper (Modolo et al., 2005).
The paper is organized as follows: a brief description of the simulation model and the parameters used are presented in Sect.2, while Sect. 3 describes the influence of the solar EUV flux on the position and the shape of the bow shock (Sect.3.1) and the magnetic pile-up boundary (Sect.3.2).Results are compared to the observations made by the Phobos 2 and Mars Global Surveyor spacecraft.

Simulation model
Inertial lengths and gyrofrequencies of ions and electrons are the characteristic scales to describe the macro and microphysics of the interaction processes.So for instance, waves have been identifed upstream of the Martian BS at the local proton gyrofrequency and is associated with the pickup Mars neutral hydrogen exosphere (Brain et al., 2002), the thickness of the shock overshoot is about few ion inertial length (Tátrallyay et al., 1997), while electrons seem to play a major role in the Martian aurorae recently observed with the Mars-Express spacecraft (Bertaux et al., 2005).In order to describe completely the plasma environment it is necessary to resolve the smaller length and time scales.Due to the computer limitations, a resolving on electron inertial lengths is not possible yet in global simulations.Simplifications and assumptions must be introduce to model the interaction.Since the plasma is governed by the ion dynamic and the typical length scale of the system is about few ion inertial length, a fully description of the ion species is essential but electrons can be characterized by a fluid, assuming that the spatial and temporal variations are greater than the typical electronic scales, this assumption defines the hybrid models.
A 3-D and multi-species hybrid model have been used to investigate the solar wind interaction with the Martian neutral environment.Ions are characterized by a gathering of macroparticles, describing the full dynamics of each ions species, while electrons are described by a massless fluid ensuring the conservation of the charge neutrality of the plasma, and contributing to the electronic pressure and the electric currents.The electronic temperature is determined by the equation of state P e =n e k B T e and the electronic pressure by assuming that the motion of the electron fluid is adiabatic with a polytropic index equal to 2. Solving the Faraday's law yields a time evolution of the magnetic field.The simulation model is based on the algorithm developed and tested by Matthews (1994).
The coordinate system used is such that the X axis points away from the Sun (X=V sw / V sw ), where V sw is the solar wind velocity, the Y axis is toward the motional electric field, defined by E conv = −V sw ×B IMF where B IMF is the Interplanetary Magnetic Field (IMF), and Z completes the right-handed system.The angle between V sw and B IMF is a free parameter taken to 90 • for this study.Simulations are computed on a 3-D uniform grid with a spatial resolution of 300 km (equals to twice the inertial length of the solar wind protons).The dimension of the simulation box is defined by −2.85≤X≤+2.85Martian radii (R M ) and −6.6≤Y, Z≤+6.6 R M .The model includes a fully absorbing obstacle with a radius of 3400 km (ions which penetrate into the obstacle boundary are stopped).No crustal fields are presently taken into account and no particular condition on the electromagnetic field is imposed inside the obstacle.Open boundaries are used in the solar wind direction and periodic boundaries in the perpendicular directions (except for planetary ions which are escaping freely from the simulation domain).
The solar wind is described by two ion species (H + and He ++ ions) while H + and O + ions represent the planetary plasma.Although alpha particles are a small fraction of the solar wind ions they contribute up to 20% of the solar wind momentum and can influence the overall dynamics of the plasma.In order to avoid of statistical problems the same number of macroparticles is used to describe protons and alpha particles, and we adjust the weight allocated to each macro-particle to preserve the physical density and the properties of the plasma.
The number of macroparticles per cell representing protons and alpha particles is equal to 2 for each species.The proton density is taken to 2.3 cm −3 and alpha particles contribute 5% of the solar wind density.The undisturbed velocity distribution function of the both ion species is assumed to be Maxwellian with a bulk velocity of 400 km s −1 and temperatures are fixed to 5 eV for protons and 22 eV for alpha particles.We take the interplanetary magnetic field to 3 nT.With such parameters the Alfvén and the sound velocity are close to 40 km s −1 ; correspondingly, the Alfvén and the sound Mach number of the solar wind are almost equal to 10.
The Martian neutral environment is described by two coronae of atomic oxygen and hydrogen with a spherical symmetry.The density profile of atomic hydrogen is given by the simplified Chamberlain's expression.The density and the temperature at the exobase are the key parameters which govern the density profile and we have adopted the values given by Krasnopolsky (1993a,b).Oxygen profile density is derived from Monte Carlo simulation computed by Kim et al. (1998), for the supra-thermal oxygen population, and from the atmospheric model for cold oxygen atoms (Krasnopolsky and Gladstone, 1996).The neutrals are ionized by solar EUV, electron impact and charge exchange processes.A more detailed description of the model as well as density profiles for the hydrogen and oxygen coronae are given in (Modolo et al., 2005).The solar EUV flux does not only change the photoionisation frequency, but also influences the extension of the coronae.To investigate the effect of this factor on the Martian environment, the simulations carried out for two extreme values of the solar cycle activity are compared.

Simulation results
The influence of the solar EUV flux on the plasma environment of Mars has been discussed in the previous paper (Modolo et al., 2005).We have shown that the solar activity strongly affects the Martian environment and the loss of planetary ions, contributing to the atmospheric erosion.The effect of the solar EUV flux on the location and shape of two main boundaries: the bow shock (BS) and the magnetic pile-up boundary (MPB) is now investigated.The simulation program is run from time t 0 s to t 1000 s and a quasi-stationnary solution is reached around time t 450 s.The simulation results are presented at time t 520 s and are similar to snapshots of the state of the simulation.
In order to discuss the main features of the plasma boundaries, we examine the results of 2-D cuts from the simulation domain: the XY plane containing the solar wind bulk velocity V SW and the motional electric field E conv , and the XZ plane containing V SW and the interplanetary magnetic field B IMF .Figures 1a and 1b display the magnetic field intensity, for a period of minimum solar EUV flux, in the XY and XZ planes, respectively.Magnetic field lines pile-up against the obstacle on the dayside and drape around the planet producing two magnetic lobes at the nightside.Figures 1c and 1d present the total magnetic field, for a maximum solar EUV flux, in the XY and XZ plane.
Two distinct boundaries are observed in Figs.1: the bow shock, separating the solar wind and the magnetosheath, and characterized by a strong jump of the magnetic field value, and the magnetic pile-up boundary, where the a second strong jump of the magnetic field in front of the obstacle is noticed.The magnetic pile-up boundary separates the magnetosheath from the magnetic pile-up region (MPR), dominated by planetary ions.The solid and dashed black curves depict the fits to the MGS observations for BS and MPB (Vignes et al., 2000).Since Vignes et al. (2000) have assumed a cylindrical symmetry of the boundaries, we use the same fits for the both planes.
The magnetic field strength jumps by a factor 3 at the subsolar shock and climbs sharply to about of 15 nT in the magnetic barrier above the obstacle.Due to a numerical diffusion  (Trotignon et al., 2006).and a finite grid size, the simulated magnetic field intensity is lower than in the observations where the magnetic field magnitude can reach 25-30 nT in the MPR bounded by the MPB (Bertucci et al., 2003a).The numerical diffusion allows to the magnetic field lines to penetrate slightly inside the obstacle.Future simulations, with a better spatial resolution and including the crustal magnetic field, should improve this issue.By comparing the XZ maps of the total magnetic field in minimum and maximum solar condition (Figs. 1a,1b and 1c,1d), no clear differences have been found, suggesting that the global magnetic environment in the vicinity of Mars is weakly dependent of the solar EUV flux.

The bow shock
A strong asymmetry in the magnetic field intensity is observed in the XY plane (Fig. 1a), while the field distribution in the XZ plane (Fig. 1b) |B| is almost symmetric.This asymmetry is attributed to the action of the motional electric field on the planetary ions which gain the momentum mainly in the positive Y direction, i.e. the motional electric field convection, breaks the symmetry.In these simulations, the shock is found to be farther from the planet in the opposite direction of the motional electric field that is in agreement with the previous hybrid simulations computed for Mars and Venus (Moore et al., 1991;Brecht and Ferrante, 1993;Shimazu, 2001;Kallio and Janhunen, 2002).However Phobos 2 observations (Dubinin et al., 1998) and, more recently, the MGS observations (Vignes et al., 2002) showed that the shock is observed at larger distances in the hemisphere in which the motional electric field is pointed out.A similar behaviour have been observed for the Venusian BS (Alexander et al., 1986).Authors suggest that mass loading of pick-up oxygen ions may produce this asymmetry, since planetary ions are accelerated mainly along the motional electric field.For the Martian bow shock, the observed asymmetry between the two hemispheres is weak: the average difference between the distance to the shock for the two hemispheres is about 0.2 R M (Vignes et al., 2002).Let remarks that these studies are based from a single spacecraft observations, and thus an assumption of a stationary IMF is needed to compare the influence of the IMF, and the motional electric field, with the bow shock distance to the planet.Though the simulations show a strong asymmetry of the oxygen ion fluxes (Modolo et al., 2005, see Figs. 5c and 5d), the BS is found to be farther in the opposite direction.Further investigations are required to rule out the effect of a numerical contribution which can explain the difference between the observations and the simulations: a best description of the Martian ionosphere combined to a whole and realistic 3-D neutral corona in the simulations is probably necessary.
Fine and complex structure of the BS, like multiple shock waves, are displayed in Fig. 1.The origin of these features is currently unknown.The width of one structure is about 600-900 km, which corresponds to 2-3 grid cells, and 4-6 proton Larmor radius.Effect of a finite grid size is not excluded and more investigations are required.Analogous structures are recognized in the solar wind proton and alpha particle density variations which are closely linked (Modolo et al., 2005).This splitting of the BS have been also identified in bi-fluid simulations (Sauer et al., 1996) and in hybrid simulations (Moore et al., 1991;Shimazu, 2001) but was not reliably experimentally observed in space, perhaps due to the inability to separate spatial structures and temporal variations of the BS with a single spacecraft.Further examinations are planned to study the persistence of such structures to different cell size and a different solar wind composition (without alpha particles).
A reasonable agreement is found between the location of BS in the simulations and observations (Table 1).Two parameters are used: R SD which is the distance along the X axis from the center of the planet to the maximum value reached at the shock and R T D which is the distance, along the Z axis (not including the asymmetry of the shock observed along the Y axis), from the center of the planet to the maximum value reached at the shock in the terminator plane.R SD and R T D refer to the nomenclature defined by Trotignon et al. (2006).The position of the different boundaries is estimated from the simulations at more or less one grid size (0.09 R M ).The interplanetary magnetic field (IMF) in the simulation is perpendicular to the plasma flow, and thus not follow the Parker spiral.This choice is motivated by the simplifaction of the simulation analysis and the understanding of the draping of the magnetic field.Although the chosen IMF is not the ideal average value, it is not unsual to find a similar case.Short term variations of IMF are important and the magnetic field can change from one direction to the opposite direction in few minutes: see for example solar wind data measured by the ACE spacecraft (and companions papers Stone et al., 1998) available through NOAA (http://sec.noaa.gov/ace/ACErtswhome.html).The BS and MPB observations may cover a large range of IMF directions, including the case of this study.The influence of the IMF on the plasma environment in the vicnity of Mars will be discussed in a future publication.
The simulated shock is found to be slightly farther from the planet than the shock fit from the MGS observations, nevertheless we estimate that the relative deviation on R SD and R T D parameters is less than 8%.
The solar wind interaction with Mars displays Venus-like features (Cloutier et al., 1999), but some difference have also been noticed.In the case of the solar wind interaction with Venus, the location of the BS is very sensitive to the upstream solar wind parameters, the orientation of the interplanetary magnetic field, and the solar cycle and the solar EUV flux (Russell et al., 1988).On the other hand, although the location of the Martian BS has been found to be highly variable, the solar activity and the solar EUV flux does not play a decisive role on the BS variations.By comparison of the Phobos 2 data, obtained at the maximum solar activity, and the first MGS results at the minimum solar activity, the mean location of the Martian BS appears to be independent of the solar cycle (Vignes et al., 2000).The simulation results occur in agreement with the Phobos 2 and MGS observations.The subsolar and the terminator positions of the BS are roughly the same for solar minimum and maximum conditions (Table 1).
Figure 2 shows a comparison of the total magnetic field and the magnetic field components at solar minimum and solar maximum along a circular Phobos-2 like orbit (r=2.8R M ) in the XZ plane.The simulated magnetic field component are presented in the Mars ecliptic coordinate system used by Luhmann et al. (1991) for a better comparison with the Phobos 2 observations.
The BS is easily identified by a sharp jump in the magnetic field value.The BS location corresponds approximately to the intersection between the circular orbit and the terminator plane.The draping of the magnetic field around the planet is clearly displayed in Fig. 2, with the presence of two magnetotail lobes with the opposite polarity.The two lobes are almost symmetric in this plane and a cross-tail current sheet is apparent.These results are in a good agreement with the Phobos 2 observations (Fig. 2 in Luhmann et al., 1991) and with simulations made by Kallio and Janhunen (2002).Another interesting point of this figure is the behaviour of the magnetic field components with a solar EUV change.The red and blue curves, corresponding to the simulation results for the maximum and minimum solar EUV flux, respectively, are very similar.The global position and the morphology of the BS, the magnetosheath and the tail lobes are weakly dependent on the solar EUV flux.Althoug a mass loading and the production rate of oxygen ions are greater for maximum solar EUV flux, the density variations of the planetary plasma between minimum and maximum solar conditions are not sufficient to have an impact on the position and the shape of the BS.

The MPB
The MPB is a thin, sharp and permanent plasma boundary located downstream of the BS and upstream of the ionopause Vignes et al. (2000).The magnetic field lines pile up and drape around the planet with a strong enhancement (Bertucci et al., 2003b(Bertucci et al., ,a, 2005a)).This boundary can be identified by a generally sharp increase of the magnetic field, a strong diminution of variability of the magnetic field magnitude and orientation with respect to the magnetosheath, and a sharp decrease of electron suprathermal electron flux (Bertucci Fig. 2. Comparison between the total magnetic field and the magnetic field components observed by the Phobos 2 spacecraft (top) (reproduced from Luhmann et al., 1991) and the simulation results (bottom) plotted in the Mars ecliptic coordinate system.The simulated magnetic field is shown at the minimum (curve in blue) and maximum (curve in red) solar EUV flux as a function of the angle measured from the -X direction along a circular orbit with a radius of 2.8 R M in the XZ simulation-plane (similar to the Phobos 2 orbit).et al., 2005a).In the simulations the MPB is also well identified from a jump in the magnetic field intensity.In addition, close to this boundary a change in the plasma composition occurs and planetary ions dominate downstream the MPB (Lundin et al., 1991;Dubinin et al., 1996;Modolo et al., Table 2. Comparison between the fits from the Phobos2 and MGS observations (Trotignon et al., 2006) and simulated MPB locations at minimum and maximum solar EUV fluxes.Fedorov et al., 2006).Figure 1 clearly displays the MPB and, similar to the BS features, a strong asymmetry is marked in the XY plane while the XZ plane is nearly symmetric.The maximum of the magnitude of the MPB is reached in the hemisphere in which the motional electric field is pointed out.In the opposite hemisphere the magnetic signature of the MPB is almost absent and the MPB does not seem to be characterized by a jump in the magnetic field intensity.Analogous asymmetry have been reported by (Mazelle et al., 2004, see Figs. 6 and 7) in the previous hybrid simulations performed by Brecht (2002).The MGS observations show a similar asymmetry (Vennerstrom et al., 2003).Meanwhile this model does not take into account the crustal magnetic field, the position and the shape of the MPB obtained from the simulations are in a reasonable agreement with the fits from the MGS observations (Trotignon et al., 2006, cf.Table 2).The crustal magnetic field will be included in a future study, with an improved spatial resolution of 150 km accompanied with a better description of the ionosphere.Note here that Crider et al. (2002) and Dubinin et al. (2006) 1 have reported a dependence of the MPB position on the presence of crustal fields.
Another magnetic signature of the MPB associated to a change in the magnetic field topology suggested by Bertucci et al. (2003aBertucci et al. ( ,b, 2005b) ) is not so obvious in the simulations.Although a draping enhancement is identified downstream the MPB, it appears that the difference between the diagnostic proposed by Bertucci et al. (2003a) and performed in the magnetosheath and downstream the MPB is not so significant.Numerical diffusion, smoothing and poor spatial resolution limit a variability of the magnetic field in the magnetosheath, both in the orientation and in the magnitude, and thus smooth the variations in the magnetic field topology.
Table 2 compares the simulated MPB positions, both in maximum and minimum solar activity, with fits derived from the MGS observations.Similar to Table 1, R SD and R T D are the standoff distances to the MPB in the subsolar region and at the terminator along the Z axis.First, a good agreement between the simulations and observations is found on the location of the boundary.Secondly, the position of the MPB is weakly dependent on the solar EUV flux, that is in a concordance with the Phobos 2 and MGS observations (Trotignon et al., 2006).
Figure 3 displays the magnetic field intensity and the density profile in the oxygen corona along the Sun-Mars axis for both minimum and maximum solar EUV flux. Figure 3 shows clearly that the MPB position is concomitant to a change in the oxygen density profile.The lower altitude part of the density profile describes the thermal oxygen population while the upper altitude part is associated to the suprathermal oxygen atoms.The MPB is located in the transition region.The thermal oxygen atmosphere is weakly dependent on the solar EUV flux and even an increase of the photoionization frequency by a factor 3 and the corresponding increase of the production rate of oxygen ions does not seem to be sufficient to push away the MPB significantly.

Conclusions
A 3-D hybrid model has been used to characterize the two main Martian boundaries: the bow shock and the magnetic pile-up boundary.The production rates are computed locally and self-consistently in this model using the local neutral densities and ionization frequencies or cross-sections of the ionization processes, thus the mass loading is considered self-consistently.We have found that the locations of the simulated bow shock and the magnetic pile-up are in a good agreement with the fits from the Mars Global Surveyor observations for both the position and the shape.
A strong asymmetry has been noticed on both the bow shock and the magnetic pile-up boundary with respect to the motional electric field direction.Multiple shock structures have also been shown.Meanwhile the planetary plasma is very sensitive to the solar EUV flux (Modolo et al., 2005), the bow shock and the magnetic pile-up boundary are weakly dependent on this parameter that is in a concordance with the observations (Vignes et al., 2000;Trotignon et al., 2006).It is believed that the position of the "obstacle boundary" at Mars which controls the BS location coincides with the MPB.The situation with the MPB is more complicated since the origin of such a sharp boundary remains unclear yet.
Several more questions need to be investigated, and, particularly, a different asymmetry of the bow shock revealed in hybrid simulations and observations.

Fig. 1 .
Fig. 1.Maps of the magnetic field strength at solar minimum condition (panels (a) and (b)) and at solar maximum condition (panels (c)and (d)).Panels (a) and (c) display the total magnetic field in the planes containing the solar wind velocity and the motional electric field, while panels (b) and (d) show the total magnetic field in the plane containing the the solar wind velocity and the interplanetary magnetic field.The solid and the dashed black curves represent fits to the observations made by Phobos 2 and MGS, for respectively the bow shock and the MPB(Trotignon et al., 2006).

Fig. 3 .
Fig. 3. Magnetic field magnitude (black curves) and atomic oxygen density (red curves) along the Sun-Mars axis for both minimum (dashed curves) and maximum (solid curves) solar EUV flux.Vertical dotted lines indicate the subsolar positions of the bow shock and the magnetic pile-up boundary derived from the MGS observations.

Table 1 .
Comparison of the BS position between the Phobos 2 and MGS observations and simulations for solar minimum and maximum conditions.R SD (R M ) R T D (R M )