Space weathering on the Moon: Farside-nearside solar wind precipitation asymmetry Planetary and Space Science

The lunar surface is continuously under the impact of solar wind plasma, which breaks the chemical bonds of the surface material resulting in weathering of the surface and a modi ﬁ ed chemical composition. Ion impact also sputters the surface material, affecting the composition of the lunar exosphere, and it controls the electrical properties of both the lunar surface and the near surface space. We have studied the lunar farside-nearside (FN) hemispherical asymmetry of the solar wind proton impact on the lunar surface along the orbit of the Moon during fast solar wind conditions. The analysis is based on a 3D hybrid model where ions are accelerated by the macroscopic j (cid:1) B and pressure gradient forces. The derived proton impact surface map shows that the highest cumulative solar wind proton addition on the lunar surface is located on the farside while the most energetic protons precipitate on the nearside. The total ion impact rate was found to be smallest when the Moon is deep in the magnetotail. The total ion impact rate on the lunar surface varies while the Moon orbits the Earth and these longitudinal variations are caused by the magnetosphere and lunar tidal locking.


Introduction
Lunar surface is continuously bombarded by the solar wind plasma because the Moon has neither an atmosphere nor a global intrinsic magnetic field that could deflect the plasma flow around it. The impacting solar wind particles modify the properties of the surface and this process is known as space weathering. Recent observations of the optical properties of walls of lunar craters has shown, that there is an East-West asymmetry inside craters, which has been suggested to be consequence of systematic change of the properties of plasma when the Moon is in the magnetosphere (Sim et al., 2017). In addition to these longitudinal space weather effect variations, there are also latitudinal variations, which has been suggesting to be related the decreasing intensity of precipitating solar wind particles with increasing latitude (Hemingway et al., 2015). However, global latitudinal and longitudinal variations in space weathering on the lunar surface has not yet studied in detail with a global plasma simulation.
It is foreseen that the properties on the impacting charged particles on the lunar surface have both long term and short-term variations because the solar wind varies in time over all time scales. Moreover, the Moon spends about one fourth of its time downstream of the Earth's bow shock in the magnetosheath and in the magnetotail where the properties of the plasma are different to the properties in the undisturbed solar wind. In the Earth's magnetosheath, both on the dusk and dawn flanks of the magnetosphere, the Moon is exposed to shocked (i.e., subsonic and highly thermalized) solar wind plasma. Due to the tidal locking of the Moon, the precipitation of the undisturbed solar wind plasma and the precipitation of magnetospheric plasma are not hemispherically identical where, on average, the lunar farside faces more frequently toward the undisturbed solar wind than the nearside. Because of the tidal locking, the cumulative total proton precipitation rate on the lunar surface is different on opposite sides of the Moon and the instantaneous plasma precipitation depends on the position of the Moon on its orbit.
One important individual question when the lunar space weathering is analysed is what is the detailed 3D ion velocity distribution function of the impacting solar wind protons. This is an important issue because the ion velocity distribution has several significant consequences on the lunar surface, near the surface particles and electrodynamics and near the lunar space. First, the impacting solar wind protons (H þ ), alpha particles (He þþ ) and multiple charged heavy ions (e.g. O 6þ ) break the chemical bonds of the lunar surface material. This results in formation of new chemical minerals and molecules, such as hydroxyl (Liu et al., 2012). Impacting multiple charged heavy solar wind ions in turn results in soft x-ray emission (see e.g. Kallio et al., 2008 and references therein). The impact of plasma is also locally affected by lunar magnetic anomalies where the Lorentz force deflects the flow the solar wind charged particles (see e.g. Vorburger et al., 2012;Jarvinen et al., 2014) and, consequently, results in non-isotropic particle precipitation and non-isotropic weathering of the surface near a magnetic anomaly. This process has triggered interest regarding the origin of lunar swirls because impacting particles darken the surface by chemically weathering where ion impact results in the formation of nanophase iron (e.g. Hood and Schubert, 1980;Blewett et al., 2007Blewett et al., , 2011Kramer et al., 2011;Lucey et al., 2017).
Second, the impacting solar wind electrons and ions together with the photoelectrons control lunar surface charging and, consequently, the electrical properties of the near surface space. In particular, the density and temperature of the plasma affect the extent and strength of the Debye layer (e.g. Dyadechkin et al., 2015). The electric field near the lunar surface in turn affects the motion of charged particles near the surface, such as electrically charged lunar dust particles (Nitter et al., 1998). The lunar dust, in turn, can be a risk for the manned missions, both for human beings on the Moon and the technical devices on the surface (see e.g. Kallio et al., 2012).
Third, the impacting ions, mostly protons and alpha particles, release material from the lunar surface by the sputter process. This sputtering is one of the surface releases processes which form the lunar exosphere (Wilson et al., 2006;Wurz et al., 2007) as well as the exospheres of other so-called airless planetary objects, such as Mercury and asteroids. Moreover, part of the impacting protons bounce off the lunar surface as energetic neutral atoms, ENAs, which have been used as a remote sensing method on plasma properties near lunar magnetic anomalies , Vorburger et al., 2013Harada et al., 2014).
The aforementioned processes are energy dependent and they also depend on the direction of the ion velocity vector in relation of the lunar surface, i.e., on the 3D proton velocity distribution function. A difference has been observed in the plasma interaction with the lunar surface using ENAs in the magnetosheath (Allegrini et al., 2013;Lue et al., 2016) and in the plasma sheet (Harada et al., 2014). Therefore, more information about the 3D proton velocity distribution function will increase our understanding of the role and consequences of several important plasma physical processes near the Moon.
To better understand the properties of instantaneous proton precipitation at different places on the lunar surface raises another important question on the cumulative effect of the proton precipitation on the lunar surface, that is, how the fact that the Moon is tidally locked and that it spends about one third of its time in plasma regions which are disturbed by the Earth's magnetic field? For example, the surface of Jupiter's moon Ganymede shows clear global surface variations correlated to its magnetic field, which is a result of a cumulative non-isotropic ion precipitation of Jupiter's magnetosphere plasma (Khurana et al., 2007).
Modelling of the properties of the plasma near the Moon is, therefore, an important scientific task. So far, the most commonly used global magnetosphere models are magnetohydrodynamic (MHD) models (see e.g. Harnett et al., 2013Harnett et al., , T oth et al., 2005, and references therein) where the charged plasma is modelled as a fluid. However, the limitation of the MHD is that it assumes a Maxwellian velocity distribution function. Instead, simulating a general 3D ion velocity distribution (non--Maxwellian), is a challenging task because it requires a kinetic 3D magnetosphere model, which is computationally highly expensive without making some model approximations and simplifications in the simulation.
One frequently used simulation method is so-called hybrid simulation, where ions are modelled as particles while electrons form a massless charge-neutralising fluid. In this approach ions are accelerated by the Lorentz force and there are no prior assumptions about the velocity distribution functions, for example, that it is Maxwellian. Such an approach has been used widely to study different solar System objects such as Earth's bow shock and magnetosphere (e.g. Winske, 1985;Swift, 1995Swift, , 2004Omidi et al., 1998;Winske et al., 2003;Lin and Wang, 2005;Karimabadi et al., 2006). Moreover, the approach has also been used to investigate the small-sized Hermean magnetosphere and precipitation of solar wind ions on the planet's surface as well as to the lunar surface when it is in the undisturbed solar wind (Kallio et al., 2008). Furthermore, Earth's magnetosphere has been simulated by global full kinetic Particle-in-Cell (PIC) simulations where also electrons are modelled as charged particles (see e.g. Cai et al., 2015).
Increased computational capacity has made it possible to make global and local simulations where MHD is combined with kinetic models, such as PIC (e.g. Daldorff et al., 2014), hybrid (see e.g. Lipatov andSibeck, 2016, Hesse et al., 2016 andref. therein) and full kinetic PIC simulations (e.g. Peng et al., 2015). The large spatial size of the Earth's magnetosphere and the large distance of the Moon from the Earth is, however, a highly challenging computational task for kinetic models. A highly computationally extensive model, in turn, limits usage to analysis of routinely numerous different solar wind conditions. Therefore, the usage of the kinetic approach is limited to providing information for studies which require information of the 3D ion velocity distribution function.
This study uses a hybrid model where ions are accelerated by the macroscopic j Â B and pressure gradient forces, as in the MHD model, and not by the Lorentz force. The advantage of the approach is that such a hybrid model, which is referred to here as the j Â B hybrid model to distinguish it from the typical hybrid model, is that the simulation time step can be much larger compared to the hybrid model because there is no need to follow the ion gyromotion. This makes the j Â B hybrid model computationally economical. The limitation of the approach is that it does not provide information about an ion gyromotion around the magnetic field caused by the Lorentz force. However, the advantage of the j Â B hybrid model compared to a typical MHD model is that it does not assume a Maxwellian velocity distribution function and, therefore, it provides a full 3D ion velocity distribution function.
The paper is organised as follows. First, the j Â B hybrid model is presented and its pros and cons are discussed. Then the basic properties on the resulted magnetosphere are presented. After that, the properties of the precipitating solar wind ions are investigated. Finally, the value of the approach for further lunar research is discussed.

Description of the j £ B hybrid model
The self-consistent 3D j Â B hybrid model, which is used to study ion precipitation on the lunar surface, is a kinetic model where protons are modelled as a particle cloud while electrons form a massless chargeneutralising fluid. A particle cloud presents a large number of real particles. The cloud, or a parcel of a fluid, has the following properties: the mass of the ion (m i ) of which the cloud is formed, which in our case is the mass of a proton, the electric charge (q i ) of ions within the cloud, in our case the unit charge e, and the weight (w i ) which represents the number of ions which are within a cloud.
The 3D simulation box is divided into a mesh with a constant cubical Cartesian grid of the grid cell length dl. The macroscopic plasma parameters, the density of protons (n) and the bulk velocity of protons (U) are obtained by accumulating the particles within the cloud, the size of which is the same as the size of the simulation cell (dl 3 ), into the grid. The magnetic field (B) is propagated in time from the electric field (E) by using Faraday's law The electric field was obtained from the generalized Ohm's law where p ¼ n k B T e (k B is the Boltzmann's constant, T e is the electron temperature) is the electron thermal pressure. In this paper the electron temperature has been assumed to be a constant, as commonly assumed in hybrid models. In Eq (2) the electron bulk velocity (u e ) is derived from the electric current and Amp ere's law where the plasma has been assumed to be quasi-neutral, that is, the particle density of electrons (n e ) has been assumed to be equal to the particle density of ions (n). It should be noted that Eqs.
(2) and (3b) state that the electric field includes the Hall term (e j Â B). Finally, the equation of motion in the model is where f i is the force density that accelerates the ion cloud i. Equations (1)-(4) and the used numerical algorithms in the j Â B hybrid model are identical to the common hybrid model which has been used to simulate the properties of solar wind plasma near different solar wind objects, Mercury, Venus, the Moon, Mars, an asteroid, a comet and the solar wind (Kallio and Janhunen, 2003;Kallio 2005;Kallio et al., 2008;Dyadechkin et al., 2017) i.e., if we assume that the ions are accelerated by the Lorentz force. The hybrid model is a practical approach if the size of the simulation box is so small that ions can travel through the box with relatively few time steps (dt) and if the maximum magnetic field is sufficiently small so that a relatively large time step can be used to follow the ion gyromotion around the magnetic field.
However, as mentioned earlier, a 3D hybrid model of Earth's magnetosphere is computationally highly expensive. Therefore, the most commonly used approach to simulate the Earth's magnetosphere is a single ion specie Magnetohydrodynamic (MHD) model where the plasma is modelled as a fluid and where the magnetic force is the j Â B-force. This approach does not include the ion gyrotime scale as in the hybrid model and, consequently, a large time step can be used which speeds up computations.
The disadvantage of the MHD approach, however, is that it does not include kinetic effects because the velocity distribution is assumed to be Maxwellian. If the length scales of interest are larger than the gyroradius, this does not limit the applicability of the MHD approach and, in fact, it is often the only possibility for making a 3D simulation of large and highly magnetised regions, like the solar corona or the magnetospheres of the giant planets Jupiter, Saturn, Uranus and Neptune.
However, the size of the Moon can be relatively large compared with the gyroradius of protons when the Moon is within the centre of the plasma sheet where the temperature of the plasma is millions of Kelvins (Wing and Newell, 2002). From the simulation point of view fast solar wind cases are especially interesting for a kinetic model because the protons are then hotter than in the slow solar wind (e.g. Ebert et al., 2009). A global magnetosphere simulation that covers at least certain kinetics effects can, therefore, be expected to provide a more realistic description of the properties around the Moon, as well as of the properties of the ions which impact on the lunar surface resulting in various physical and chemical space weathering processes on the surface.
In this paper a 3D kinetic magnetosphere is used by adapting the previously used hybrid model (see Kallio and Janhunen, 2003;Kallio, 2005). Rather than using the Lorentz force, the force density was assumed to be the typical j Â B force and the pressure gradient as in an MHD: Therefore the equation of the motion of an ion cloud is where the value of the ion mass density ρ i ( m i ðw i =dl 3 Þ) is the same as the mass density of a cloud. It is informative to note that Eq. (6a) has a similar mathematical form as the MHD model where magnetised fluid is accelerated by the j Â B force and the pressure gradient force ρðr; tÞ where ρðr; tÞ is a scalar field that gives the mass density at a given position (r) at a given time (t). Eq. (6a) and Eq. (7) have, however, the fundamental difference that in the former the velocity is the velocity of the cloud in a Lagrangian (moving) frame and the mass density of the cloud is a constant, while, in the latter, the mass density is a time and position dependent scalar field and the bulk velocity is a time and position dependent vector field. It is also illustrative to note that the equation of motion in a hybrid model can also be expressed in a form that resembles Eq. (7) by writing is the mass density of the macroparticle particle and δn i ( w i =dl 3 ) is the number density of the particle of weight w i . The electric current δj i e δn i ðv i À u e Þ can be interpreted as the part of the total electric current j which is associated with the individual macroparticle i of the weight w i and the associated same number w i of electrons from the total electron fluid. Similarly, the pressure force δrp i ðδn i =nÞrp can be interpret to be the part of the total pressure force associated with the individual macroparticle i. It should be, however, noted that although Eq. (8b) has formal similarities with Eq. (7), it is still the common Newton's 2nd law in a hybrid model which is only reformulated. In particular, an ion having an equation of motion as in Eq. (8b) makes a gyromotion around the magnetic field and, consequently, has a characteristic ion gyroperiod time, in contrast to a cloud with an equation of motion as in Eq. (7). The advantage of the j Â B hybrid model is that the ion velocity distribution can be fully 3D as is the case in the hybrid model. For example, it is anticipated that the velocity distribution function can have different temperatures along the magnetic field than perpendicular to it because the j Â B force causes acceleration perpendicular to the magnetic field. However, the j Â B hybrid model has the same challenge as the hybrid model in terms of conservation laws: there is no guarantee that the total energy and momentum are conserved within the simulation box, unlike in a single fluid MHD model which is based on the conservation of momentum and energy. Therefore, a j Â B hybrid model run is in that sense a numerical experiment where the validity and usefulness of the run is not known beforehand, but its applicability must be investigated by making the simulations and carefully analysing their results.
In this paper, the detailed analysis of the ion precipitation to the Moon's surface was made by adding 36 "virtual Moon detectors" placed uniformly along the lunar orbit. They are spherical shells around which the radius was the same as the radius of the Moon used in the simulation (1737 km). Then the position (r i,hit ), the velocity (v i,hit ) and time (t i,hit ) were recorded when a cloud moved from outside the range of the detector to the inside the range of the detector. The set of hits, {r i , v i, t i } hit was used to make so-called ion precipitation maps where the surface of the Moon was divided into 12 Â 6 longitude Â latitude bins and the average ion impact flux (< j >, unit: #/m 2 /s) at every bin was derived by dividing the total number of protons which hit the latitude (la) longitude (lo) bin (w la;lo ) with the surface area of the bin (dA la;lo ) and the length of the time when the ions were collected (dT) < j> la;lo ¼ w la;lo ðdA la;lo dTÞ (9) The hits were also used to derive the macroscopic plasma parameters: the proton density and the proton bulk velocity U were both investigated in a way similar to the properties of plasma near the Moon in the solar wind (see Kallio et al., 2016, for details of the algorithms). The macroscopic parameters were also interpolated along the orbit of the Moon by using plasma values which were saved at every cell.
The coordinate system was the following: The x-axis points from the Sun toward the centre of the Earth, which was located at (x, y, z) ¼ (0, 0, 0). The Earth's magnetic field was modelled as an ideal magnetic dipole. The strength of the dipolar field at the surface of the Earth at the magnetic equator was 30870 nT pointing toward the þz axis. The y-axis completed the right hand coordinate system. The orbit of the Moon, on which the virtual Moon detectors were placed, was approximated by a circle of a radius 382680 km~60 R E (R E ¼ 6371 km) around the Earth on the z ¼ 0 plane.
The inner boundary of the simulation was a sphere at r ¼ R inner ¼ 7 R E inside which the electrical resistivity was assumed to be zero. Replacing the planet with a "perfectly conducting sphere" has been commonly used in hybrid models to mimic a conducting ionosphere around a planet. In the analysed j Â B hybrid model, as in typical magnetosphere MHD models, the inner obstacle, the size of which was larger than the radius of the Earth, was used to cut away the near planet region where the magnetic field was high and a small time step should be used to follow the ions accurately. Furthermore, the ion was taken away from the simulation if it hit the inner boundary.
The size of the simulation box was [x; y; z] ¼ [À100, 50; À65, 65; À65, 65] R E . The grid size, dx, was 1.0 R E and the time step, dt, was 1 s. The initial average number of particles per cell was 30. The simulation was run until it reached a quasi-stationary solution. Particles were propagated by a leap-frog method and by using the Boris-Buneman schema (see Kallio and Janhunen, 2003, for details). The weight (w i ) of the particle cloud was 8.62⋅10 25 .
In this paper the so-called "open" magnetosphere case was investigated, when the interplanetary magnetic field, IMF, upstream of the bow shock was [B x , B y , B z ]~[0, 0, À4] nT, i.e., purely southward orientation. Generally speaking, from the kinetic simulation point of view, events when the solar wind is hot are especially interesting because of the large ion gyroradius. Therefore, in this study properties of the solar wind were chosen which mimic fast solar wind conditions when the solar wind velocity and temperatures are high and the density is low. The solar wind parameters upstream of the bow shock were the following: the solar wind density, n sw~5 cm À3 , the solar wind velocity U sw~[ 600, 0, 0] km/s, and the solar wind temperature T sw~1 0 6 K. The electron temperature, which is used to derive the thermal pressure in Eq. (6a), was also assumed to be 10 6 K.
3. Magnetosphere during a fast solar wind pass Fig. 1 gives two 3D overviews of the analysed j Â B hybrid model run. One can see in Fig. 1a how the magnetosphere causes the solar wind flow to deflect around the Earth on the z ¼ 0 plane, and how the magnetic field lines near the y ¼ 0 plane form a classical "open" magnetosphere magnetic field morphology. An open magnetopause is formed where the total magnetic field increases downwind of the shock. The magnetic field lines form magnetic tail lobes where the magnetic field is predominantly along the x-direction.
From the ion acceleration point of view the fundamental physical parameter is the electric current (j) which results in the j Â B force acceleration. The acceleration depends both on the total electric current and on the orientation of j with respect to the magnetic field. Fig. 1b is given to illustrate the electric current system and how it is related to the magnetic field structures as well as electric current regions on the y ¼ 0 and z ¼ 0 planes. Note that on the dayside the electric current makes two classical magnetospheric Chapman-Ferraro type electric current loops, one on the northern hemisphere (z > 0) and another on the southern hemisphere (z 0). On the nightside the electric current forms two toroidal currents which enclose the magnetic field lines at the magnetic tail lobes, as is anticipated to exist in the Earth's magnetotail. These tail currents form the cross tail current sheet between the magnetic tail lobes of different magnetic polarity and they close by flowing around the magnetic tail lobes. Fig. 2 presents the run in more detail. Fig. 2a shows how the magnetic field is increased at the bow shock and the magnetic field lines become draped downstream of the bow shock in the magnetosheath. The originally dipolar magnetic field becomes stretched out tailwards on the nightside and a low magnetic field region is formed near the magnetic equator (z ¼ 0) which separates the magnetic tail lobes where the magnetic field points away from the Earth (z 0 hemisphere) and toward the Earth (z > 0 hemisphere). As already mentioned, the basic morphology of the magnetic field lines shown in Fig. 2a represents a well-known open magnetosphere situation. Note that because the lunar orbit is modelled in this work as a circle in the z ¼ 0 plane, the magnetic field is weak along the lunar orbit on its path through the magnetotail, as can be seen in Fig. 2b. The density of the solar wind protons is given in Fig. 2c and d. The bow shock can be seen as enhanced density of the solar wind protons. Downstream of the bow shock the plasma density is increased in the magnetosheath. The high density region in the magnetosheath ends at the magnetopause where the low density magnetosphere is located. In addition to the magnetosheath a second, relatively weak, local density enhancement can be found at the plasma sheet near the z ¼ 0 plane. It is also important to note the following consequences of the parameters given in Fig. 2 when the ion precipitation towards the Moon is investigated. The orbit of the Moon is in the simulation on the z ¼ 0 plane and the density of the solar wind protons in this location is given in Fig. 2d. One can see that the Moon faces increased plasma density at the bow shock, enhanced density on the magnetosheath and a relatively low plasma density within the plasma sheet.
It is informative to make a closer look at the electric current because a strong electric field can result in a noticeable acceleration of ions. In addition, the values of the electric current in Fig. 1b are not clearly visible because the colour planes were behind the electric current lines. Fig. 2e and f shows clearly the increase of the electric current first at the bow shock. Then the electric current is strong at places where the tail current crosses the noon-midnight (y ¼ 0) plane (Fig. 2e). The third high electric current region is the cross tail electric current near the equator (z ¼ 0) plane. There is also a strong electric current region near the inner obstacle, which is at 7 R E . It should be noted that details of this near-Earth electric current system can be affected much by the fact that the model does not include an ionosphere and, therefore, neither magnetosphereionosphere currents. By contrast, a relatively large spherical volume around the Earth was cut away from the simulation box and replaced by a perfectly conducting sphere, as discussed above. For the same reason it can be anticipated that the closer the parameters of interest are to the inner obstacle, the more they can be affected by the simplified assumptions used in the model. However, the orbit of the Moon is so far away from the Earth (~60 R E ) that these limitations of the inner magnetosphere modelling are not expected to reduce considerably the realism of the values along the lunar orbit.
Although the purpose of this study was not to make magnetospheric research it is interesting to note how the magnetic field lines become highly aligned with the x-axis (Fig. 2a) near the x-axis at x~600⋅10 3 km where the cross tail electric current is also increased (Fig. 2e and f). These are signatures which one would expect to find at a location where a magnetic reconnection can to occur. In future, however, this has to be investigated in detail by using an even larger simulation box to minimise possible artefacts caused by the boundaries of the simulation box.

Properties of impacting solar wind protons
Figs. 1 and 2 gave some of the basic plasma and field properties in the magnetosphere. The focus of the paper is, however, to study how the magnetosphere affects the ion impact on the Moon, which is tidally locked to the Earth and which always faces the same hemisphere toward the Earth. Fig. 3 illustrates how the revolution of the Moon affects its illumination conditions. It also helps to interpret proton precipitation maps because similar to sunlight, the solar wind protons flow predominantly anti-sunward in the far magnetosphere. Therefore, areas on the Moon which are illuminated by sunlight are also areas where one would expect the highest proton precipitation flux to occur at the given portion along the orbit. The figure also illustrates the developed moving Moon centred MES (Moon-Earth-Sun) coordinate system where the x MES axis points always from the Earth to the centre of the Moon, the y MES axis is opposite the orbital motion of the Moon in the model circular orbit and the z-axis completes the right hand coordinate system. It is worth to note that in this study our own MES coordinate system was created to optimise the simulation usage, although the MES coordinate system has certain similarities with Selenographic coordinates. However, in the analysed simulation the orbital plane and the radius of the Moon as well as the tilt angle of the Earth are highly idealised and, therefore, simulation results are presented in the own MES coordinates.
At every five positions of the Moon depicted in Fig. 3 small white background inserts show the impact positions of protons on the lunar surface, i.e., ion impact maps, in the latitude-longitude coordinates. In the inserts the latitude is the angle between the position of the hit (x, y, z) in the MES coordinates where the latitude is þ90 at the lunar North pole (þz MES ), 0 at lunar equator (z MES ¼ 0) and À90 at the lunar South pole (-z MES ). Longitude gives the angle between the position of the hit (x, y) on the MSE xy-plane and the -x MES axis so that the longitude is [À180 , 0 ] on the y MES > 0 side and [0 , 180 ] on the y MES < 0. This means that the point (latitude, longitude) ¼ (0 , 0 ) always points directly toward the centre of the Earth and it is at the centre of the nearside hemisphere. The point (0 , AE180 ) is, instead, at the centre of the farside hemisphere and, consequently, is never visible from the Earth. In addition, the longitude range [À90 , 90 ] are on the nearside hemisphere while the longitude ranges [À180 , À90 ] and [90 , 180 ] are on the farside hemisphere.
In the five examples of the ion impact maps given in Fig. 3 every point corresponds to a hit of an ion cloud on the lunar surface and, consequently, the more dots per area element there are on the map, the higher the proton precipitation flux is at the given position on the lunar orbit. The latitude-longitude maps show quantitatively that proton impacts on the lunar surface occur often on regions which are also illuminated by the sunlight. This is not, however, strictly true because the magnetosphere deflects the bulk flow of the solar wind away from the Earth. Moreover, the solar wind can be highly heated downstream of the bow shock and, consequently, the velocity of protons can differ noticeably from the direction of the solar wind flow.
As mentioned before, the properties of the impacting protons were investigated by using two different datasets: macroscopic parameters on the simulation cells and recorded proton impacts on the lunar surface. Interpolated plasma and field values along the lunar orbit on the Earth's nightside are given in Fig. 4. One can see that when the Moon enters from the solar wind into the magnetosheath at phase angle ϕ~10 the proton density (n), the proton thermal pressure p (¼ n k B T, where k B is the Boltzmann constant) and the total magnetic field (B) are increased while the proton bulk velocity (U) is decreased. Statistical fluctuations in the proton temperature (T) are too large to identify bow shock crossing from that parameter. Parameters change in the opposite way when the Moon is moving out of the magnetosheath into the solar wind at ϕ~170 . One can also see that the Moon faces the smallest plasma densities, bulk velocities and the magnetic fields along the Moon's orbit near the centre of the tail around ϕ~90 . The temperature and thermal pressure of protons are instead, highest near the centre of the tail. Note also that the magnetic field is very small within the tail in Fig. 4 because the lunar orbit goes exactly through the cross tail current sheet and the plasma sheet.
Although the purely z-directed IMF that is used does not represent a nominal or statistical solar wind condition for the calculation of the magnetosphere and the magnetotail, it is however informative to put the derived plasma and field parameters along the lunar orbit in a wider perspective. To do that, Fig. 4 presents also plasma parameters based on a 3D MHD simulation where a full Earth year is simulated by using observed parameters as modelled input parameters (see Kallio and Facsk o, 2015, for details of the used MHD runs). During one year the Earth faces many different IMF situations with non-zero x and y components and different magnetic dipole tilt angles. Therefore, the MHD statistics include various kinds of full 3D magnetosphere configuration cases, tilted and rotated cross tail currents sheets and plasma sheets, positions of the bow shock and magnetosphere as well as closed and open magnetosphere cases. For example, it can be anticipated that the analysed open magnetosphere case, which includes clear reconnection, can result in higher plasma acceleration than on the more typical IMF cases. This might be one reason why the plasma velocities in the analysed case exceeds the average statistical MHD velocities (Fig. 4b). One should also recall that the plasma sheet is a narrow layer (cf. Fig. 2c) and in the analysed IMF case it is on the z ¼ 0 plane where also the lunar orbit is modelled to be. In the reality, however, the plasma sheet can be twisted or tilted and the lunar orbit does not remain always within the plasma sheet. Nevertheless, the similarities between the analysed j Â B hybrid model run and the statistical MHD values, in terms of absolute magnitude and general trends of how parameters change along the lunar orbit, suggests that the conclusions shown later in the paper can be considered to represent also features which occur generally in the magnetosphere and not only during the analysed pure southward fast solar wind case.
The second data set, collection of impacting protons into the particle detectors, provides further details about the properties of the impacting protons. The macroscopic parameters derived from the particle detectors are shown in Fig. 5a -5c. It should be noted that the particle detector data has more statistical fluctuations than the average macroscopic data which has been shown Fig. 4. This occurs because the macroscopic parameter derived from particle detectors includes only limited amount of particles, as illustrated in five inserts in Fig. 3. The macroscopic parameters, instead, are derived by using "particle clouds" which have the same size as the size of the simulation cell and, therefore, at every time step plasma parameters at a given cell can be contributed by particles within nine cells. Moreover, the used time averaging decreases statistical fluctuations even more and provides relatively "smooth" macroscopic plasma parameters. For example, in Fig. 5a the particle density at waxing gibbous looks to be higher than that at waning gibbous although in the average data (Fig. 4a) no such asymmetry can be seen. The origin of this asymmetry is most likely statistical fluctuations. However, particle data are valuable comprehensive data because they provide information which can not obtained from macroscopic parameters. The total particle flux (Fig. 5d) and the total energy flux (Fig. 5e) on the lunar surface have a global minimum near the centre of the tail.
However, both the total energy density and the average proton kinetic energy (Fig. 5f) density have their global maximum values around the global minimum value. In addition, local minimum values in the average proton kinetic energy can be seen downstream of the bow shock at around ϕ~20 and at ϕ~160 . Note also how the magnetosphere deflects the solar wind flow downstream the bow shock, which can be seen as a positive sign of U y when the Moon moves through the bow shock into the magnetosheath at ϕ~20 -40 and a negative sign of U y when the moves from the magnetosheath back to the solar wind at ϕ~140 -160 (Fig. 5c). Deeper in the tail at ϕ~60 -120 the sign of U y indicates that the flow is converging into the centre of the tail implying that ions are filling the plasmasheet.  The normalised total cumulative proton flux < j > along a lunar orbit in the model is shown in Fig. 6a. It shows that the highest cumulative flux is on the farside opposite the Earth (latitude ¼ 0 , longitude ¼ AE180 ). The minimum cumulative flux, which is about 70% of the maximum flux, is located at the nearside closest the Earth (latitude ¼ 0 , longitude ¼ 0 ). This longitudinal asymmetry, or the farside-nearside (FN) hemisphere asymmetry, would not exist without the magnetosphere and shadowing of the Earth. In a hypothetical situation, where the Earth would not have an intrinsic magnetic field (or its size would be very small) disturbing the flow of the solar wind plasma, the Moon would be all the time in the undisturbed solar wind and the cumulative particle flux would not have longitudinal variations but only latitudinal variations, as shown in the simulation shown in Fig. 6b.

Cumulative proton impact flux
The farside-nearside asymmetry in the impact flux results from the fact that the particle flux to the farside during the time when the Moon is in the undisturbed solar wind (Fig. 6c) is larger than the particle flux to the nearside during the time when the Moon is in the disturbed solar wind (Fig. 6d) in the magnetosheath and in the magnetotail. The Earth's magnetosphere affects therefore on the solar wind protons similarly as a diverging optical lens would have an effect on the sunlight. Because of the tidal locking of the Moon's rotation, the same region on the lunar surface is always facing toward the attenuated or enhanced solar wind flow at the nightside of the Earth resulting in the farside-nearside cumulative proton precipitation asymmetry.
In addition to the proton impact flux, another important parameter when the space weather effects are considered is the proton energy flux and its energy dependence. For example, lunar magnetic anomalies reflect the lowest energy protons before they hit the surface (Bamford et al., 2012;Lue et al., 2014) and, therefore, their possibility to result in space weathering effects near the surface. Fig. 7 gives the energy distribution function, f(E) [# bin /(dE bin <E bin > m 2 s)], along the lunar orbit derived from the virtual Moon Detectors No 0-18, i.e., detectors on the Earth's nightside from ϕ ¼ 0 to ϕ ¼ 180 (c.f. Fig. 3). The values are derived by dividing the energy range 0 eV-5 keV into constant energy bins dE bin ¼ 100 eV, calculating the number of precipitating protons (# bin ) the energy of which is at energy bins, and dividing the result with the centre energy of the bin (<E bin >). The figure shows that while the Moon faces in the solar wind a relatively beam-like proton flow, the energy distribution function becomes wider, i.e., hotter, downstream of the bow shock. At the same time, the bulk velocity decreases. Moreover, there are regions where the precipitating protons have relatively low energies as well as high energies in comparison to the undisturbed solar wind.
An interesting question is how protons with different energies are precipitating spatially on the surface and, especially, do the low and high energy protons have spatial cumulative anisotropies? This was investigated in Fig. 8 which shows the cumulative ion impact flux on one lunar orbit at three different energy ranges: (i) low energy protons at energies  less than 1500 eV, (ii) nominal energy protons in the energy range 1500 eV-2500 eV and (iii) high energy protons at energies above 2500 eV. These energy ranges were chosen because, according to Fig. 7, the number of protons in the energy ranges E < 1500 eV and E > 2500 eV relative to the energy range 1500 eV < E < 2500 eV is higher in the Earth's nightside than in the undisturbed Earth's dayside. Therefore, it can be anticipated that also spatial non-homogenous cumulative proton impact patterns can be formed because the tidal locking of the Moon always ensures that the same face is pointing toward the Earth at the same position on the orbit. Fig. 8 shows that the cumulative flux of protons with different energies is not distributed isotopically on the surface. As one can see by comparing Fig. 8b and Detector panels No 5 and No 13 in Fig. 5, these two maximum flux surface regions are formed by the proton precipitation on the Waxing Gibbous phase (ϕ~45 ) and the Waning Gibbous (ϕ~135 ) phase between the bow shock and the centre of the tail. The cumulative proton impact pattern flux of the nominal energy protons (Fig. 8b) resembles the total ion impact map shown in Fig. 6 because most of the precipitating protons are within the nominal energy range. The highest energy protons, in turn, hit predominantly around the centre of the nearside (Fig. 8c). This occurs because the highest average energy of the protons is found around the centre of the tail (see Fig. 5f) when the centre of the nearside faces against the outflowing magnetospheric plasma. Fig. 8a represents the low energy ion impact regions. This is shown for completeness although these low energy ions do not form such a clear ion impact patterns on the surface than the medium and the high energy ions. However, the low particle flux region (the region between contour values 1 and 3 in Fig. 8a) of the low energy ions looks to be located in large areas around the centre of the nearside hemisphere.

Discussion
We have investigated the properties of the solar wind protons impacting on the lunar surface by the 3D j Â B hybrid model. The general scientific goal was to understand the global effects of space weathering on the surface of the Moon over the full Earth orbit, as well as the properties of the protons and electromagnetic fields near the Moon and above its surface. Especially, the primary goal was to analyse the precipitation self-consistently of protons by taking full account of kinetic effects and the full 3D proton velocity distribution function.
The simulation results in maps of precipitating protons, which can be of value as an input in dedicated specific space weathering models. The j Â B hybrid model was built in the simulation platform that has been used earlier to study the interaction of solar wind plasma with the lunar surface in 3D (Kallio, 2005), the Debye layer near the lunar surface (Dyadechkin et al., 2015), the properties of plasma near lunar magnetic anomalies (Jarvinen et al., 2014) and the motion of charged dust particles near the lunar surface and around the Moon (Dyadechkin et al., 2015). Therefore, the compiled proton impact maps can be directly used in such kind of studies and to obtain more detailed information of the surface charging, the electric field near the surface, the role of particle sputtering in the formation of the exosphere, the properties of Energetic Neutral Atoms (ENAs) from the surface, the motion of electrically-charged lunar dust particles, and the relation between lunar magnetic anomalies, their effect on the local plasma, and swirls on the surface. The model contains a full 3D ion velocity distribution function and, therefore, it is foreseen that it can deepen our understanding of how flowing plasma enters lunar craters (Zimmerman et al., 2011) and how the lunar terminator region is electrically charged by the particle precipitation (Farrell et al., 2007).
From the permanently shadowed craters near the lunar poles, and the need to understand the cumulative effects of proton precipitation and impact on the lunar surface, also arises the question of what are these long-term cumulative effects when the time scale is hundreds of millions years? The proton impact map analysed in this study indicates that the tidal locking of the Moon combined with the magnetosphere results in the cumulative farside-nearside (FN) hemispherical proton impact asymmetry. During the history of the Moon the strength of the FNasymmetry has, however, been different to what it is at the moment. For example, the distance of the Moon to the Earth has been smaller in the past. Moreover, the Earth's intrinsic magnetic field, as well as the properties of the solar wind, has changed during the history of the Moon. One can, therefore, envisage that the effects of the magnetosphere on proton precipitation could have been stronger in the past than these days because the Moon has spent less time in the solar wind.
Furthermore, possible changes in the angle of the Moon with respect to the ecliptic plane and lunar orbital plane are also expected to play a role in the cumulative precipitation effects of protons. It is also worth to note that the time a dust grain spends on the lunar surface is a function of its size. For sub-micron grains that time is around 10,000 years and becomes longer with larger size (see Szalay and Hor anyi, 2016, for the most recent residence time estimations). The residence time of a large centimeter-size regolith grain on the lunar surface can be in the order of 100 Ma to 1-2 Ga (see e.g. Lorenzetti et al., 2005). Thus, visible and cumulative effects on the surface cannot be older than these residence times. Solar wind includes also ions where the mass and electric charge state is larger than for protons (e.g. Wurz et al., 2001;Wimmer-Schweingruber et al., 2006). Therefore, inclusion of He þþ and possibly heavy solar wind ions, especially O 6þ ions, into the model would also be a valuable extension because they can enhance the total sputter yield (e.g. Killen et al., 2012). A detailed simulation study of the cosmogenic cumulative effects of proton precipitation on the Moon would, therefore, require models for the Earth's intrinsic magnetic field, properties of the Sun in time and the astronomical properties of the Moon (the distance and the tilt angle).
An interesting question is also what is the role of ions that are escaping from the Earth's ionosphere and atmosphere and which hit the lunar surface? The motion of ionospheric particles in the magnetosphere can be made by a test particle simulation where the ion motion is given by the electric and magnetic fields. In a more comprehensive approach, the model would be self-consistent such that the ionospheric particles could affect the magnetospheric configuration. Technically, similarly as in the hybrid model, it is simple and computationally economical to add different ions species into the j Â B hybrid model, such as planetary ions or solar wind alpha particles. Also the inclusion of multiple charged heavy ions, such as O 6þ from the solar wind, provides interesting possibilities because the charge exchange generates soft X-rays on the surface (e.g. Kallio et al., 2008) and which can be used for plasma diagnostic near solar system objects, such as near comets, the Earth, Jupiter and Saturn.
The j Â B hybrid model yielded a magnetosphere configuration that reproduced the major well-known general magnetic field configurations and properties of plasma. On the other hand, there are many reasons why the analysed space weathering case is just the first step toward a more advanced j Â B hybrid model. For example, the role of the magnetosphere-ionosphere coupling and the role of the radius of the inner boundary have to be investigated, such as varying the radius of the obstacle and the electrical conductivity inside the inner obstacle. The Earth's intrinsic magnetic field was also assumed to have a zero-tilt angle, which is a useful test case to investigate asymmetries or possible artefacts in the simulation. However, later a dipole tilt angle corresponding to the analysed time instance has to be included to make one-to-one comparison between the simulated results and observations.
Although the analysed pure southward IMF case has many characteristic features that are classical in magnetospheric physics, the analysed fast, very hot solar wind like event occurs relatively seldom. Moreover, although the pass of a fast solar wind is highly interesting in terms of basic space plasma physics and aurora physics, it does not represent a typical solar wind situation. For example, the bow shock in the analysed fast solar wind case is further from the Earth as it typically is (see e.g. Verigin et al., 2001) and runs with more typical solar wind conditions should be made in the future. Although the derived cumulative proton impact maps are capturing the basic features that occur also during nominal solar wind conditions, more detailed study will be needed for regular solar wind conditions, e.g. such as simulating the magnetosphere over several months and even a year, as done recently based on a 3D MHD simulation (Kallio and Facsko, 2015).
It should also be noted the orbit of the Moon was modelled in this study as a circle in the ecliptic plane for simplicity. However, when the long-term cumulative effects are considered one should take into account the evolution of the lunar orbit. Lunar orbital cycles affect how much time the Moon stays on average at different regions in the Earth's magnetosphere and in the magnetotail and, for example, in the plasma sheet during its approximately 18.6 year cycle (Hapgood, 2007). In the future the derived ion impact maps should also be compared with plasma and field observations measured near the Moon, such as plasma and ENA observations made by Kaguya (e.g. Harada et al., 2010;Nishino et al., 2017), Chandrayaan-1 (e.g. Barabash et al., 2009;Futaana et al., 2008Futaana et al., , 2010Lue et al., 2014) and ARTEMIS (e.g. Poppe et al., 2012Poppe et al., , 2017 missions.

Conclusions
The impact of solar wind protons on the lunar surface was investigated using a 3D hybrid model where protons are accelerated by the j Â B force. The study showed that in the analysed fast solar wind case the total cumulative proton impact rate is higher on the farside of the Moon than on the nearside. Spatial asymmetry was also found for the impact energy: protons with the smallest energies hit the lunar surface predominantly near the farside-nearside terminator while protons with the highest energies hit most frequently near the centre of the nearside. Space weathering effects depend on the flux and the energy of the protons and, because of the tidal locking of the Moon and the existence of the Earth's magnetosphere, they are spatially non-isotropic over the hemisphere, which deflects the flow of the solar wind around it.