Energetic protons at Mars : interpretation of SLED / Phobos-2 observations by a kinetic model

Mars has neither a significant global intrinsic magnetic field nor a dense atmosphere. Therefore, solar energetic particles (SEPs) from the Sun can penetrate close to the planet (under some circumstances reaching the surface). On 13 March 1989 the SLED instrument aboard the Phobos2 spacecraft recorded the presence of SEPs near Mars while traversing a circular orbit (at 2.8 RM). In the present study the response of the Martian plasma environment to SEP impingement on 13 March was simulated using a kinetic model. The electric and magnetic fields were derived using a 3D self-consistent hybrid model (HYB-Mars) where ions are modelled as particles while electrons form a massless charge neutralizing fluid. The case study shows that the model successfully reproduced several of the observed features of the in situ observations: (1) a flux enhancement near the inbound bow shock, (2) the formation of a magnetic shadow where the energetic particle flux was decreased relative to its solar wind values, (3) the energy dependency of the flux enhancement near the bow shock and (4) how the size of the magnetic shadow depends on the incident particle energy. Overall, it is demonstrated that the Martian magnetic field environment resulting from the Mars–solar wind interaction significantly modulated the Martian energetic particle environment.


Introduction
The active Sun is a source of energetic particles which can potentially impact planetary bodies in the Solar System, particularly in the inner heliosphere.Significant solar flares can be associated with the acceleration of particles up to several hundred MeV/nucleon (in some instances up to a few GeV/nucleon).Such solar energetic particle (SEP) events are mostly composed of protons with about 10 % He and <1 % heavier elements.
A simplified impulsive/gradual scheme of SEP events has been introduced to distinguish between two types of SEP events (Cane et al., 1986).In this scheme impulsive events are defined to be of relatively short duration (<1 day) and to display a high proton content.Gradual SEP events on the other hand are of longer duration (days), have higher fluxes, display a wider spread in longitude and are associated with fast coronal mass ejection-driven shocks.These enhancements typically display a rapid rise in proton fluxes on a time scale of tens of minutes to an hour.They attain a maximum that is sometimes followed by a second, occasionally higher, intensity peak at energies <50 MeV which is recorded when the interplanetary shock reaches the observer (Reames, 1999).Overall, the particle profile recorded depends on where the observer is located relative to the moving shock source.
SEPs pose a hazard to technical devices and to life (for example on manned space missions).The effects of SEPs are most pronounced near objects that do not have significant intrinsic magnetic fields, which could shield their near-space environments from energetic particles.Moreover, if an object has a tenuous atmosphere, an energetic particle can penetrate deep into that atmosphere and may even reach the surface of A schematic illustration of the attitude of the SLED telescopes, the orbits of Phobos-2 and the main plasma regions.The figure shows representative encounters of Phobos-2 with Mars as indicated by plotting a circular orbit (C).Also represented are the bow shock (BS), the magnetosheath, magnetic pileup boundary (MPB), the magnetotail and the direction of the interplanetary magnetic field (IMF).The data analysed in this paper come from one of the circular orbits traversed.The conical shape illustrates the field of view of the SLED instrument while traversing a circular orbit near the terminator plane.
the body.This is the case at Mars (see e.g.Leblanc et al., 2002), since this planet does not have a significant global intrinsic magnetic field and has an atmospheric pressure which is only about one per cent of that of the atmospheric pressure at the Earth.
Pioneering measurements of energetic particles (range ∼30 keV to a few tens of MeV) were recorded at Mars in 1989 by the SLED (Solar Low Energy Detector) instrument on the Phobos-2 mission (McKenna-Lawlor et al., 1992).More recently, the MARIE (Mars Radiation Environment Experiment) aboard NASA's Mars Odyssey spacecraft, which was launched in 2001 (Zeitlin et al., 2004;Luhmann et al., 2007), recorded energetic particles with kinetic energy >30 MeV over approximately one year.Shortly, energetic particle measurements will be obtained on the Martian surface by the Radiation Assessment Detector (RAD) on the Mars Science Laboratory (MSL) mission, which was launched in 2011 (Hassler et al., 2006) as well as around Mars by the solar energetic particle (SEP) experiment on the MAVEN (Mars Atmosphere and Volatile EvolutioN) mission, which is scheduled for launch in 2013.

The present study
The purpose of the present paper is to extend an earlier novel analysis of the SLED energetic particle measurements based on a global hybrid model of the Mars-solar wind interaction (McKenna-Lawlor et al., 2012).The model was found, us-ing representative upstream conditions, to reproduce the socalled "magnetic shadow" behind Mars when viewed along the upstream direction of the interplanetary magnetic field, IMF, and was qualitatively in agreement with the SLED observations (McKenna-Lawlor et al., 2012).Overall, the simulations suggested that the configuration of a magnetic shadow depends on the solar wind density and velocity, and on the magnitude and direction of the interplanetary magnetic field.
In the present study the response of the Martian plasma environment to SEP impingement on 13 March 1989 was simulated using a kinetic model.The electric and magnetic fields were derived using a 3-D self-consistent hybrid model (HYB-Mars) in which ions are modelled as particles while electrons form a massless charge neutralizing fluid.Plasma and magnetic field data measured contemporaneously with the particle measurements aboard Phobos-2 were used as inputs to the model.A case study of particle circumstances on 13 March was selected for investigation since the IMF remained stable over a relatively long time during the in situ measurements.
The paper is organized as follows.First, in Sect. 2 energetic particle observations recorded by the SLED instrument aboard Phobos-2 along a circular orbit at 2.8 R M are introduced.In Sect. 3 the basic properties of the model are described.In Sect. 4 the simulation results are presented.Then the spatial distribution of energetic ions near Mars is investigated.After that the flux of energetic protons is derived along the orbit of Phobos-2 and the simulated fluxes are compared with the observed fluxes.Finally, similarities and differences between the modelled and observed values are discussed.

Phobos-observations
The energetic particle observations analysed in this paper were made by the SLED instrument on-board the Phobos-2 mission.The SLED instrument featured two telescopes (named Te1 and Te2), both of which viewed in the same direction, namely 55 • westward of the Sun-Mars line -which is approximately in the nominal direction of the interplanetary magnetic field at Mars (Fig. 1).Both telescopes had a 40 • apex angle.The difference between Te1 and Te2 is that Te2 was equipped with an additional Al foil which absorbed protons with energies <350 keV while allowing the detection of 35-350 keV electrons.The count rate differences between Te1 and Te2 allow proton and electron fluxes to be differentiated (see the Appendix).The data analysed in this paper were obtained along circular orbits of Mars at r = 2.8 R M .During the circular orbits the Phobos-2 spacecraft was near the orbital plane of the Phobos moon which has an inclination to the ecliptic of 26 • .
The interaction geometry of Phobos-2 with the Mars environment is shown in Fig. 1 for 13 March 1989.The inbound bow shock crossing (BS IN ), when Phobos-2 entered from the solar wind into the magnetosheath while in circular orbit, and the outbound BS crossing, (BS OUT ), when Phobos-2 exited from the magnetosheath to the solar wind, are highly asymmetric from the viewpoint of the SLED measurements.In the IMF condition shown in Fig. 1, the BS IN crossing is in the region of a perpendicular BS and the BS OUT crossing is in the region of a parallel BS.Therefore, in the former case a clear BS crossing is anticipated to be seen, whereas in the latter case the BS crossing is anticipated to take place under turbulent conditions and to be not so well defined.Moreover, when SLED made measurements near BS IN , the fieldsof-view (FoVs) of the SLED telescopes were first generally along the surface of the BS, whereas later in the magnetosheath the FoVs were toward the surface of Mars.Near BS OUT the FoVs were toward the solar wind or, more precisely, toward the foreshock region.These geometries have implications with regard to the interpretation of the data, as will be discussed later in Sect. 4.
Phobos-2 performed 114 circular orbits around Mars in February-March 1989.During March the SLED instrument recorded several particle flux enhancements in association with major ongoing solar flares (see e.g.McKenna-Lawlor Table 1.Energy channels of two telescopes (Te1 and Te2) of the SLED instrument aboard the Phobos-2 spacecraft (McKenna-Lawlor et al., 1992).The abbreviations in parentheses are used later in this paper when referring to the energy channels.[−600, 0, 0] km s −1 The interplanetary magnetic field (B sw ) [−4, 4, 0] nT et al., 2012, Fig. 1).For the present paper one of these events was chosen for detailed analysis using the following criteria.First, and most importantly, there had to be magnetic field measurements available for the analysed SLED event.Second, the IMF had to be as stable as possible because the measurements would be compared with a simulation run made using constant upstream IMF conditions.The stability of the IMF was defined by determining the IMF just before BS IN and just after BS OUT .The time difference between BS IN and BS OUT was over 3 h and in all of the circular orbits the IMF varied to some extent between BS IN and BS OUT .Third, in order to ensure that SLED collected as large a sample of the energetic protons as possible, the direction of the IMF was required to be relatively close to the spiral angle of 55 • toward which the SLED telescopes pointed.An additional advantage was gained if the velocity and/or the density of the solar wind could be evaluated from contemporaneous Aspera/Phobos-2 observations.
The top of Fig. 2 shows the SLED particle data from Te1, and the bottom of Fig. 2 presents corresponding magnetic field data recorded aboard Phobos-2 while in circular orbit on 13 March 1989.In Fig. 2, electrons recorded by the instrument have been subtracted from the ion counts using the algorithm described in the Appendix.
A characteristic feature of the IMF in the analysed orbit is that it is quite stationary in the solar wind before BS IN but highly fluctuating after BS OUT .This implies that sometime between the inbound and outbound BS crossing the properties of the IMF changed substantially.In the present paper the IMF used in the simulation ([−4.0,4.0, 0.0] nT) was determined just before BS IN and, thus, it is anticipated that the accuracy of the modelled energetic proton fluxes decreased from BS IN to BS OUT along the orbit of Phobos-2.Therefore, in what follows, detailed comparisons between the modelled and observed particle fluxes are focused on the region near BS IN .
One can identify several plasma and magnetic field regions in Fig. 2. In the magnetic field a clear crossing of the (perpendicular) bow shock took place at about 00:20 UT.Before that the spiral angle in the range [0 • , 360 • ], represented in Fig. 2 by the angle between the IMF vector projected on the z = 0 plane and the −x-axis, where the angle increases clockwise when viewed on the +z-axis, is close to the FoV of the SLED telescopes (55 • ).Moreover, the IMF clock angle (the angle in the range [0 • , 360 • ] between the IMF vector projected on the x = 0 plane and the +z-axis (where the angle increases clockwise when viewed from the Sun)) is close to 90 • , implying that the IMF was near the orbital plane when the spacecraft was near BS IN .Crossing of the bow shock at ∼00:20 UT can be identified clearly in the SLED data by flux enhancements in the low energy channels.After the inbound BS crossing, Phobos-2 entered the so-called magnetic shadow region where a decrease in the particle flux occurred in all the SLED energy channels (see McKenna-Lawlor et al., 2012, for a detailed discussion of the magnetic shadow at Mars).

Description of the model
The kinetic model developed in this study is based on a selfconsistent global hybrid plasma model (HYB-Mars) of the Mars-solar wind interaction.The hybrid model gives the global electric field (E) and the global magnetic field (B) around Mars, which are used to trace energetic particles by the Lorentz force.The HYB-Mars plasma model is described in detail elsewhere (Kallio et al., 2010), and here only its basic properties are summarised.A new feature in the HYB-Mars model particular to this paper is the high energy particle tracing mode, in which the injection of energetic protons into the simulation is implemented after a stationary global solution for the Mars-solar wind interaction is reached.The density of the injected SEP protons is relatively low compared to that of the solar wind, and, thus, they do not affect the global structure of the Martian-induced magnetosphere.
Special emphasis is placed here on describing how energetic protons are injected into the simulation.An important feature in the HYB-Mars high energy particle tracing mode is that the incident direction of the energetic protons can be arbitrarily chosen.This feature has to exist in the model because energetic protons move predominantly along the direction of the IMF, which is typically not along the Mars-Sun line (as is the case for the main ∼1 keV solar wind proton population).The three-dimensional (3-D) velocity distribution of the energetic protons impacting on Mars is not known, and several trial velocity distribution functions were investigated here when interpreting the SLED observations.

Hybrid model HYB-Mars
In the hybrid simulation, ions are modelled as particles while electrons form a massless, charge neutralizing fluid.The approach is self-consistent, as is the case in magnetohydrodynamic (MHD) models.The advantage of the hybrid model for studying ions and kinetic effects is that the ion velocity distribution can be fully 3-D.In the HYB-Mars model ions are accelerated by the Lorentz force and the equation of motion of ions is (1) Here m i , q i , v i and r i are the mass, electric charge, velocity and position of an ion, E is the electric field and B is the magnetic field.The model automatically includes finite ion gyroradius effects and does not contain any pre-assumption regarding the velocity distribution function of the ions.It is recalled for comparison that, in an MHD model, the velocity distribution function is assumed to be Maxwellian.The 3-D hybrid model used in this paper, HYB-Mars, was used in several earlier studies of the Mars-solar wind interaction (for example, Kallio et al., 2010, and references therein).The model contains planetary O + and O + 2 ions.The O + ions originate from the spherically symmetric oxygen corona and from the ionosphere.The O + 2 ions come from the ionosphere.The model does not contain a self-consistent ionosphere, and so ionospheric O + and O + 2 ions are simply injected from a spherical shell located 210 km above the surface of Mars.Ionospheric O + and O + 2 ions are each assumed to have an initial temperature of 10 5 K, and the initial bulk velocity of these ions is taken to be zero.On the dayside the oxygen emission flux from the obstacle, Q [s −1 m −2 ], depends on the solar zenith angle (SZA) as Q ∼ cos(SZA).On the nightside Q is constant.The simple ionospheric model adopted implies that emission of ionospheric ions in the hybrid simulation is axially symmetric around the Mars-Sun line.The model also incorporates H + ions from the exosphere that originate from the spherically symmetric hydrogen corona.Martian magnetic anomalies are not included in the model.The lack of a self-consistent ionosphere and magnetic anomalies is not, however, anticipated to cause significant uncertainties in the results presented in this paper, because the properties of energetic protons are studied relatively far from the planet at r = 2.8 R M (R M = 3390 km is the value for the Mars radius used in the simulation).
The coordinate system of the model is as follows: The solar wind flows along the −x-direction; the z-axis points to the north and is perpendicular to the orbital plane of Mars.The y-axis completes the right-handed coordinate system.The dimensions of the simulation box in the present study were −4.0 A hierarchically refined grid is employed in the HYB-Mars model.The grid size is 678 km when r > 3.5 R M , and 339 km below 3.5 R M .The time step is 0.02 s from t = 0 s up to t = 300 s when the hybrid simulation has reached a stationary solution.From t = 300 s onwards, energetic proton injection is started, accompanied by a reduction in the time step to 0.005 s (three lowest energy populations) and 0.001 s (two highest energy populations) to account for higher ion velocities.These small time steps ensure that a high energy ion does not move over a whole grid cell during one time step (which would result in artificial density fluctuations in the simulation box).Injection was continued until enough energetic protons (2 × 10 6 per ion population) were collected, which was at around t = 315 s (for an account of the various ion populations involved, see Sect.3.2).
The upstream parameters used in the case study for 13 March 1989 are based on Phobos-2 observations.The bulk velocity (U sw ) was approximated to be [−600, 0, 0] km s −1 based on Aspera/Phobos-2 measurements.The interplanetary magnetic field (B sw ) was taken to be [−4, 4, 0] nT, based on MAGMA/Phobos-2 measurements, as discussed in Sect. 2. The density of the solar wind (n sw ) cannot be determined accurately from the ASPERA observations (see Kallio et al., 1994, for a discussion of the inaccuracies of the ASPERA density and bulk velocity measurements).The most appropriate value for the solar wind density was determined by making several runs at different density values (3 cm −3 , 4 cm −3 , 5 cm −3 and 6 cm −3 ) and comparing the simulated magnetic field along the orbit of Phobos-2 with the observed magnetic field.The simulated magnetic field showed similar behaviour at all four densities, and the changes in the magnetic field values between the runs were not significant.The smallest modelled density of the solar wind resulted in the most distant bow shock.Further, the position and magnetic features of the simulated inbound bow shock crossing (BS IN ) were in best agreement with the observed position and features of BS IN .Therefore, n sw = 3 cm −3 was adapted in the case study as the density of the solar wind.Table 2 summarises the upstream parameters used in this paper in the hybrid model.

HYB-Mars: high energy particle tracing mode
The injection of energetic protons into the HYB-Mars simulation was started at t = 300 s when the global solution had reached a stationary state.The simulation run at t > 300 s is referred to in this paper as the high energy particle tracing mode, because (1) the density of the energetic protons was kept low as compared with the solar wind, and thus the energetic protons did not affect the self-consistent simulation, and (2) because the time step was set to be sufficiently small that the trajectories of the energetic protons could be modelled accurately.
The energetic particles in the solar wind were accelerated by the Lorentz force (Eq. 1) in a similar way to the acceleration of all the ions in the HYB model.The incident energetic protons moved predominantly along the direction of the IMF, and they were injected into the simulation box on the x-face (x = 3.0 R M ), on the y-face (y = −4.0R M ) and on the z-faces (z = ±4 R M ).The time steps in the high energy particle tracing mode were reduced to 0.005 s and 0.001 s in order to ensure that the energetic ions did not jump over the smallest grid cells.All the energetic protons had a statistical weight factor, w i , which indicates how many real protons correspond to a simulated macroparticle (i).
One of the most critical decisions to make in constructing the model concerned the question as to what kind of incident upstream 3-D velocity distribution should be used www.ann-geophys.net/30/1595/2012/for the ∼30 keV < E <∼ MeV protons.SLED observations showed that the energy spectra, F [(cm 2 s ster keV) −1 ], can be modelled using a power law F = αE −K , where α is the amplitude and K is the spectral index (Afonin et al., 1991).In this paper the velocity distribution function characterized by a continuous energy spectrum is referred to as showing "energy scattering".The term "full scattering" is used to describe the distribution in a situation where protons do not flow exactly along the IMF but have both energy scattering and angular scattering (as will be described in detail below).The third tested velocity distribution is referred to as the "energy beam" distribution.In this case the velocity distribution is taken to be made up of several discrete ion beams which do not display either energy or angular scattering in velocity space.In all three methods the injection of the energetic proton population is spatially homogenous and temporally stationary.Figure 3 illustrates these velocity distributions.
The beam distribution contains six energy beams corresponding to the threshold energies of energy channels 1-5 of the SLED instrument (cf.Table 1): 30 keV, 50 keV, 100 keV, 200 keV, 600 keV, and 3.2 MeV.These energy values were chosen to probe the limiting behaviour of particles in different energy channels.The 100 keV beam was included to obtain a better energy resolution at the low energies where there were significant changes in beam behaviour with respect to different energies.The chosen energies provide a rough estimate of the energy spectrum dependence within a channel, when the upper and the lower limiting energy populations are combined into counts for the corresponding channel.The ions in the beams were injected with equal, negligibly small densities.During analysis, the beams were weighted according to a power law dependence (∼ E −1.5 ) as discussed later.The collected counts were then converted into simulated channel counts as a weighted average (by normalizing to upstream fluxes).
In the energy scattered distributions, the energy spectra [(cm 2 s ster keV) −1 ] were taken to follow ∼ E −1.5 , which roughly corresponds to the energy spectrum obtained from the SLED data measured in the solar wind before the bow shock crossing in the analysed orbit using a least squares fit.To compare with earlier research, McKenna-Lawlor et al. (1992) found similar SEP energy spectra upstream of the interaction region, with spectral indices around −1.75.As can be seen in Fig. 3b, all protons are on the same line (l SEP ) in velocity space as in the beam distribution.The orientation of the l SEP line in velocity space is given by the vector U SW + w || B SW /|B SW | where w || is (as will be discussed in detail later in Sect.3.2) the ion velocity in the moving solar wind frame of reference along the IMF.The continuous beam consists of five proton populations, corresponding to the (1-5) SLED energy channels.The individual populations follow a power law spectrum, with diminishingly small densities.The populations were afterwards combined and binned according to the individual SLED channels to account for acceleration/deceleration of particles in different energy regimes.
In the full scattered distribution model, protons no longer move along the l SEP line in velocity space but display angular scattering around l SEP .The distribution is implemented identically in the energy scattered distribution, notwithstanding the angular distribution.The simple pitch angle distribution function made for this paper comprises a Gaussian distribution, centred on the streaming axis and cut off at the 45 • pitch angle full width at half maximum 30 • .The choice of notable angular spread is motivated by the analysis of angular distributions of SEPs (Reames et al., 2001), particularly those associated with (small) gradual events which show a significant spread in the SEP magnetic azimuth distribution.Albeit the analysis was made using WIND spacecraft data measured in the vicinity of Earth, it turned out that the assumed full scattered distribution produced the most credible result using our model (as will be seen in Sect.4.2).Note that in Fig. 3c the same number of macroparticles is included in every SLED energy channel, and, therefore, the density of points in Fig. 3c varies from energy channel to energy channel.
To summarize, the high energy protons in our study were modelled using three different velocity distribution functions: 1. a "beam" distribution (Fig. 3a) composed of tight, discrete energy beams.
2. an "energy scattered" distribution (Fig. 3b), i.e. a unidirectional beam in velocity space with a continuous energy distribution.
3. a "full scattered" distribution (Fig. 3c), i.e. a population continuous in its energy and pitch angle distributions.

Trajectories of energetic particles
Energetic protons rotate, i.e. perform gyromotion, around the magnetic field lines of the IMF.The important feature in an energetic particle simulation is that the particle generator has to produce a particle population in the undisturbed solar wind whose properties remain unchanged from point to point.Otherwise the properties of the energetic protons in the region of space analyzed would depend on the positions where the particles were injected into the simulation.
In order to reach this goal, the energetic particle generator has to produce protons for which the velocity distribution function remains everywhere similar in the undisturbed solar wind, where E and B are constant.Without paying special attention to the energetic particle generation function, artificial variations in the velocity distribution can be introduced which in turn manifest themselves as non-constant plasma bulk parameters (density, bulk velocity, particle flux).In theory it is always possible for an energetic proton in a constant electric and magnetic field to move exactly along a straight line in space (i.e. to follow a trajectory along a line).This can be seen by noting that the Lorenz force, F L , is zero for a charged particle q i having a velocity, v i , equal to where B • SW ≡ B SW /|B SW | and w || is the ion velocity in the solar wind rest frame of reference along the IMF.The convective electric field E SW (≡ −U SW × B SW ) is zero in the solar wind rest frame.This means that one can obtain straight trajectories by launching ions along the IMF in the moving solar wind reference frame.The velocities of those protons remain unchanged in the undisturbed solar wind, and, consequently, they produce a particle population for which the density and bulk velocity are each constant in the solar wind.
The straight trajectories are, however, not aligned with each other for different w || values.Therefore, none of the straight trajectories are exactly along the direction of the IMF.The larger is w || , the more the direction of a line trajectory can be along the IMF.Moreover, the difference between a line trajectory and the direction of the IMF is largest at the lowest energies.Energetic protons in the beam distribution (Fig. 3a) and in the energy scattered distribution (Fig. 3b) are therefore straight lines in velocity space, and the directions of these lines are not exactly aligned with the direction of the IMF.
In the full scattered distribution (Fig. 3c), energetic protons execute gyromotion around the IMF because they have a non-zero velocity component, w ⊥ , perpendicular to the direction of the IMF in the rest frame of reference of the solar wind.In velocity space, an ion with a non-zero w ⊥ performs circular motion around the point w || B • SW in the solar wind rest frame and this ion is in the plane which is normal to the vector B • SW .This means that, if we follow the positions of the ions shown in Fig. 3c in time, they perform gyromotion around the B • SW vector.A time stationary velocity distribution can be obtained if the velocity distribution is symmetric around B • SW , as is the case for the fully scattered distribution shown in Fig. 3c.In a non-symmetric distribution the number of protons at different gyrophases would not be the same and the gyration of the ions in time would result in a "wobbling" velocity distribution accompanied by, for example, non-constant density and bulk velocity.

Magnetic field environment
The magnetic field is the key physical parameter required to understand the motion of energetic particles at Mars. Figure 4 illustrates the basic morphology of the near Mars magnetic field used in the analysed hybrid model run.Figure 4 shows magnetic field lines which are connected to the circular orbit of Phobos-2.The starting points for the field line tracings were chosen to illustrate that part of the orbit when the spacecraft moved from the solar wind to the deep magnetotail.
One can identify in Fig. 4   remained relatively unaltered after the crossing.Note that, as already mentioned, details of the morphology of the magnetic field near BS IN can be complicated because the field lines were very close the surface of the bow shock.Some of the magnetic field lines can, therefore, sometimes be in the solar wind and sometimes in the magnetosheath.In the hybrid model these field lines cannot be modelled very accurately because of the relatively large grid size (339 km) near the bow shock.
Point No. 3 is in the draped magnetic field line region which was reached when Phobos-2 moved from the magnetosheath into the magnetotail.At that point the magnetic field lines were still connected to the undisturbed IMF on the dayside.At point No. 4, Phobos-2 was deep in the magnetotail/magnetic tail lobe region where the magnetic field lines were highly draped and not connected to the upstream solar wind in the region shown in Fig. 3.
The simulated magnetic field derived along the Phobos-2 orbit is compared with the observed magnetic field in Fig. 5.The comparison shows the time range (00:00-02:00 UT) when the IMF can be assumed to have been most closely stationary during the analysed orbit.Figure 5 shows that the model can reproduce several general trends and events in the data: (1) a stationary IMF in the solar wind, (2) the crossing of the bow shock at ∼00:20 UT, (3) an increase in the mag-netic field after the BS crossing in the magnetosheath and (4) entry of the spacecraft into one of the magnetic tail lobes where B x is highly negative.
A detailed quantitative comparison shows that the modelled magnetic field does not, however, exactly reproduce the observed data.As discussed in more detail later in Sect.5, the differences found may either be associated with the limitations of the simulation or with the variations in the upstream parameters after Phobos-2 crossed the bow shock.In spite of the differences, similar trends in the model and in the observations shown in Fig. 5 suggest that that the HYB-Mars simulation can be used to provide a, relatively good, first-order approximation of the Martian magnetic field environment.

Energetic protons
It is useful to study in detail the motion of a set of individual energetic protons at Mars in order to interpret the simulated particle fluxes.Figure 6 shows the trajectory of several protons that were launched along the IMF in the solar wind rest frame toward Mars.The initial energies of the test particles (E = 30 keV, 50 keV, 200 keV, 600 keV, 3.2 MeV) were chosen to give information concerning ions in different energy channels of the SLED instrument (cf.Table 1).The electric and magnetic fields used in determining the Lorenz force were derived from the analyzed HYB-Mars simulation run.In Fig. 6 it is seen that the Martian magnetosphere causes more deviation in the directions of ∼30 keV and 50 keV ions (green lines) than in the directions of higher energy ∼600 keV and ∼3.2 MeV ions (red lines).Note also the large gyroradii of the energetic ions in comparison to the size of Mars and to the planetary solar wind interaction region.Further, in the simulation, some of the lower energy energetic ions (green lines) were reflected back into the solar wind when they encountered the Martian magnetosphere.It is noted that an earlier SEP simulation showed that protons with energies exceeding about 50 keV can hit the Martian exobase practically undisturbed (Leblanc et al., 2002).On the other hand, strong attenuation of the particle flux due to the influence of the Martian magnetosphere occurs at lower energies.
Let us now study the dynamics of the energetic ions near Mars.The magnetic field controls the motion of the particles seen in the highest energy channels of SLED (cf.Table 1) in channel 4 (0.6-3.2 MeV), channel 5 (3.2-4.5 MeV) and channel 6 (>30 MeV).The magnetic field changes only the direction of the ion velocities, but not their magnitudes.Only the electric field can change the kinetic energy of a charged particle moving in an electromagnetic field.The relative change of energy of the most energetic ions studied here was quite small compared to their incident energies, because the electric field cannot accelerate ions very much within the simulation box.An order of magnitude estimate for the maximum change of energy of an ion within the box can be obtained by considering an ion moving through the whole box along the z-axis.Then the energy change of the ion would be ∼65 keV (=600 km s −1 × 4 nT ×8 R M ).In practice, ions do not cross such a large potential drop so that the actual change in energy is less than this.The electric field can, however, result in notable relative energization of the protons in the low energy channels 1 (30-50 keV) and 2 (50-200 keV), particularly in the case of channel 1.
Figure 7 shows how energetic protons were accelerated and scattered near Mars along the orbit of Phobos-2 in the analysed simulation run.The horizontal axis represents the azimuthal, or longitudinal, angle of the Phobos-2 spacecraft along its orbit.The vertical axis defines the angle β between the ion velocity vector and the nominal SLED look direction d SLED (i.e. on the xy-plane with a spiral angle of 55 • ).That is, an ion moves exactly along the direction of the d SLED vector if the β angle is zero; it moves perpendicular to the d SLED vector if the β angle is 90 • and in the opposite direction to d SLED if the β angle is 180 • .Thus, the d SLED direction is a fixed direction within the model.The initial velocity distribution of the energetic protons was taken to be a beam (cf.Fig. 3a), and the initial velocities of the six proton beams were set at 30 keV, 50 keV, 100 keV, 200 keV, 600 keV and 3.2 MeV in order to study the motion of ions in the different energy channels of the SLED instrument (cf.Table 1).Note that the scattering was largest in the lowest energy channels.Moreover, a clear increase in energy can be seen in the low energy channels but not in the high energy channels.The angular and energy scattering of protons at BS IN at longitude ∼100 • and at BS OUT at longitude about 240 • are different because of the different geometries and, consequently, different magnetic morphology, at these bow shock crossings (cf.Fig. 4).It should be recalled that the ion trajectories are straight lines in the solar wind in the beam velocity distribution model.Therefore, as discussed in Sect.3.2.1 ("Trajectories of energetic particles") the ions do not flow exactly along the IMF even in the undisturbed solar wind; the β angle is not zero there, especially, in the case of the lowest energy channels.Note also that there are particles in the solar wind at high β angles (>100 • ), which are protons scattered upwind from the Martian magnetic environment, as seen earlier in Fig. 6.
Figure 7 includes also important information relevant to the interpretation of SLED data.The SLED telescopes have a 40 • apex angle, which means that, if SLED measured the ions in the situation shown in Fig. 7, it would only have detected ions at angles β ≤ 20 • .This implies that even if SLED had recorded all the ambient energetic protons when the spacecraft was in the solar wind, it could not have detected all these protons after the bow shock crossing.Note also that a significant portion of a SEP population might not hit the SLED aperture even in the undisturbed solar wind, but could yet be deflected into the aperture by the solar wind-planet Note that the red and yellow trajectories are almost straight lines and that they end when the ions hit the side wall of the simulation box.Note also that the colour of the initially green trajectories becomes grey or white when the trajectory is within or behind the strong magnetic field region.
interaction region; e.g. 30 keV and 50 keV beams closely straddle the β = 20 • line in the upwind region.A percentage of protons out of the FoV of the SLED telescopes depended on the energy of the particles, as can be seen by comparing the six panels in Fig. 7.
Before a detailed comparison between the simulated and observed particle fluxes along the analysed Phobos-2 orbit is made, it is useful to give an overview of the modelled particle fluxes at the orbital distance where the spacecraft made its measurements.From the data analysis point of view, a useful feature of the Phobos-2 circular orbit is that one can study the particle flux by deriving its values on a spherical shell having the radius of the circular orbit of the spacecraft.In this way one can estimate how the particle flux along the orbit of the spacecraft is sensitive to the direction of the IMF clock angle.Note that if the clock angle changes during a circular orbit, as is the case in reality, the particle flux on the spherical shell correspondingly rotates in the simulation around the x-axis because the properties of Mars and its planetary ions are axially symmetric in the hybrid model.The different IMF spiral angle cases cannot, however, be modelled by a simulation run made for fixed IMF conditions.Overall, such an overview helps visualization of the size of the regions where the flux is higher or lower than its undisturbed value, as well as indication of the sensitivity of the results to the IMF clock angle.
Figure 8 shows the simulated particle flux of the energetic protons in three SLED energy channels and the total magnetic field at r = 2.8 R M , i.e. the orbital radius of Phobos-2 at the time period considered.Only protons with β < 90 • are considered, i.e. protons moving towards the SLED instrument -back-scattered protons are ignored.The top panels provide a 3-D presentation where the particle fluxes are shown at the surface of the spherical shell.On the bottom panels the same data are shown in a 2-D latitude-longitude map.In the test particle simulation protons were recorded when passing the shell until the hit count reached two million hits per ion population.Afterwards the collected particles were binned according to the energies of the SLED channels.Spatial binning was achieved by dividing the shell using a rectangular 80 × 80 grid in latitude and longitude; that is, the size of the latitude × longitude grid was 2.25 • × 5 • .The accumulated particle flux was finally cosine-corrected for the local cross-section of the detector surface with respect to the direction of the particles.
One can identify a low particle flux region in all three energy channels behind Mars in Fig. 8.The low particle flux region is usually referred to as the magnetic shadow (see e.g.McKenna-Lawlor et al., 2012).Note that the size of the magnetic shadow decreases with increasing particle energy.Moreover, enhanced particle flux can be observed near the edges of the magnetic shadow.Previous simulations showed that this feature is a consequence of the fact that the directions of the high energy particles are not distorted as much as the directions of the low energy ions and that this is a consequence of finite ion gyroradius effects (McKenna-Lawlor et al., 2012).
Comparisons between the simulated particle fluxes and the SLED telescope 1 observations, from which electrons have been removed (see the Appendix), are displayed in Fig. 9.The measured particle fluxes are compared with several test particle simulations using three different velocity distribution models and two different detector aperture models.The velocity distribution models (the beam, the energy scattered, and the full scattered models) are provided earlier in Fig. 3.The two detector aperture models were the FoV of SLED and a FoV of 2π.All simulated energy fluxes were normalised independently so that the simulated flux in the solar wind was equal to the observed particle flux.
Figure 9 shows that the full scattered model provides better agreement with the measured data than the beam and energy scattered models -which result in an unrealistically high flux increase near the bow shock in channels 2-3.Moreover, those two models result in a flux increase near the bow shock in channels 4-5, although this flux is actually seen to decrease in the measured SLED data.In the two fully scattered simulation cases, the fluxes collected using the FoV of the SLED telescopes looking toward the nominal spiral angle of 55 • resulted in a slightly more realistic behaviour of the fluxes than when the particle flux was collected from 2π space.Note also that there is an energy dependence in the solar wind flux measured by SLED and reproduced by the model which occurs prior to the bow shock.This can be seen in Fig. 9 between 00:00 UT and 00:20 UT where one can see a gradual rise in the flux in channel 1 as the bow shock is approached.The amount of increase reduces in channels 2 and 3 while becoming flat in channel 4. In channel 5, in contrast, the flux decreased on approaching the bow shock.Overall, the angular scattered case resulted in a "smoother", more diffuse solution, which reproduced the SLED observations better than an initial distribution lacking angular scattering.

Discussion
In this paper a kinetic model has been presented and used to interpret the energetic proton environment at Mars under the extreme (SEP) conditions monitored in situ by the SLED instrument aboard Phobos-2 on 13 March 1989.
It turned out that even a fully scattered simulation case did not reproduce all the observed features near BS IN .In particular a decrease in the simulated flux in channels 3-5 started later in the simulated than in the observed case.This may be an indication that the fully scattered model is still an oversim-plified description of the true velocity distribution function in the analysed time frame.
Moreover, in the model it was assumed that the velocity distribution is time independent.The magnetic connection between Phobos-2 and the IMF near BS IN is very sensitive to the shape and position of the bow shock (cf.Fig. 4).Energetic particles follow more closely the magnetic field lines than do the low energy ions.The ion flux near BS IN is therefore anticipated to be very sensitive to the 3-D morphology of the magnetic field and its draping at the point where Phobos-2 enters the magnetosheath in the perpendicular bow shock region.Small changes in the position and shape of the bow shock are, therefore, anticipated to exert the strongest influence on high energy ions.
The properties of the bow shock in the model may include several inaccuracies associated with numerical issues (e.g. the finite grid size) or the upstream parameters utilized.The method used to filter electrons from the SLED data is also relatively simple.The attitudes of the FoV of the SLED telescopes are not exactly known because of the inaccuracy associated with the attitude of the Phobos-2 spacecraft (see e.g. as compared with the spatial size of the flux change region near BS IN .One of the biggest sources of inaccuracy is related to the 3-D velocity distribution of the energetic protons which is not known.The analysis indicated, however, that the best agreement between the simulated and measured fluxes was observed by using the most physically justifiable velocity distribution of the three models tested, which incorporated the effects of both energy and angular scattering in velocity space.
Overall, the study suggests that the model developed can be used to study the effects of energetic proton events at Mars.The data set used in this paper comes from 1989, but in principle the same tools could be used to model more recent observations at Mars by Mars Odyssey or the forthcoming observations by the MSL and MAVEN missions.The new high energy particle tracing mode implemented in the HYB model platform can in addition be used to study other plane-tary plasma interactions, for example, the solar wind interaction under extreme conditions at Mercury and Venus.

Summary
The global hybrid Mars (HYB-Mars) model was developed to study energetic protons near Mars.The model makes it possible to derive the electric and magnetic field for given upstream parameters and to include energetic particles in the simulation.The hybrid model has been upgraded by incorporating a so-called high energy particle tracing mode in the global simulation where energetic ions can be injected into the simulation using various manually inserted velocity distribution functions.Also, several diagnostic tools have been developed and incorporated into the HYB model platform in order to support one-to-one comparisons between the simulated particle fluxes and the observed particle fluxes.In particular, a virtual SLED instrument has been included in The upgraded model can reproduce several of the main characteristic features measured in situ at Mars under conditions when the simulation was considered to be most suited for making comparisons with measured data.The interval selected for these comparisons was a time period near the inbound bow shock crossing of Mars by Phobos-2 on 13 March 1989 when the upstream parameters were relatively stationary (so that the Mars-solar wind interaction could be modelled using constant upstream parameters).
Comparison with the SLED/Phobos-2 observations showed that the model successfully reproduced several of the observed key features: (a) a particle flux enhancement near the bow shock recorded in the low SLED energy channels, (b) a particle flux decrease near the bow shock in the SLED high energy channels, (c) the formation of a magnetic shadow, (d) the dependence on particle energy of the bow shock increase recorded at low energies and the size of the magnetic shadow.

Cleaning of electrons from the SLED data channels 1 and 2
In order to make one-to-one comparisons between simulated proton flux and SLED data, it is necessary to take into account the fact that SLED measured both ions and electrons (cf.Table 1).Comparison of counts in Te1 and Te2 gives a possibility to estimate the contribution of electrons to the SLED proton counts.
The accurate removal of electron counts from the SLED data is a complicated task, and the following, relatively simple and straightforward, method is used in this paper.An electron filtering algorithm was used to remove electrons from telescope 1 counts, resulting in a significant reduction in the overall counts and a prominent shadow in the data of telescope 1 channel 1 (which measured at the lowest energies recorded by the instrument).The algorithm uses ion data from the telescope 1 channels to remove ions from the channels of telescope 2. The remaining (electron) counts in the telescope 2 channels are then used to remove electrons from the channels of telescope 1 thereby creating a kind of ladder that reaches from the highest-energy channels to the lowest.The effect of different energy windows in the different channels and of the energy spectrum is estimated using a power law dependence on energy, assuming that both ions and electrons share the same, constant-in-time spectrum.
The algorithm uses the following estimator for obtaining reduced channel b counts c b reduced from channel a (possibly reduced) counts c a and channel b counts c b : Here A is the normalization constant (both ions and electrons have their own constants) and the indices a and b correspond to the number of the energy channel in a given telescope.Equation ( A1) is used to determine the value of A, which is then used in Eq. (A2).The energy windows from which counts c a and c b are collected are [E a1 , E a2 ] and [E b1 , E b2 ], respectively.The estimator takes into account the different sizes of the energy windows and also allows some extrapolation of the energy windows, as is the case in the first step.Negative reduced counts are set to zero.The value of the spectral index K in Eq. ( 3) was assumed to be a constant −1.5 during the whole measurement period for both ions and electrons (see also the discussion in Sect.3.2).All ions are assumed to be solar wind protons.Equation ( 3) is used recursively in order to remove electrons from energy channels 1 and 2 in telescope 1 in four steps (cf.Table 1 for the energy channel notation adopted): Step An example of the difference between the original SLED count rates and the reduced SLED count rates can be seen in Fig. 2. Note how the magnetic shadow becomes clearly visible in the lowest energy channel and slightly more visible in the second lowest channel.However the filtering of electrons does not significantly affect the higher energy channels 3-6, and so Fig. 2 shows only the original count rates for these channels.
Fig.1.A schematic illustration of the attitude of the SLED telescopes, the orbits of Phobos-2 and the main plasma regions.The figure shows representative encounters of Phobos-2 with Mars as indicated by plotting a circular orbit (C).Also represented are the bow shock (BS), the magnetosheath, magnetic pileup boundary (MPB), the magnetotail and the direction of the interplanetary magnetic field (IMF).The data analysed in this paper come from one of the circular orbits traversed.The conical shape illustrates the field of view of the SLED instrument while traversing a circular orbit near the terminator plane.

Fig. 2 .
Fig. 2. Measurements made along the circular orbit of Phobos-2 on 13 March 1989.Top: SLED/Phobos-2 count rates in Te1 in six energy channels.The energy ranges and the measured particles of Te1 are the following (cf.Table 1): channel 1 (34-51 keV, ions and electrons), channel 2 (51-202 keV, ions and electrons), channel 3 (202-609 keV, ions and electrons), channel 4 (0.6-3.2 MeV, ions), channel 5 (3.2-4.5 MeV ions), channel 6 (>30 MeV ions).The dashed and dotted lines show the raw data in channels 1 and 2. The solid lines show complementary data when the electron background has been removed from channels 1 and 2 (see the Appendix).No significant background was present in the higher energy channels.Middle: the corresponding magnetic field data recorded by MAGMA/ Phobos-2 are given in MSO coordinates.Bottom: the spiral angle, labelled as "Parker Angle", and the clock angle.The two vertical lines show the crossings of the bow shock: Phobos-2 crossed the bow shock (BS IN ) and entered the magnetosheath at around 00:30 UT and came back to the solar wind (BS OUT ) at around 03:40 UT.Note that the upstream conditions changed at around 01:40 UT.

Fig. 3 .
Fig. 3. Illustration of the incident upstream velocity distributions for the energetic protons tested in this paper.The dots correspond to macroparticles in the upstream SEP flow, with their colour corresponding to macroparticle weight.The energetic proton flux is modelled using three different initial velocity distribution functions for the high energy SW protons: (a) the beam distribution model with six energy beams corresponding to the SLED energy channels, (b) the energy scattered distribution model based on SLED data and (c) the fully scattered distribution model which includes both angular and energy scattering.The angular scattering is not known and in the runs it was assumed to be as shown in panel (c).Note that in panel (c) the distribution of macroparticles within an ion population goes as ∼ E −1.5 (applicable to panel (b) as well) and according to a Gaussian distribution in pitch angle, and therefore this is not a uniform distribution.

Fig. 4 .
Fig. 4. The simulated magnetic field around Mars and several magnetic field lines which are connected to the orbit of Phobos-2.The "cloud" in the figure is a volumetric plot of the magnitude of the magnetic field.The magnetic field lines are coloured in order to help the eye to distinguish different lines from each other.The view point is on the +z-axis.Note that the draped magnetic field lines are three-dimensional in the simulation and the figure only shows their projection on the xy-plane.The unit on the coordinate axis is 1000 km.See the text for a description of how the morphology of the magnetic field changes when Phobos-2 moves from the solar wind (point 1) deep into the magnetotail (point 4).

Fig. 5 .
Fig. 5. Comparison of the magnetic field modelled by the HYB-Mars simulation (dashed/dotted lines) and the magnetic field (solid lines) measured by MAGMA/Phobos-2 on 13 March 1989.Phobos-2 was initially in the solar wind, then crossed the bow shock (BS IN ) and entered the magnetic tail lobe.

Fig. 6 .
Fig. 6.Illustration of the motion of energetic protons near Mars.The "cloud" shows the strength of the magnetic field in 3-D (see the colour bar, top left).The circular white sphere represents the surface of Mars.The lines show trajectories of several energetic protons.The colour of a line represents the initial energy of an ion: The lowest energies (∼30 keV and 50 keV) are shown in green; the yellow trajectories represent ∼200 keV and 600 keV protons and the red trajectories ∼3.2 MeV protons (see the colour bar, bottom left).Note that the red and yellow trajectories are almost straight lines and that they end when the ions hit the side wall of the simulation box.Note also that the colour of the initially green trajectories becomes grey or white when the trajectory is within or behind the strong magnetic field region.

Fig. 7 .
Fig. 7. Modelled scattering of high energy solar wind ions in the Martian magnetosphere.The six panels show the energies of the studied energetic protons whose initial energies were (from top to bottom) 30 keV, 50 keV, 100 keV, 200 keV, 600 keV and 3.2 MeV.A dot corresponds to a proton collected near the Phobos-2 circular orbit (within a 100 km radius) on 13 March 1989.The vertical axis shows the angle β which is the angle between the ion velocity vector and the nominal SLED look direction.The colour of the dot gives the energy of the ion.Note that all the panels have their own colour scale which is shown on the right side of the panels.The horizontal axis shows the longitude of the position of the Phobos-2 spacecraft, and the angle is measured from the positive x-axis on the Martian orbital plane, with positive values toward the positive y-axis.Therefore, Phobos-2 was at the subsolar point at longitude 0 • and in the middle of the magnetotail at longitude 180 • .BS IN took place at longitude ∼100 • and BS OUT at longitude ∼240 • .

Fig. 9 .
Fig. 9. Comparison of the simulated energetic proton fluxes and the particle fluxes measured by the SLED/Phobos-2 instrument along the circular orbit of Phobos-2 on 13 May 1989.The SLED ion data come from Telescope 1, from which the electrons were removed using the method described in the Appendix.The five panels display the normalized particle fluxes in five SLED energy channels.The plots show derived fluxes based on three different velocity distribution models: the beam, the energy scattered, and the fully scattered models (see Fig. 3 for details).Two "fully scattered" fluxes are shown in cases where the particle fluxes were collected from 2π space (green solid line) and when the fluxes were collected from the field-of-view (FoV) of the SLED telescopes, looking toward the nominal spiral angle of 55 • at Mars.The fully scattered model results in the best agreement with the data.For comparison, see Fig. 5 for contemporaneously measured and for simulated magnetic fields.

cc
1: Use ion counts in Te1ch4 to remove ion counts from Te2ch3 to derive counts associated with electrons c Te2Ch3 (e − ) = c Te2Ch3 − c Use the electron counts obtained from Te2ch3 to remove electrons from Te1ch3 and obtain get pure ion counts from Te1ch3 c Te1Ch3 (H + )=c Te1Ch3 −c Te2ch3 e − the ion counts obtained from Te1ch3 to remove ion counts from Te2ch1 and Te2ch2 c Te2Ch2 (e − ) = c Te2Ch2 − c Te1Ch3 (H + ) Te2Ch1 (e − ) = c Te2Ch1 − c Te1Ch3 (H + ) 400 keV 350 keV E k dE 609 keV 202 keV E k dE (A2d) Step 4: Use the electron counts obtained from Te2ch1 and Te2ch2 to remove electron counts from Te1ch1 and Te1ch2 c Te1Ch1 (H + )= c Te1Ch1 −c Te2Ch1 (e − ) Te1Ch2 (H + )=c Te1Ch2 −c Te2Ch2 (e − )

Table 2 .
Solar wind parameters used in this paper in the hybrid model.
four magnetic field line regions which are individually labelled by the numbers 1-4.Point No. 1 is in the undisturbed solar wind where the magnetic field lines are in the undisturbed IMF.The Martian bow shock was crossed near point No. 2. This is the place of BS IN .Because this is a crossing of an almost fully perpendicular bow shock region, the direction of the magnetic field www.ann-geophys.net/30/1595/2012/Ann.Geophys., 30, 1595-1609, 2012