Modelling the solar wind interaction with Mercury by a quasi-neutral hybrid model

Quasi-neutral hybrid model is a self-consistent modelling approach that includes positively charged particles and an electron fluid. The approach has received an increasing interest in space plasma physics research because it makes it possible to study several plasma physical processes that are difficult or impossible to model by self-consistent fluid models, such as the effects associated with the ions’ finite gyroradius, the velocity difference between different ion species, or the non-Maxwellian velocity distribution function. By now quasi-neutral hybrid models have been used to study the solar wind interaction with the non-magnetised Solar System bodies of Mars, Venus, Titan and comets. Localized, two-dimensional hybrid model runs have also been made to study terrestrial dayside magnetosheath. However, the Hermean plasma environment has not yet been analysed by a global quasi-neutral hybrid model. In this paper we present a new quasi-neutral hybrid model developed to study various processes associated with the Mercury-solar wind interaction. Emphasis is placed on addressing advantages and disadvantages of the approach to study different plasma physical processes near the planet. The basic assumptions of the approach and the algorithms used in the new model are thoroughly presented. Finally, some of the first three-dimensional hybrid model runs made for Mercury are presented. The resulting macroscopic plasma parameters and the morphology of the magnetic field demonstrate the applicability of the new approach to study the Mercury-solar wind interaction globally. In addition, the real advantage of the kinetic hybrid model approach is to study the property of individual ions, and the study clearly demonstrates the large potential of the approach to address these more detailed issues by a quasi-neutral hybrid model in the future.


Introduction
The Hermean magnetosphere is anticipated to have certain unique features compared with other planetary magnetospheres.The atmosphere of Mercury is very tenuous, and the inner edge of the exosphere is the surface of the planet.The planet provides an example of a so-called surface-bounded atmosphere (Killen and Ip, 1999).In these circumstances the sputtered surface particles can escape from the surface practically without collisions with other atmospheric particles.By now, six neutral species have been positively identified: hydrogen (H), helium (He) and oxygen (O) by Mariner 10 and afterwards sodium (Na), potassium (K) and calcium (Ca) by Earth-based observations (see, Killen andIp, 1999, andref. therein, andBida et al., 2000).The relative importance of various processes responsible for emitting and accelerating particles from the surface and the role of the ions in the global Mercury-solar wind interaction are not well known.
The planet's intrinsic magnetic field is strong enough to form a magnetosphere around the planet.In magnetospheric physics, Mercury provides an example of a "miniature" magnetosphere, whose linear dimensions are only about 5% of those of the Earth (Russell et al., 1988).How the smaller spatial scale and the missing ionosphere affect Mercury's magnetospheric processes compared to the magnetospheric processes at Earth is unclear.How effectively can the magnetosphere shield the planet against the solar wind particles, how much can the magnetosphere store energy in its magnetotail, and how the interaction between the magnetosphere and the planet's surface takes place (which can take place because of the very tenuous atmosphere) are still open questions.The formation of the intrinsic magnetic field and its strength have also been under discussion ever since the Mariner 10 flyby in 1974.
Modelling of the Hermean plasma environment is a challenging problem because a physically realistic model should include processes between very different regions: the solar wind, the magnetosphere and the planetary surface.The modelling approaches that have been used so far to study near Hermean space can be divided into two categories: analytical models and fluid models.Analytical models are obtained by scaling down models that describe plasmas and/or fields at the terrestrial magnetosphere.As an example, the motion of ions in the Hermean magnetosphere has been studied by Tsyganenko's (Lukyanov et al., 1997), and Luhmann-Friesen's magnetic field model (Delcourt et al., 2002), and by an Earth's magnetotail model (Ip, 1997).Tsyganenko's (Luhmann et al., 1998) and Toffoletto-Hill's (Killen et al., 2001) magnetic field model have also been used to analyse how the Hermean magnetosphere is connected to the interplanetary magnetic field.
Although analytical models can give valuable new information about the physical processes near Mercury, they do not treat plasma and fields self-consistently.Spreiter and Stahara's gas dynamic (GD) model (see, for example, Spreiter andStahara, 1980, 1992) is the most well-known fluid model that has been used to analyse how the solar wind flows around a planetary body.The model solves gas dynamic parameters (particle density, bulk velocity and thermal pressure) in the planetary magneto-ionosheath.However, the model does not include magnetic field effects.The magnetic field, occasionally referred to as the convective magnetic field, is derived afterwards from the gas dynamic parameters, by assuming that the magnetic field is frozen into the gas dynamic flow.Recently, modern computational facilities have made it possible to model the Mercury-solar wind interaction self-consistently by a global three-dimensional (3-D) magnetohydrodynamic (MHD) model (Gombosi, 2000;Kabin et al., 2000;Ip and Kopp, 2002).
A quasi-neutral hybrid model may be regarded to provide a new third global modelling approach.As an MHD model, it describes the solar wind-planetary body interaction selfconsistently.A multi-ion hybrid model can include several ion species, such as the solar wind protons and planetary ions, say, Na + or K + ions, which are moved by the Lorentz force.Electrons form a massless charge neutralizing fluid.The advantage of the approach is that it makes it possible to study several plasma physical processes in the vicinity of the planet that are difficult or impossible to model by selfconsistent MHD models, such as the effects associated with the ions' finite gyroradius, the velocity difference between different ion species, the large ion loss cone, or the non-Maxwellian velocity distribution function.In particular, particles in the hybrid model can hit the surface of the planet, which makes it possible to study the effects and the role of the particle sputtering process, which have received notable interests (see, for example, Ip, 1997, andKillen et al., 2001, andreferences therein).So far 1-D quasi-neutral hybrid models have been used to study comets (see, for example, Omidi and Winske, 1986;Puhl-Quinn and Cravens, 1995), a 2-D model to study the terrestrial magnetosheath (Lin, 2001) and Venus (Terada et al., 2002) and 3-D global hybrid models to study Martian (see, for example, Kallio andJanhunen, 2001, 2002), Venusian (Brecht and Ferrante, 1991;Shimazu, 1999) and Titanian (Brecht et al., 2000) plasma environments.However, to the authors knowledge, a quasi-neutral hybrid model has not yet been applied to Mercury.
The paper is organised as follows.The basic assumptions of a hybrid model approach and the algorithms used in the developed new hybrid model for Mercury are first thoroughly presented.An emphasis is placed on addressing the advantages and disadvantages of the approach to study various Hermean plasma physical processes.The applicability of the algorithms used and the numerical solution is demonstrated by presenting an example of the magnetic field and plasma parameters in a "closed magnetosphere" case.Finally, the development of the hybrid model approach during the forthcoming years is addressed.

Description of an Adaptive Quasi-neutral Hybrid Model
The state of the quasi-neutral hybrid model consists of the position, r, velocity, v and weight, w, of each ion and the gridded magnetic field.

Ions and electrons
The interaction between the solar wind and Mercury is modelled by a new, adaptive quasi-neutral hybrid simulation code.The model can include solar wind protons (H + ) and several planetary ions, such as Na + and K + ions.Electrons are modelled as a massless fluid.The ions are moved by the Lorentz force: where m i , x i , q i and v i , are the mass, the position, the electric charge and the velocity of an ion i (i = H + , or a planetary ion).E, B and U e are the electric field, the magnetic field and the bulk velocity of electrons, respectively.In the model the ions are modelled as cubical "clouds" with homogenous charge density and the same size as the grid cell of volume V.A macroparticle gives a contribution to all cells that it overlaps.In such a volume-weighting scheme the charge assignment is performed as follows: if is the volume of a macroparticle j of type i (a proton or a planetary ion) of weight w j i which overlaps a grid cell at the point r = r k of a volume V , then the grid cell contains (3) ions due to the macroparticle j (see Hockney and Eastwood, 1988, p. 142-146, for details on the plasma cloud approach and the "cloud-in-cell" used, CIC charge assignment procedure).The ion density, n i , and the ion bulk velocity, U i , in the analysed grid cell is, therefore, obtained from the equations where the sum j is over all macroparticles of type i (i = H + or a planetary ion) which overlap the analysed grid cell with a volume V , and v j i is the velocity of the macroparticle j of type i.
Electrons are a massless fluid and their momentum equation is written as where σ is the electric conductivity.In our model, σ can be any pre-described function which depends on the position r.
If ideal conductivity is used (σ = ∞), Eq. ( 6) simplifies to an equation that states that the magnetic field is "frozen into" the electron fluid: Equations ( 6) or ( 7), are used to eliminate the explicit electric field term from the model equations, for example, from the Lorentz force, as is done on the right-hand side of Eq. ( 2).It should be noted that in Eqs. ( 6) and ( 7), zero electron temperature and, consequently, zero electron pressure is assumed for simplicity.The non-zero electron pressure term could also be implemented, but then one would have to specify how the electron pressure is derived from the known parameters, for example, by assuming that the electron fluid is adiabatic or that the electron temperature is constant.The effect of various electron temperature models could then be examined, but such an undertaking is beyond the scope of the present study.
The density n e and the (bulk) velocity U e of the electron fluid are obtained from the equations U e = (en e ) −1 (−j In Eq. ( 8) the sum i is over all ion species which implies quasi-neutrality: the electrons are connected electrically to the ions in such a way that their density is the same as the total charge density of the ions.The rightmost equation in Eq. ( 8) describes the run analysed in this paper which includes H + and Na + ions.In Eq. ( 9) the sum i is over all ion species, with the electric charge q i and it follows from the definition of the electric current density j .The electric current density is derived from Ampére's law: where the displacement current term is neglected.

Propagation of the magnetic field
The time evolution of magnetic field follows from Faraday's law.In this paper the electric conductivity is assumed to be infinite everywhere and, therefore, the magnetic field is propagated by the equation

Main algorithm
Equations ( 1), ( 2) and ( 11) are solved numerically as follows.The dynamical variables are the ion positions and velocities (x i , v i ) and the magnetic field B which is a grid quantity.The velocities are stored on whole time levels (n) and the positions and the magnetic field on half time levels (n − 1/2).The leapfrog algorithm to propagate the system from (x By the same pass, accumulate the density ) and the momentum density per unit mass p n from v n i .

Propagate the provisional (predictor step) magnetic field
B n * from B n−1/2 by using B n−1/2 and u n e * .It is provisional because B n−1/2 is used, which is not fully timecentred: 4. Compute the current density * /e). 5. Propagate the corrected (corrector step) magnetic field The electron flow u n e is not time-centred, but that cannot be avoided.
In the developed hybrid model the magnetic field can be solved by the leapfrog algorithm with a predictor-corrector scheme by doing steps 1-6 above, or without a predictorcorrector step by doing easily steps 1, 2, 5 and 6.The run analysed in the paper was made without the predictorcorrector but, instead, by using a small time step.
The particle mover is split into a position mover and a velocity mover.The positions and velocities are propagated by following a time-reversible second-order accurate Buneman scheme as follows (Hockney and Eastwood, 1988, p. 112).The position mover is, The velocity update formula where a = (q i /m i ) t and B n+1/2 must be evaluated at the particle position x n+1/2 i .By defining the vector ω = (1/2) a B n+1/2 , the update operation can be written in a slightly more compact way as 2.4 Advanced features: adaptive grid and particle splitting and joining The developed hybrid simulation model has many similarities with earlier models that have been used to study Mars and Venus (Brecht et al., 1993;Shimazu, 1999), except that the model includes an adapted hierarchically refined cubic grid and employs a particle splitting and joining technique.The use of an adapted grid (see Kallio and Janhunen, 2001, for an example of a spatially symmetric refined grid structure) gives a dramatic increase in performance, making it possible to run 3-D simulations of Mercury overnight in a personal computer.However, to realise this performance benefit in a particle simulation, one also has to allow for particle splitting and joining.Without particle splitting, the rather uninteresting undisturbed solar wind region would consume the most computational effort.
In the Mercury hybrid simulation the ions are split and joined when needed, to maintain about 10-20 ions per grid cell in all cells.There are many proposed methods for an adaptive control of the number of particles (Coppa et al., 1996).We adopted a simple method where each particle carries a weight w i , a number telling how many physical particles the macroparticle represents.Particle splitting is a relatively straightforward task because an ion is split into two by halving the weight and introducing a small random spatial displacement, so that the particle trajectories start to slowly diverge.Particle joining is a more complicated procedure because if two ions were joined to produce one ion, the momentum and the energy could not, in general, be conserved.Therefore, we select three ions for each joining process and produce two resultant ions.The three ions are selected among those in the grid cell that have the smallest velocity vector deviation squared, weighted by the ion weights w i .Of course, only ions of the same kind are joined.
The splitting algorithm works in detail as follows.The cell size is variable in the simulation box.There is a parameter N opt D (typically about 10) which specifies how many particles is "optimum" in a grid cell.One wants the number of particles in a cell to be in a certain range, 7...13, for instance.If there are fewer than 7 particles in a cell, one of them is split into two.The split one is the one with the largest weight.After splitting, the daughter particles are otherwise identical but have a small, randomly directed spatial displacement.The magnitude of the displacement is a certain fixed percentage of the local cell size, 10%, say.The displacement is perpendicular to the velocity of the particle(s).Such a randomly directed displacement vector can be generated by first generating any unit vector a which is linearly independent of the velocity vector v, and then generating another unit vector b = (v × a) 0 .A randomly directed unit vector which is perpendicular to v is then given by a cos φ + b sin φ, where φ is a random angle in the range 0...2π.Multiplying this by the desired magnitude gives the displacement.
The joining algorithm works, in turn, as follows.If n > N opt D is for some time, then joining takes places.To enable momentum and energy conservation, a join must produce at least two particles, so that, for example, the processes 3 → 2 and 4 → 2 are permissible but 2 → 1 is not.If we consider only the 3 → 2 process, the desired number of triples n is the same as the number of excess particles n − N opt D but truncated, i.e. n = min(n − N opt D , floor(n/3)).The triples should be selected so that they are close together.The smallest weights should be selected first for joining.
Consider now a cell which has an excess number of particles.We first find a pair of particles that minimises (w 1 + w 2 )(v 1 −v 2 ) 2 .We obviously also must require that the joined particles belong in the same species.When the optimal pair has been found, we look for a third partner that minimises w , where v CM is the centre of mass velocity of the three candidate particles.
When the three joined particles are found, we go to the centre of mass (CM) coordinates with respect to both position and velocity.We require that the momentum, kinetic energy, angular momentum and centre of mass position of the two new particles are the same as those of the deleted three old particles.The two particles have 12 degrees of freedom (each is described by a position and velocity vector) and the listed conservation laws fix 10 of these parameters.We fix the rest of the parameters by requiring that the centre of masses of both resulting particles are in the plane formed by the three original particles and that the velocity of the resultant particles (in the CM frame) is perpendicular to the angular momentum vector.This gives a pair of particles that is not too far from each other while still giving angular momentum conservation.In some rare cases one of the produced particles may fall outside the cell, even though the original three are inside.In such a case we cancel the whole joining step (this is not strictly necessary, but we do it for clarity).No random numbers are needed in the join process.
2.5 Input models: intrinsic magnetic field, planetary ions, solar wind and resistivity The hybrid model for Mercury uses the same equations and algorithms as its predecessor, the hybrid model for Mars (Kallio andJanhunen, 2001, 2002).The main difference between the Mercury and Mars hybrid models is that at Mercury the magnetic field includes a temporally constant magnetic field component caused by the Hermean intrinsic magnetic field.As a first approximation the intrinsic magnetic field is modelled by a magnetic dipole.The origin, orientation and the strength of the Hermean intrinsic dipolar magnetic field can be freely varied.It should be noted that it is also possible to include several magnetic dipoles in order to mimic possible magnetic multipoles.
The model includes two types of ions: the solar wind protons (H + ) and planetary ions.The density (n sw ), the bulk velocity (U SW = [U x , U z , U z ]), the temperature of the solar wind protons, as well as the strength and the orientation of the IMF can be varied freely.The model makes it possible to include several planetary ion species.At the moment, runs that included five planetary ion species have been tested.The planetary ions are originating from two sources: (1) from the exosphere ("a corona source") and (2) from the surface ("a surface source").Ions from the exosphere can be produced by photoionization and by charge exchange processes between ions and planetary neutrals.The exospheric density profile can be an arbitrary 3-D function.Ions from the surface can be emitted according to a 2-D (any latitudinal and longitudinal dependence) ion emission model.The surface source is used to model the ion production from the exosphere in cases when the neutral density scale height is small compared to the grid size, i.e. of the order of tens of kilometres, as well as to mimic the emission of ions from the surface itself caused by particle precipitation and solar photons.
Finally, the model makes it possible to use any 3-D electric conductivity model, if needed.However, it is not evident which conductivity value should be most appropriate at different sites outside and within the planet.For example, one may anticipate that adding resistivity to the magnetopause and to the cross tail current sheet can result in more intensive reconnection than without a resistivity.It is also noteworthy that the model does include certain numerical resistivity and, consequently, reconnection, although the resistivity is put to zero in Eq. ( 6).Therefore, in the run presented in this paper, as in the majority of the test runs, the conductivity has been assumed to be infinite, i.e. the resistivity η (= 1/σ ) is put to zero both within the planet and outside it.

Boundary conditions
The model includes two types of boundaries: the six faces of the simulation box and the surface of the planet.The solar wind protons were launched at a face x = const.An ion is taken away from the simulation if it hits the simulation box.The boundary conditions for the planetary body provide a more challenging problem.The surface can act as a source for planetary ions by the implemented surface source model, as pointed out before.In the analysed hybrid model run the surface is treated as a fully absorbing surface for the ions, i.e. an ion is taken away from the simulation if its radial distance becomes smaller than the radius of the planet.The density and velocity of ions within the planet are therefore zero.Furthermore, the velocity of electrons is also put to zero inside the planet.

Zero-divergence magnetic field and conservation laws
The magnetic field B is stored in the hybrid model on cell faces, which preserves ∇ • B = 0 within a round-off error when B is propagated by Faraday's law (see Dai and Woodward, 1998).The developed model does not, therefore, include unphysical magnetic sinks or sources (monopoles).
In the hybrid model mass is automatically conserved because the ions are modelled as particles.However, there is nothing that could automatically cause an exact conservation of momentum and energy within the simulation box.The motion of the macroparticles and, consequently, also the electric currents and magnetic fields (Eq.11), are changed by the Lorentz force (Eq.2), but nothing forces the change in the total momentum or energy to remain unchanged.The hybrid model may, therefore, include non-physical energy sinks or sources.Therefore, the runs have to be carefully diagnosed, in order to make a optimum choice of the free numerical parameters, especially the time step.For example, in the analysed run the total energy was monitored, in order to make sure that the total energy in the simulation box does not increase faster than the energy which flows into the box from the solar wind, and that the total energy remains practically unchanged during the time (100 s < t < 320 s) when the solution presented in this paper was analysed.

Limitations of the hybrid model approach
Although the developed hybrid model has been shown to provide a powerful tool to study non-magnetic planet Mars (Kallio andJanhunen, 2001, 2002), it is not self-evident without extensive testing that the same approach is valid for the magnetised planet Mercury both because of (1) the approximations used and (2) the numerical algorithm.

Excluded physical processes
The quasi-neutral hybrid model assumes quasi-neutrality of the electric charge.Such an assumption is valid if the spatial scale L of the studied physical process is larger than the plasma Debye length λ D (≡ ε 0 k T /(n e e 2 )), i.e.L > λ D .For example, L > λ D limits the approach to model plasma sheath effects near the planet which may be caused by photoelectrons (see discussion in Grard, 1997).The model also assumes massless electrons (m e = 0), and, therefore, the approach does not include electron particle effects or waves which depend on the electron inertia.Therefore, the approach cannot describe finite electron gyroradius effects that may play a role, say, at the cross tail current sheet.Waveelectron interactions also cannot be modelled by the quasineutral hybrid model.
It should also be noted that 3-D quasi-neutral hybrid model is typically operated with the time step of the order of the ion gyroperiod, which restricts its possibility to include highfrequency phenomena.The hybrid model is also computationally expensive which limits the possibly to study potentially slowly evolving processes with a 3-D quasi-neutral hybrid model.How these assumptions and limitations restrict the approach to describe basic the features of the global Mercury-solar wind interaction depends on the importance of the excluded processes.

Numerical scheme and implemented constraints
As noted before, the accuracy of the conservation of the total momentum and energy in a quasi-neutral hybrid model has to be studied case by case.A non-energy conservative scheme may be numerically highly unstable because a possible numerical instability has potentially unlimited energy growth.
Building up the Mercury-solar wind interaction by launching solar wind protons into the simulation box provides an example of a period that needs extra caution, to avoid the sudden impact of the solar wind to the planetary magnetic field, which would trigger a numerical instability.In reality, such an instability can be damped over the course of time, but in the hybrid simulation it may continue to grow without limit due to the artificial energy sources.In particular, Eq. ( 9) shows that if the plasma density, n e , drops temporarily to a low value (which may take place simply due to the natural temporal variations), then the electron velocity U e and, consequently, the electric field, can be suddenly increased.Such a strong electric field situation can cause a sudden acceleration of ions (Eq.2) and the growth of the magnetic field (Eq.11), which, in turn, can result in an increase in the Lorentz force (v i × B) and electric current, etc.
There are many ways to restrict the possible increase in the energy.One may include numerical diffusion in Eq. ( 11) that introduces spatial smoothing of the magnetic field.The increase in the magnetic energy may also be controlled by introducing a minimum plasma density, n min , in Eq. ( 9), which will be used if n e < n min .One may also introduce a maxi-mum electron velocity value in Eq. ( 9), U max e , which will be used for the electron velocity if U e > U max e in a grid cell, or some other restriction can be employed.One can also use the solution obtained from GD or MHD models as a initial state from which the hybrid simulation is initiated (Shimazu, 1999;Terada et al., 2002).
At the run presented in this paper below, the smooth growth of the system is obtained by restricting the magnitude of the electric field E in Faraday's law (Eq.11) to a maximum electric field value, E max , where |E| will not be exceeded in Faraday's law, i.e. |E| = min(|E|,E max ) in Eq. ( 11).The smaller the E max value is used, the more slowly the initial dipolar magnetic field is deformed by the solar wind.It should be noted that E max is used only in Faraday's law when B is propagated (Eq.11) but that there are no restrictions for the electric field in the Lorentz force (Eq.2).Therefore, E max does not cause explicit constraints for the acceleration of ions, although it can affect the acceleration of ions indirectly by restricting the time variability of B and, consequently, also the magnitude of the v i × B force.Another restriction on E in Faraday's law is that E is set to zero, if the total ion density is smaller than a predetermined minimum value, n min .Such a constraint was introduced in order to avoid electric field fluctuations and possible sudden particle accelerations, if the total ion density in the simulation becomes low.
To summarise, in the analysed run presented in Sect. 3 the interaction between the planet and the solar wind was forced to develop smoothly, while the solar wind protons fill the simulation box by introducing two constraints: E max of 0.001 V m −1 and n min of 0.1 cm −3 .In the analysed run the maximum magnetic field growth rate, max( B/ t), can thus be estimated to be ∼E max /L ∼ (10 −3 V m −1 )/ 305 10 3 m = 3.3 nT s −1 , where L is the grid size used, 305 km, and, consequently, the total magnetic field increase during the 320-s long run is estimated to be ∼1000 nT.The maximum rate of current density increase, max( j/ t), can similarly be estimated to be ∼ (L µ o ) −1 max( B/ t) ∼ 8.6 × 10 −9 Am −2 s −1 , and the possible maximum current density at the end of the run is estimated to be ∼ 3 × 10 −6 Am −2 .
The constraints used can be anticipated to cause some artificial smoothness in the growth of the magnetic field and electric currents.These constraints may especially affect the development of the magnetic field in the magnetic lobes where the particle density is low.The analysed run may nonetheless be anticipated to improve our insight into the properties of the plasmas and fields near Mercury, because the solution qualitatively resembles many runs using different constraints.

Applications and advantages of the developed hybrid model approach to study Mercury
The applicability of the hybrid model approach to study Mercury results mostly from the small size of the planet and its relatively weak intrinsic magnetic field.The spatial size of the Hermean interaction region with the solar wind is comparable to or smaller than the interaction regions formed at Mars, Venus and Titan, which have been intensively studied by quasi-neutral hybrid models (Brecht and Ferrante, 1991;Shimazu, 1999;Kallio and Janhunen, 2001;Brecht et al., 2000).The absolute number of macroparticles and the grid cells used for Mars (Kallio andJanhunen, 2001, 2002) may thus be anticipated to be sufficient to obtain adequate spatial resolution for Mercury as well.
The weakness of the Hermean intrinsic magnetic field has several implications, which make it an object suitable for a hybrid model approach.First, the CFL (Courant-Friedrich-Levy) conditions can be easily fulfilled.This condition states that the value of the grid size, L, divided by the time step, T step , must be larger than the maximum speed of information, c max , in the simulation, L/ T step ≤ c max , for stability reasons.In practice, the test runs suggest that the time step T step of ∼0.01 s is small enough to obtain a stable solution with a good spatial resolution (a few hundred kilometre) at Mercury.It should also be noted that a time step of ∼0.01 s is also much smaller than the minimum gyrotime of protons of T c (H + ) ∼0.1 s (= gyrotime in a 600 nT magnetic field) and, therefore, the particle mover scheme (Eq.12) follows the trajectory of ions accurately.No such large time step could have been used if the intrinsic magnetic field was much stronger, such as in the case near the Earth, since the strong magnetic field would have yielded a high wave velocities and small gyrotimes.
Another time scale of importance is Lx box /U SW , where Lx box is the length of the simulation box in the x direction.The resulting time, that could be referred to as a filling time, T fill , gives an estimation for the time scale for solar wind protons to fill the simulation box.One may state that the parameters in the hybrid simulation should be chosen in such a way that the simulation could be run at least a few times T fill , in order to give the system enough time to stabilise and the resulting physical processes to evolve.In the analysed run, T fill = 8R M /430 km s −1 = 45 s and, thus, the run time of 320 s is seven times T fill .
From the computational cost point of view, Mercury is thus an object that is anticipated to be well suited for even 3-D hybrid simulation study, with a good spatial resolution.Furthermore, from the physical point of view, the relatively weak intrinsic magnetic field and spatially limited magnetosphere are anticipated to result in an environment at which a finite gyroradius and other kinetic effects may start to play an important role, i.e. especially for the case when heavy planetary ions are being considered (see, for example, Ip 1987), for several reasons.First, the bulk motion of different ion species may differ noticeably from each other.The velocity distribution may also be anticipated to be non-Maxwellian because of the large ion loss cone.A quasi-neutral hybrid model can be used to study all of these effects because ion species move independently from each other and no assumptions of a Maxwellian velocity distribution or similar temperature or bulk velocities of different ion species have been made.Also, in the hybrid model ions can precipitate onto the Fig. 1.Morphology of the magnetic field in the quasi-neutral hybrid model in a pure northward IMF case.Magnetic field lines are connected to Mercury (red lines) and to the solar wind (blue lines).The magnetic field lines are calculated by using the following starting points: 20 starting points along two lines which go via the points (1.1,0,0) R M and (−1,0,+−4) R M and 30 points along the two lines via the points (−1,0,0) R M and (−1,0,+−4) R M .Note that, although the starting points are on the XZ plane the magnetic field lines are not on the XZ plane.The small square on the left bottom corner shows the size of the grid which was used at the field line tracing.
planetary surface.Ions can also be emitted from the surface by a given temperature model.The implemented planetary ion sources may thus be anticipated to give valuable indirect information about various physical processes and questions possibly associated with kinetic effects, such as the role of planetary ions for magnetospheric dynamics (Ip, 1987), the electric conductivity outside and inside of the planet and, thus, the current system (Glassmeier, 2000), as well as particle sputtering (Ip, 1997;Killen et al., 2001).
It should finally be noted that the test runs made without splitting and joining techniques have shown that such a technique has to be used in order to obtain a solution in the nightside where the plasma density can be several orders of magnitude smaller than in the solar wind.One could not have been able to obtain good spatial resolution in the nightside by using a constant weight w i for all ions, without substantially increasing the number of ions per grid cell.
3 Case Study: a closed magnetosphere

Input parameters and the used constraints
Figures 1-6 give an example of the plasma and field parameters resulting from a run in which the following solar wind parameters were chosen to mimic the situation at Mercury in its pericenter: n = 76 cm −3 and U SW = [−430, 0, 0] km s −1 .
The magnitude of IMF at the pericenter at the 0.31 Astronomical Unit (AU) has been estimated to be about 46 nT and the spiral angle, the angle between the IMF and the Mercury-Sun line, to be about 17 • (Slavin and Holzer, 1981).In this paper, however, a "closed" magnetosphere case was studied by adopting an IMF that was directed along +z-axis: B IMF = [0, 0, 10] nT.Such an IMF gives a rough estimation for the IMF transverse component (the IMF component perpendicular to the solar wind velocity vector) without the IMF x-component.This choice has been made to study a case where the magnetosphere does not include a noticeable north-south asymmetry caused by a large IMF x-component (see Kallio and Janhunen, 2003, about the effects of the IMF x-component to the Hermean magnetosphere in the quasineutral hybrid model).One consequence of the chosen IMF is that in the analyzed run the Alfvén Mach number, M A , is 17.Finally, in the simulation the solar wind protons launched into the simulation box drawn from a Maxwellian velocity distribution have a thermal energy kT = 19 eV.If one adopts a polytropic index, γ , of 5/3, the flow corresponds to sonic Mach number, M s (= U SW / √ γ kT /m) of 7.7.The electron temperature has been assumed to be zero everywhere (cf.Eq. 7).Thus, the analyzed case represents more of a superalfénic and supersonic case, as is the case in reality, because the Mach numbers at the pericenter have been estimated to be M A = 3.9 and M s = 5.5 (Slavin and Holzer, 1981).
The run included two simple, spherically symmetric, exponentially decreasing neutral density profiles with artificially chosen scale heights of 1500 km and 400 km.The analysed run includes only one planetary ion species with an m/q ratio of 23 amu/e.It should be noted that although the m/q ratio used is the same as Na ions, the implemented ions are not strictly speaking supposed to represent only Na ions, but instead, to mimic, at least approximately, the possible effects of various planetary neutrals.The total photoionisation rate was chosen to be 10 24 s −1 .No charge exchange processes were included.The model includes a 64 × 64 × 64 equal size grid cell of 305 km.The time step was 0.005 s.The average number of ions in a grid cell was chosen to be 15, i.e. about 4 million ions in total.Resistivity η was put to zero, both outside and inside of the planet.The smooth growth of the system (see discussion in Sect.2.8.2) was obtained by using E max of 1 mV m −1 and n min of 0.1 cm −3 .The intrinsic Hermean magnetic field was modelled by a magnetic dipole moment aligned with the −z axis, with field [0,0,300] nT at the magnetic equator.The size of the simulation box was −4 × R M < x, y, z < 4 × R M (R M = 2440 km).The computation time corresponds to 320 s in real time.The plasma parameters shown in Figs.3-6 show mean values calculated from ten instantaneous values at 10-s intervals from physical time 160 to 250 s, in order to illustrate average plasma parameters.
3.2 Models for the shape of the bowshock and the magnetopause The plasma and field parameters shown below are superimposed to 2 curves referred to as the bowshock and the magnetopause in the paper.The curves are obtained by using empirical Earth bowshock and magnetopause shapes.The two curves are included in order to help the eye to see the possible asymmetries, to see how the changes in the plasma parameters are related to the overall interaction regions, and to help to make a comparison between different plasma parameters.
The shape of the bowshock is modelled by the formula (Slavin and Holzer, 1981) Here, r is the radial distance from the focus point, r focus to the bowshock, e is the eccentricity, is the angle between the a point of a bowshock and the x-axis and L is the distance of the bowshock from the focus point on the x = x o plane.The terrestrial bowshock can be modelled by choosing e = 1.16,L = 23.3R E and x o = 3 R E (Slavin and Holzer, 1981).In the analysed Mercury run Eq.( 20) was found to give a relatively good approximation for the average position of the bowshock, if L and x o were scaled down by a factor of 6.2.The values at Mercury used (L = 3.76 R M and x o = 0.48 R M ) correspond to the sub-solar bowshock distance of 2.22 R M .
The magnetopause is modelled by an empirical formula (Shue et al., 1997) where r o is the subsolar distance from the centre of the planet, is the angle from the x-axis and α is the level of flaring.Modification of the shape of the Earth's magnetopause caused by "erosion" or "reconnection" have been modelled by using various functional forms for r o = r o (P dyn., B z ) and α = α(P dyn., B z ), where P dyn. is the dynamic pressure for the solar wind and B z is the IMF zcomponent (see references in Yang et al., 2002).Earlier, the possible role of erosion has been considered in regards to the applicability of Mariner 10 nightside magnetopause crossings for the prediction of the subsolar distance of the Hermean magnetopause (Slavin and Holzer, 1979).In the analysed hybrid simulation run, r o and α were chosen to provide an approximation for the average shape of the magnetopause at the analysed run at time, t = 160 − 250 s: r o = 1.37 R M and α = 0.8.

Plasma and magnetic field in the hybrid model
Figure 1 illustrates how the initially dipolar magnetic field is deformed due to the solar wind in the analysed hybrid model run.The position of the magnetic tail lobes, the cross tail current sheet and the magnetic cusps are clearly seen.Although not seen in Fig. 1, it should be noted that although the IMF is purely northward, there is a dawn-dusk asymmetry in the magnetic field.Such an asymmetry would not have arisen in an MHD model for the same input parameters.The dawn-dusk symmetry can be broken by the convective electric field, which in the solar wind points to the −y direction, and by finite ion gyroradius effects.Due to the asymmetry, the magnetic field lines shown in Fig. 1 are not entirely in the XZ plane.They also have a non-zero B y component, although the points from which the field line tracing were started are on the XZ plane.
Figure 2 shows the normalised magnetic field vectors [B x , B z ]/ B 2 x + B 2 y + B 2 z on the noon-midnight meridian plane.Although a detailed study of the morphology is beyond the scope of the present study, such a presentation of the magnetic field makes possible a more quantitative analysis of the field properties than the 3-D field lines given in Fig. 1.The solid blue line in Fig. 2 illustrates a boundary at which the magnetic field changes its direction noticeably near the terminator plane and in the nightside.Such a change in direction is expected to be formed in a closed magnetosphere case at the magnetopause behind the magnetic cusps where the draped magnetic field in the magnetosheath is almost antiparallel to the magnetospheric magnetic field.The site of the sudden changes in the direction of B in Fig. 2 can thus be interpreted as the position of the magnetopause.The shape of the transition region can be approximated by Eq. ( 21), using the subsolar distance, r o , of 1.37 R M and α = 0.8, as noted in Sect.3.2.
The density of the solar wind protons near the planet on the magnetic equator and on the noon-midnight meridian plane is shown in Fig. 3.When the solar wind reaches the planet, the density of the solar wind protons first increases, but then starts to decrease closer to the planet and on the nightside.The maximum proton density is about four times larger than the undisturbed density in the solar wind.The position where the increase takes place is modelled in Fig. 3 by a line which is derived from Eq. ( 20) by using e = 1.16,L = 3.76 R M and x o = 0.48 R M .Note that the density is not axially symmetric around the x-axis, but is higher in the magnetosheath on the XZ plane than on the XY plane.The dependence of the solar zenith angle (SZA) is also different on these two planes: on the XY plane (Fig. 3, top) the density decreases with increasing SZA in the magnetosheath, but on the XZ plane (Fig. 3, bottom) the maximum densities are highest close to the magnetic cusps, i.e. at SZA ∼65 • (see Fig. 1).A slight dawn-to-dusk asymmetry can also be seen on the XY plane, with the bowshock being farther away from the planet on the duskside than on the dawnside.
Figure 4 shows in more detail the low solar wind proton density regions.On the magnetic equator the density is lowest near the planet on the nightside (Fig. 4,top ).The density of the solar wind protons in the magnetosphere is typically ∼0.1-1 cm −3 , i.e. over a hundred times smaller than upstream of the bowshock of 76 cm −3 .The proton density is slightly asymmetric in the magnetosphere, with the density near the planet being slightly smaller in the dusk than in the dawn side.On the noon-midnight meridian plane the most tenuous proton regions are within the magnetic tail lobes (Fig. 4, bottom ).The increase in the density near the planet at the magnetic cusps is also clearly seen.
The bulk speed of the antisunward solar wind protons, |U x |, first decreases at the bowshock, but then starts to increase in the magnetosheath with increasing SZA (Fig. 5).The maximum decrease in the solar wind speed near the xaxis is about 1/4 of its speed in the solar wind.The flow in the magnetosphere is mostly antisunward (U x < 0).Positive U x regions can be found everywhere around the planet on the XY -plane.On the XZ plane these U x < 0 regions are located on the nighside near the magnetic lobes.Finally, Fig. 6 illustrates how the bulk velocity has small U y and U z components near the x-axis on the dayside, but that these transverse velocity components increase in the magnetosheath with increasing SZA. and analysed.Although the details of the solution have been found to vary from test run to test run, the observed physically relevant basic features in these runs imply that the quasi-neutral hybrid model can be used to study effectively many important features near Mercury, such as the formation and characteristics of the bowshock, magnetopause and the magnetotail.Near future plans are to use the developed model to study the precipitation of solar wind protons and the planetary ions to the surface, properties of the planetary ions near the planet, the role of the resistivity of the planet and to study how the upstream parameters control properties of the Hermean magnetosphere.

Summary
The quasi-neutral hybrid model provides many unique possibilities to study self-consistently how the solar wind interacts with a small planetary body, such as the role of finite gyroradius effects, different ion species and plasma-surface interaction processes.The approach is anticipated to be especially well suited for studying weakly magnetised bodies, whose overall spatial size is much smaller than the size of the terrestrial magnetosphere, because finite gyroradius effects then start to play an important role.Nowadays the computational resources make it possible to conduct 3-D quasi-neutral hybrid model runs with a relatively good spatial resolution.
The paper represents the first attempt to use a global 3-D quasi-neutral hybrid model to study how the solar wind interacts with Mercury.The analysed run illustrates how a bowshock and a magnetopause are formed in front of the planet, as well as the formation of the Hermean magnetotail.In general, although much more testing has to be done in the future to analyse the role of different upstream conditions, planetary neutral models, as well as numerical schemes, the analysis suggests the potential of the quasi-neutral hybrid model approach to address several important plasma physical processes near Mercury.

Fig. 2 .
Fig. 2. Normalised magnetic field B x and B z components on the XY -plane: [B x , B z ]/ B 2 x + B 2 y + B 2 z .The blue solid line shows the boundary where the direction of the magnetic field changes suddenly.The run is the same as shown in Fig. 1.