© F. Lavorenti et al. 2022

1. Introduction

Mercury is one of the least explored planets of the Solar System. Over decades of space exploration, only two missions have been devoted to the innermost planet of the Solar System.

The NASA Mariner10 mission in the 1970s provided a snapshot of the Hermean environment with its three flybys (Russell et al. 1988). These observations showed the presence of a planetary magnetic field (Ness et al. 1974) and of a structured plasma environment (Ogilvie et al. 1977). In particular, Mariner10 was able to perform measurements of the core of the electron distribution function (in the range 13.4–687 eV) during most of its orbit inside the Hermean magnetosphere (Christon 1987). However, Mariner10 was not able to observe ions due to a technical failure (Ogilvie et al. 1974).

Three decades later, the NASA MESSENGER mission deeply extended our knowledge of the Hermean environment during its four years of orbital observations (Solomon & Anderson 2018). This wealth of in situ observations enabled the first systematic studies of the interaction between the solar-wind plasma and Mercury’s magnetized environment (Raines et al. 2015 and references therein). These studies evidence the highly dynamical character of Mercury’s plasma environment and shed light on the tight coupling between the solar-wind, magnetosphere, exosphere and surface of Mercury. The mission addressed several plasma processes occurring at the global planetary scale (of the order of 2400 km) and down to the ion kinetic scale (of the order of 100 km) (Boardsen et al. 2012; Raines et al. 2014; Gershman et al. 2014, 2015; Schmid et al. 2021). However, given the instrumental constraints of the mission, MESSENGER did not address the plasma processes ongoing at the electron scale. The instrumental suite onboard MESSENGER provided observations of electrons above ∼10 keV, thus missing the core of the distribution function, and it could not provide plasma wave observations below the typical ion frequency range (of the order of a few Hz).

The ESA/JAXA BepiColombo mission has been designed to continue building upon the legacy left by MESSENGER (Benkhoff et al. 2021). BepiColombo is now cruising to Mercury (arrival expected in 2025) and performed its first Mercury flyby on 1 October 2021. The novelty of this mission resides in (i) its two-spacecraft nature (BepiColombo is composed of one Mercury Planetary Orbiter, MPO, and of one magnetospheric orbiter, nicknamed Mio) and in (ii) its advanced instrumental suite for environmental studies enabling observations down to the electron scale (Milillo et al. 2020). BepiColombo is the first mission able to provide a simultaneous multi-point picture of the Hermean environment from the global planetary scale down to the electron scale (of the order of 1 km). In order to both analyze and optimally plan such novel observations, global numerical models of the Hermean environment with resolutions down to the electron scale are needed.

To date, global numerical models of Mercury’s plasma environment have been mostly limited to magnetohydrodynamic (MHD; Kabin et al. 2000; Ip & Kopp 2002; Yagi et al. 2010; Pantellini et al. 2015; Jia et al. 2015, 2019), multi-fluid (Kidder et al. 2008; Benna et al. 2010; Dong et al. 2019), and hybrid (meaning kinetic ions and fluid massless electrons) (Kallio & Janhunen 2003; Trávníček et al. 2007, 2009, 2010; Richer et al. 2012; Müller et al. 2012; Herčík et al. 2013; Herčík & Trávníček 2016; Fatemi et al. 2018, 2020; Exner et al. 2018, 2020) models. These models appropriately reproduce the plasma dynamics down to the ion kinetic scale and have been extensively and successfully used in support of the Mariner10 and MESSENGER missions. However, to interpret the forthcoming new observations of the BepiColombo mission, more refined models able to self-consistently include electron kinetic physics are required.

Recently, a first attempt to locally embed electron kinetic physics in a global MHD simulation was carried out by Chen et al. (2019). Their model can be used to study the role of electrons in a precise subset of the global magnetosphere (for instance, the magnetotail); on the other hand, it cannot reproduce dynamical processes encompassing the global magnetosphere system such as, for example, the global electron circulation around the planet (Walsh et al. 2013). In order to overcome this limitation and to self-consistently include both ion and electron physics, in this work, we study the interaction between the solar wind and the Hermean plasma environment using a global fully-kinetic model. A similar approach was presented by Lapenta et al. (2022) using the results of a global hybrid simulation to initialize a fully-kinetic one. In that work, the authors found that the feedback of kinetic electrons does not affect the large-scale structure of the magnetosphere. Nonetheless, on the small scale, they found (i) thinner magnetospheric boundaries with higher currents and (ii) more efficient energization processes in the magnetosphere. The results of Lapenta et al. (2022) point out the importance of a kinetic model for electrons in global simulations of planetary magnetospheres, a goal that is today attainable thanks to the increasing computational power of current High Performance Computing (HPC) facilities.

Although only a few in situ electron observations are available at Mercury, some hints on their global dynamics have been discussed in past works using Mariner10 (Christon 1987) and MESSENGER data (Ho et al. 2012, 2016; Baker et al. 2016; Dewey et al. 2018). One of the most significant outcomes of these observations is the presence of a quasi-stable, high-energy (up to tens of keV) electron population inside the Hermean magnetosphere, mainly observed in the night post-midnight sector (local time 0–6 h). According to those observations, the Hermean environment is significantly populated by electrons with energies up to tens of keV. In this work, we investigate the origin of these energetic electrons and their loss mechanism inside the Hermean magnetosphere by means of global, fully kinetic simulations.

The paper is organized as follows. In Sect. 2, we present our fully kinetic model, and in Sect. 3 we analyze the large-scale structure obtained from our simulations and validate our model against nominal bow shock and magnetopause shapes and positions at Mercury. In Sect. 4, we focus on our simulations results regarding electron energization and circulation in the Hermean magnetosphere, and in Sect. 5 we build and discuss synthetic electron energy spectra obtained from our simulations along Mariner10 and BepiColombo trajectories during their first Mercury flybys. Finally, in Sect. 6, we discuss the results and conclude the paper.

2. The model: Fully kinetic global simulations

The simulations are done using the semi-implicit, fully kinetic particle-in-cell (PIC) code iPIC3D (Markidis & Lapenta 2010). It solves the Vlasov-Maxwell system of equations for both ions and electrons by discretizing the distribution function of both species using macro-particles. We present two different simulations, namely RunN and RunS, characterized by a purely northward or southward interplanetary magnetic field (IMF), respectively. In the following, we use the Mercury-centered Solar Orbital (MSO) reference frame, where the x-axis points from the planet center to the sun, the z-axis is anti-parallel to Mercury’s magnetic dipole, and the y-axis points from dawn to dusk.

The simulations use a three-dimensional cartesian geometry, a simulation domain (L_x, L_y, L_z) = (−9 : +6,±6,±6) R (here, R is the planet radius) in MSO coordinates, divided into (N_x, N_y, N_z) = (960,768,768) grid cells. The simulations are initialized with 64 macro-particles per cell per species. The time step is dt = 1.4 ms, much smaller than the electron gyro-period (τ_ce ≡ 2π/ω_ce ≡ 2πm_ec/eB_SW = 31.5 ms). The grid spacing is dx = dy = dz = 0.015 R = 1.5 ρ_e, where ${\rho }_{\mathrm{e}}\equiv c\sqrt{T_{\mathrm{e},\mathrm{SW}}{m}_{\mathrm{e}}}/e{B}_{\mathrm{SW}}$ is the electron gyro-radius. One important feature of implicit PIC codes is that, even if the grid spacing is of the order of the electron gyro-radius ρ_e ≲ dx < ρ_i, the code correctly captures the sub-grid electron kinetic physics as long as the condition 0.1 < v_thedt/dx < 1 is fulfilled (Brackbill & Forslund 1982). In our case, v_thedt/dx = 0.18. However, since sub-grid oscillations on the electron scale tend to accumulate energy at the grid scale, a smoothing technique is required to avoid numerical instabilities at the grid scale (see Tóth et al. 2017 for details on the filter used; in our case, the filter is applied two times per cycle with a factor of α = 0.5).

The simulation setup includes (i) the solar-wind plasma (uniform density, magnetic field, and temperature), (ii) the spherical planet centered at x, y, z = (0, 0, 0) with radius R, and (iii) the dipolar magnetic field centered in x, y, z = (0, 0, 0.2) R with intensity 200 nT/R³.

Given that we simulate an infinite system using a finite numerical box, particular care must be taken when imposing the external boundary conditions. Because the planet interaction region is sufficiently far (on the order of 5 planet radii) from the external boundaries of the box, we populate the last cells of the boundaries with solar-wind plasma (an exception is made for the antisunward boundary from which plasma is only allowed to exit the box). To smooth out the electromagnetic fields fluctuations close to the boundaries, we employ a linear simple absorption layer (LSAL) method (Berendeev et al. 2018) on the ten outermost cells of the simulation domain. As discussed in Berendeev et al. (2018), this method avoids wave reflection and ensures numerical stability at the boundaries. With this method, the magnetic and electric field values in the outermost boundary cells are equal to their corresponding solar-wind values, B_SW and E_SW. Since in the solar wind the plasma frozen-in condition holds, the electric field is given by E_SW = −V_SW × B_SW. A divergence cleaning routine for the magnetic field ensures that div(B) = 0 in the absorbing layer.

Particular care must also be taken in handling the plasma interacting directly with the planet. This is especially true for models, such as fully kinetic ones, that do not impose quasi-neutrality. In order to avoid the generation of spurious net charge density on the planet surface, the macro-particles falling into the planet are removed from the simulation using a charge-balanced method. This method ensures that (i) the same amount of positive and negative charged macro-particles is removed from the simulation box at the planet surface boundary and that (ii) the excess of electron flux at the planetary surface is expelled back radially. This inner boundary condition mimics the sub-grid interaction between precipitating plasma and the planet surface ongoing on length scales of few Debye lengths.

In both simulations, the solar-wind plasma is initialized with a density of n_SW = 30 cm⁻³, velocity of V_SW = ( − 400,0,0) km s⁻¹, magnetic field of B_SW = (0,0,±20) nT, and temperature of T_i, SW = T_e, SW = 21.5 eV. The ram pressure is $P_{\text{ram}}\equiv n_{\text{SW}}{V}_{\text{SW}}^2=8.2$ nPa, and the plasma beta is $\beta \equiv 8\pi n_{\text{SW}}{T}_{\text{SW}}/B_{\text{SW}}^2=1.3$. The only parameter that is different between the two runs (RunN and RunS) is the direction of the magnetic field B_SW. These parameters are representative of those expected in the proximity of Mercury at aphelion (James et al. 2017; Sarantos et al. 2007) and similar to those used in Aizawa et al. (2021).

Once the interaction between the solar-wind flow and the planetary magnetic field reached a quasi steady-state, we stopped the simulation. Such a timescale corresponds approximately to the solar-wind ballistic time taken to cross the simulation box, T ≈ L_x/V_x, SW ≈ 10 s (corresponding to approximately 7000 numerical cycles).

In order to keep computational resources manageable, we chose to reduce (i) the ion-to-electron mass ratio m_i/m_e = 100 (from its real value of 1836), (ii) the electron plasma-to-cyclotron frequency ratio ω_pe/ω_ce = 17.8 (from its real value of 83), and (iii) the planet radius R = 230 km (from its real value of 2440 km). With these values, we obtain a hierarchy of scale lengths R = 10 ρ_i = 100 ρ_e that – although compressed – maintains a sufficiently large separation between planetary, ion, and electron kinetic scales. The first and second rescalings represent a standard procedure in plasma fully kinetic simulations (see, e.g., Bret & Dieckmann 2010). The third rescaling is analogous to what was proposed and validated in Tóth et al. (2017) for the case of the Earth and also discussed extensively in Markidis et al. (2021). Moreover, this planet rescaling approach was already adopted in past works using global hybrid models in support of MESSENGER observations (Trávníček et al. 2007, 2009, 2010) and recently in a global fully kinetic model (Lapenta et al. 2022).

In the next section, we validate our scaled-down model against the nominal shapes and positions of the magnetosphere boundaries. A broader discussion on the impact of these rescalings on our results is given in Sect. 6.

3. Overview of the large-scale structure: Magnetosphere boundary identification

3.1. Validation of the global, fully kinetic model

We start by validating our model by showing that the large-scale structure of the system typically obtained by fluid and hybrid models is correctly reproduced, as shown in Fig. 1. Indeed, both simulations show the formation of a bow shock standing in front of the planet, the magnetosheath with increased density and magnetic field amplitude, the magnetopause with its strong current, J_y, and the magnetosphere cavity with very low density and high magnetic field. Such a large-scale configuration is shown in Fig. 1 for both simulations using the ion density n_i [cm⁻³] (left panels), the magnetic field amplitude |B| [nT] (central panels), and the plasma current J_y [nA m⁻²] (right panels). In Fig. 2, we show a cut along the subsolar line, Y_MSO = Z_MSO = 0.

Fig. 1. Overview of large-scale structure and magnetosphere boundaries in our simulations. Panels a–f show the results of RunN. Panels a–c show dipolar YMSO = 0 cuts of the ion density ni (a), magnetic field amplitude |B| (b), and plasma current density along the y-axis Jy (c). Panels d–f show equatorial ZMSO = 0 cuts of the same quantities for RunN. Panels g–l show the results of RunS in the same format as those of RunN. All plots represent a time of t = 15.9 R/Vx. The average bow shock (dashed black line) and magnetopause (solid black line) profiles found by Winslow et al. (2013) using MESSENGER observations are added. The white region around the planet corresponds to grid cells with zero macro-particles.

The observed bow shock and magnetopause are compared to their nominal shape and position obtained from MESSENGER data (Winslow et al. 2013). The authors modeled the bow shock using a paraboloid model (with parameters X₀ = 0.5 R, e = 1.04, p = 2.75 R; Slavin et al. 2009) and the magnetopause using the Shue et al. (1997) model (with parameters R_ss = 1.45 R, α = 0.5). The resulting bow shock and magnetopause profiles are shown in Fig. 1 using dashed and solid black lines, respectively. Overall, the bow shock and magnetopause found in our simulations are in good agreement with their nominal shapes. As expected, a better agreement is observed in the equatorial plane (see Fig. 1), whereas some discrepancy is found in the dipolar plane due to the broken assumption of cylindrical symmetry assumed in these relatively simple analytical models (the magnetic field along the z-axis breaks the symmetry).

We compute the bow shock and magnetopause stand-off distances as the maximum of the current density, |J(x)|,  and the point where the magnetic pressure equals the solar-wind ram pressure, respectively. The stand-off positions of the bow shock and magnetopause obtained from our simulations are in agreement with the nominal values reported by Winslow et al. (2013), as shown in Fig. 2. The relative discrepancy for the bow shock and the magnetopause is on the order of 1% and 10%, respectively. Interestingly, the discrepancy for the magnetopause is larger in both runs when compared to the one for the bow shock. Nonetheless, similar discrepancies of the order of 10% have been found by Aizawa et al. (2021) using global hybrid and fluid models. We conclude that our fully-kinetic model satisfactorily reproduces the large-scale structures of the Hermean magnetosphere. This validates our model at large scales and paves the way for further analysis of the small-scale electron physics in the next sections.

Fig. 2. Cut of results of our two simulations RunN (left panels) and RunS (right panels) along the subsolar line (x-axis). Panels a–b: proton and electron densities. Panels c–d: magnetic field components and amplitude. Panels e–f: current density along YMSO of protons, electrons, and total. Panels g–h: magnetic Pmag ≡ B2/2μ0, ram $P_{\text{ram}}\equiv m_{\text{i}}{n}_{\text{i}}{V}_{x,\text{i}}^2+m_{\text{e}}{n}_{\text{e}}{V}_{x,\text{e}}^2$, and total pressures. These cuts were performed at the same time as Fig. 1 (t = 15.9 R/Vx). We highlight the bow shock and magnetopause stand-off distances obtained from our model (solid red vertical lines) and the ones obtained using the average model of Winslow et al. (2013) constrained by MESSENGER observations (dashed red vertical lines).

3.2. Impact of electron physics at large scales: Boosting the charge current at the boundaries

Past works using hybrid or fluid models found a total current density at the bow shock and magnetopause on the order of 100 nA m⁻² (Janhunen & Kallio 2004; Benna et al. 2010; Exner et al. 2020; Aizawa et al. 2021). From our fully kinetic simulations, we instead find significantly higher values of the current density on the order of 1000 nA m⁻², as shown in Figs. 2e,f. Such higher values are due to the presence of a dominant electron current at the magnetic field boundaries, as shown in Figs. 2e,f.

Past global fluid (Benna et al. 2010; Aizawa et al. 2021) and hybrid (Janhunen & Kallio 2004; Exner et al. 2020; Aizawa et al. 2021) simulations treated electrons as a massless neutralizing fluid. Under such assumptions, the electron current was computed from Ampère’s equation using the curl of the magnetic field and the proton current. Differently from those works, in our model the electron current is computed self-consistently from their distribution function evolving under the Vlasov-Maxwell system of equations. Since the magnetic field amplitude on both sides of the magnetospheric layers in our work is consistent with past simulations, we conclude that the width of these layers (bow shock and magnetopause) tends to be overestimated in global fluid and hybrid models with massless electrons. This result highlights the key role of electrons as current carriers in the magnetic boundaries of the Hermean magnetosphere, and the importance of an appropriate modeling of this species even in global simulations to satisfactorily address large-scale planetary boundaries. In other words, even if the location of the magnetospheric boundaries is well-reproduced with models using fluid massless electrons, a more physically correct treatment of the electrons is required to draw conclusions on the physics at the magnetopause at and below the ion scale.

3.3. IMF configuration and magnetic reconnection: Feedback from small to large scales

By changing the direction of the IMF, a strong reconfiguration of the global structure of the magnetosphere is observed, as shown in Figs. 1 and 2 where we show the results of the two different IMF configurations. Such global reconfiguration is primarily driven by magnetic reconnection at the nose (around X_MSO ≈ 1.5 R) and in the tail (around X_MSO ≈ −2 R) of the magnetosphere. It is much more active in RunS (run with southward IMF) as compared to RunN (run with northward IMF). These regions are highlighted in Figs. 1i and 2f for RunS by the plot of the current J_y. The position of the neutral line in the tail in RunS is in good agreement with MESSENGER observations (Poh et al. 2017).

Magnetic reconnection has a strong impact on the global dynamics of the magnetosphere since it rapidly changes the large-scale magnetic field topology, thus allowing for an efficient plasma injection in the magnetosphere and fast conversion of magnetic to particle energy (Vasyliunas 1975; Yamada et al. 2010; Treumann & Baumjohann 2013). Magnetic reconnection is intrinsically a multi-scale phenomenon strongly coupling the ion and sub-ion kinetic scale dynamics. Therefore, a precise description of the physics at the electron scale is important to correctly capture the features observed in laboratory and space plasmas. This is usually achieved by means of fluid or hybrid models with inertial electron closures (Wang et al. 2015; Ng et al. 2017; Finelli et al. 2021; Fadanelli et al. 2021; Jain et al. 2022) or by means of fully-kinetic models (Pritchett 2001a,b; Divin et al. 2007; Haggerty et al. 2015) as in this work.

Magnetic reconnection drives electron acceleration and heating, thus feeding the global magnetosphere circulation pattern with suprathermal electrons. We focus on such a strong impact of magnetic reconnection on the electron dynamics in the next section.

4. Electron heating and circulation in the magnetosphere: Global feedback of magnetic reconnection

As outlined in Sect. 1, the role of the electrons in the interaction between the solar wind and the Hermean magnetosphere is still a matter of debate. In this section, we show that electrons are strongly energized and partially trapped in the Hermean magnetosphere under southward IMF conditions. We focus on the results of RunS since electron heating and injection in the magnetosphere are much more efficient for southward IMF conditions, as shown in Fig. 3 by comparing the two runs. In the case of northward IMF – although electrons are heated efficiently at the magnetopause boundary, their flux across the magnetopause is negligible; in turn, virtually no electrons interact directly with the planetary dipole. A more detailed comparison between the two runs is presented in Sect. 5.

Fig. 3. Time evolution of electron temperature in nose and tail reconnection regions. Dipolar YMSO = 0 cut of simulations RunN (top panels) and RunS (bottom panels) at different times (time advancing from left to right). We show the color map of the electron temperature Te = (2T⊥, e + T∥, e)/3 in units of [eV] in logarithmic scale. On top of that, we superpose the in-plane magnetic field lines (gray solid lines). The white region around the planet corresponds to grid cells with a number of macro-particles too small to reliably compute the value of Te.

The electron dynamics in RunS is mainly driven by (i) the occurrence of magnetic reconnection causing strong acceleration and heating and (ii) the planetary dipole magnetic field causing the trapping and drifting of electrons inside the magnetosphere. In the following, we discuss these two mechanisms separately.

First, magnetic reconnection affects plasma parcels flowing through the reconnection region by (i) accelerating the plasma (increase of the plasma bulk flow velocity) up to the Alfvén speed V_A, in computed from the inflow plasma parameters¹ and (ii) heating the plasma (increase of the thermal speed of the particles) by an amount proportional to the inflow magnetic energy $\mathrm{\Delta }T\approx 0.1m_{\text{i}}{V}_{\text{A},\text{in}}^2$ (Phan et al. 2014; Shay et al. 2014; Haggerty et al. 2015). All in all, the magnetic field energy of the plasma in the inflow acts as an energy reservoir for the kinetic energy in the outflow.

In RunS, the two reconnection sites at the nose and tail of the magnetosphere present different electron temperatures in the outflows. In Fig. 3, by showing the evolution of the electron temperature T_e = (T_e, ∥ + 2T_e, ⊥)/3 in the dipolar Y_MSO = 0 plane (the reconnection plane), we highlight the circulation of solar-wind electrons from the nose of the magnetopause to the tail finally entering the magnetosphere cavity. The electron temperature increases with respect to its upstream value, T_e, SW ≈ 20 eV, up to 100–200 eV in the nose, reaching larger values on the order of few keVs in the tail. This asymmetry in electron heating between nose and tail reconnection regions is due to the different inflow plasma parameters; while in the nose the magnetic energy is $m_{\text{i}}{V}_{\text{A},\text{in}}^2\approx 1$ keV, in the tail it is $m_{\text{i}}{V}_{\text{A},\text{in}}^2\approx 10$ keV, given the lower density and higher magnetic field amplitude in the lobes as compared to the magnetosheath. Thus, more magnetic energy is available for conversion into kinetic energy in the tail compared to the nose. This simple estimation explains why we observe electrons with higher energy in the tail (by around a factor of 10) with respect to the nose of the magnetosphere.

Magnetic reconnection allows for the injection of solar-wind electrons into the Hermean magnetosphere. As shown by the time evolution in Fig. 3, from the time t ≈ 10R/V_x, the nightside part of the magnetosphere starts to be populated with high-energy electrons ejected from the reconnection site in the tail around X_MSO ≈ −2R (such a planetward flow of plasma is also called a substorm, in analogy with Earth; Christon 1987). Once inside the Hermean magnetosphere, the motion of such an electron’s substorm is driven by the dipolar structure of the magnetic field. Similarly to Earth, in this region electrons bounce back-and-forth along closed magnetic field lines with period τ_B. Electrons also move around the planet following a longitudinal drift motion with a period of τ_i due to curvature drift. These motions have very a different timescale, τ_ce ≪ τ_B ≪ τ_i (here, τ_ce is the electron gyro-period as defined in Sect. 2), thus allowing us to treat them separately using the adiabatic approximation (Northrop 1963).

Using the adiabatic approximation, the bouncing motion of electrons along a given magnetic field line is described by the conservation of energy:

$$\begin{array}{c}\hfill E=\frac{1}{2}{m}_{\mathrm{e}}{v}_{\Vert }^2(s)+\mu B(s),\end{array}$$(1)

where $\mu =m\langle v_{\perp }^2\rangle /2B$ is the first adiabatic invariant. Here, brackets (⟨⟩) indicate the average over the gyro-motion, and s is the curvilinear coordinate along magnetic field lines. Similarly to the problem of a harmonic oscillator, the inversion points s_m are defined as those where the kinetic energy goes to zero, and thus

$$\begin{array}{c}\hfill B(s_{\mathrm{m}})=\frac{E}{\mu }=\frac{B(s_0)}{{sin}^2{\theta }_0},\end{array}$$(2)

where θ₀ is the initial particle pitch angle, θ₀ ≡ tan⁻¹(v_0,⊥/v_0,‖), computed at s₀ (the point where the magnetic field line crosses the magnetic equator). From Eq. (2) and using the value of the magnetic field amplitude at Mercury’s southern pole, B_south (on the order of 200 nT), we obtain the following loss-cone equation for electrons:

$$\begin{array}{c}\hfill |{\theta }_0|<{sin}^{-1}\sqrt{\frac{B(s_0)}{B_{\mathrm{south}}}}.\end{array}$$(3)

This equation defines a family of loss cones associated with different magnetic field lines (parameterized by s₀) varying with radius and local time (hereafter LT). Thus, spatial variations of the magnetic field amplitude B(s₀) in the equatorial plane induce variations in the electron pitch-angle distribution around the planet.

In RunS, we observe such a variation of electron distribution function around the planet as a function of local time. This is shown in Fig. 4 using three-dimensional spatial distributions of electrons with T_e, ⊥ > 500 eV (orange) and T_e, ∥ > 500 eV (purple). Fewer energetic electrons are observed in the dayside sector compared to the nightside one, and the few observed there have pitch angles close to 90°. This means that most of the electrons do not complete a full-drift orbit around the planet creating a stable continuous belt as on Earth, but rather bounce along magnetic field lines in the nightside before drifting dawnward and eventually falling on the planet surface. In the following, we refer to this partial circulation pattern as nightside-trapping.

Fig. 4. Three-dimensional view of high-energy electron population in simulation RunS. The three panels show three different views of the same quantities. Grid cells with values of perpendicular (parallel) temperatures greater than 500 eV are shown in orange (purple). We show an equatorial cut of the electron density in grayscale. Magnetic field lines are shown around the neutral point in the tail, and these lines are colored according to the value of Bx (blue and red for sunward and antisunward directions, respectively). Both quantities are computed at a time of 15.9 R/Vx, which is the same as in Figs. 1 and 2. Mariner10 and BepiColombo first Mercury flyby trajectories are shown using gray and blue solid lines, respectively.

To investigate such peculiar property of the Hermean environment, we look at the variation of the profiles B(s) with local time (LT) and radial distance, as shown in Fig. 5. The compression of the dipole magnetic field by the solar wind stretches the field lines on the nightside (blue-black lines in Fig. 5) and compresses the field lines on the dayside (yellow-orange lines in Fig. 5). Such a dayside compression increases the magnetic field amplitude B(s₀) at the subsolar equator, which in turn increases the loss-cone angle for electrons according to Eq. (3). The loss-cone angles resulting from Eq. (3) are reported in Figs. 5a–c for varying LT (0 corresponds to local midnight and 12 corresponds to local noon) and radial distance (from 1.1 to 1.5 R). For any radial distance close enough to the planet, the loss-cone angle steadily increases by moving from nightside (angle on the order of 45°) to dayside (angle on the order of 70°). This explains the nature of the electron nightside trapping observed in RunS. On top of that, the values of the loss-cone angles reported in Fig. 5 point out the higher probability of having trapped particles in the outer shells (≈1.5 R) rather than very close to the planet (≈1.1 R). This interpretation of adiabatic nightside trapping is consistent with the spatial distribution of high-energy electrons around the planet shown in Fig. 4.

Fig. 5. Magnetic field amplitude along magnetic field lines |B(s)| [nT] computed at different local times (LT). LT 12 corresponds to subsolar longitude (yellow) and LT 0 corresponds to local midnight (dark blue). Panels (a)–(c) correspond to different equatorial distances from the center of the planet: 1.1 R (a), 1.3 R (b), and 1.5 R (c). Resulting values for the loss-cone angle from Eq. (3) are also shown in each panel. Plots done at a time of 15.9 R/Vx, which is the same as in Figs. 1–4.

The loss-cone mechanism creates an inhomogenous distribution of high-energy electrons inside the Hermean magnetosphere (in the range R ≈ 1.1−1.5) with anisotropic energy distribution function (T_e, ⊥ > T_e, ∥). Such an electron population is expected to be the target of observation by satellites, able to measure in situ electron energy distribution functions (hereafter eEDFs). In the next section, we address this point by showing synthetic eEDFs sampled along Mariner10 and BepiColombo first Mercury flyby trajectories in our simulations.

5. Comparison between observations and simulations

The simulation results discussed so far are of particular interest for the interpretation of in situ electron observations at Mercury. To ease comparison with observations, we sampled the synthetic electron energy distribution functions (eEDFs) from our two simulations along Mariner10 and BepiColombo first-Mercury-flyby trajectories, as shown in Fig. 6.

Fig. 6. Synthetic electron energy distribution functions (eEDFs) sampled along Mariner10 (a)–(d) and BepiColombo (e)–(h) first Mercury flybys in our simulation boxes using northward IMF (left panels) and southward IMF (right panels) conditions. The plots are done using simulation results at time = 14.7R/Vx. The energy ranges and resolutions used to compute the eEDF are the same as those of the instruments PLS (onboard Mariner, top) and MEA1 (onboard BepiColombo, bottom). The horizontal axis shows the time from closest approach (CA) in minutes, and t = 0 is highlighted using vertical red dashed line. On the time axis, we highlight the bow shock and magnetopause crossings obtained by the two spacecraft in situ. These are highlighted using gray areas and vertical black dashed lines. The electron density was obtained by integrating the synthetic eEDF in energy at each given time.

To compute the energy spectra in Fig. 6 from our simulations, we select the macro-particles lying in the neighborhood of the spacecraft trajectory. In particular, we use a range of ±4 cells (equal to ≈ ± 0.06 R) in all directions in order to minimize the particle noise and mimic the effect of particle precipitation on the instruments in an acquisition time, δt. Moreover, we assume that the field of view of the instruments spans the whole 4π solid angle in an acquisition time of δt. Although this assumption is usually partially verified by the instruments, the main features of the spectra in Fig. 6 are also expected to hold in the case of limited field of view. To ease the comparison with in situ observations, the spectra in Figs. 6a,b have the same energy range and resolution of the Mariner10/PLS instrument (Ogilvie et al. 1974) and those in Figs. 6e,f of the BepiColombo/MEA1 (Saito et al. 2021) instrument operating in solar-wind mode during the flyby. Moreover, in Fig. 6 we report the bow shock and magnetopause crossing times obtained by the two spacecraft in situ at Mercury as discussed in Russell et al. (1988) (panels a–d) for Mariner10 and in André & Aizawa (priv. comm.) for BepiColombo (panels e–h). In the following, we discuss these two spacecraft’s flybys modeling results separately.

Mariner10 synthetic eEDF and electron densities obtained from our simulations, shown in Figs. 6a–d, present clear signatures of (i) inbound bow shock crossing at a time of t ≈ [ − 20, −18], (ii) inbound magnetopause crossing at a time of t ≈ −7, (iii) outbound magnetopause crossing at a time of t ≈ [6, 7], and (iv) outbound bow shock crossing at a time of t ≈ [12, 14] (times are given in minutes to closest approach). These crossing times are consistent with those observed by Mariner10 (Russell et al. 1988), and the latter are and shown in Fig. 6 using vertical black dashed lines and gray areas. A more precise comparison of the crossing times between observations and simulations would require an ad hoc simulation initialized with the upstream solar-wind parameters observed by Mariner10. A further step of this kind will be addressed in future works. Here, we mainly focus on the qualitative features of the eEDF in the Hermean magnetosphere. Inside the magnetosphere (t ≈ [ − 7, +7]) in both simulations the plasma density is strongly depleted, and we observe different signatures in the two runs.

In RunN, electrons with energies up to a few keV are observed around the magnetopause (the mechanism accelerating these electrons, however, remains unclear) and no electrons are observed inside the magnetosphere.

In RunS, electrons with energies up to tens of keV (well above the instrumental cutoff of 687 eV) are encountered around and after closest approach (t ≈ [0, +7]). These high-energy electrons are produced by magnetic reconnection in the tail, as discussed in Sect. 4.

This result supports the interpretation of Mariner10/PLS observations outlined in past works (Ogilvie et al. 1974; Christon 1987), stating that Mariner10 observed substorm injected electrons around the closest approach with energies well above the instrumental cutoff. It can be further used to interpret the recent observations by BepiColombo on Mercury.

BepiColombo synthetic eEDF and electron densities obtained from our simulations, shown in Figs. 6e–h, present similar signatures to those of Mariner10. However, two main differences are observed. First, for BepiColombo the inbound bow shock crossing is less sharp than for Mariner10. This is due to the large distance from the planet of the spacecraft at the encounter of the inbound bow shock (see the trajectory in Fig. 4). BepiColombo’s first flyby trajectory makes plasma observations more asymmetric between inbound and outbound as compared to Mariner10. Second, signatures of high-energy electrons around closest approach are fainter for BepiColombo than for Mariner10. This effect can be understood again in terms of trajectories. From Fig. 4, we note that BepiColombo passes beneath the high-energy region of the nightside, while Mariner10 passes right through such a region (shown in purple and orange in Fig. 4).

We stress that a more precise, quantitative comparison between observations and simulations would require a good knowledge of the upstream solar-wind parameters that define the precise values of, among others, the time of the bow shock and magnetopause crossings, the energy of the electrons inside the magnetosphere, and the density of electrons in the magnetosheath. Thus, the results of this section are limited to the qualitative main features observed by Mariner10 and BepiColombo first Mercury flybys.

A comparison of the synthetic eEDF between RunN and RunS demonstrate that by changing only one parameter (the IMF direction), a sudden change in the energy and density of electrons inside the magnetosphere is induced, as shown in Fig. 6. We suggest that the Hermean environment responds to such changes in the solar wind by a corresponding reconfiguration as from RunN to RunS or viceversa. The characteristic timescale for such a reconfiguration is expected to be comparable to the Dungey cycle of τ ≈ 10 R/V_x ≈ minutes. Thus, it is only in the case of stable solar-wind conditions on timescale longer than a few minutes that observational signatures in the eEDF inside the magnetosphere can be linked to the upstream IMF direction.

6. Discussion and conclusions

This work represents a first step towards global, fully kinetic modeling of planetary magnetospheres in which the electron dynamics is included self-consistently from the global planetary scale down to the gyro-radius scale. The high computational needs of this model impose strong constraints on the possible scale separation between planet, ion, and electron scale length, as discussed at the end of Sect. 2. Scaling down the ion-to-electron mass ratio and the plasma-to-cyclotron frequency ratio remains a classic procedure in fully-kinetic simulations, the effects of which have been extensively discussed in past works (Bret & Dieckmann 2010; Le et al. 2013; Lavorenti et al. 2021). On the other hand, using a smaller planet radius of ≈230 km instead of the realistic Mercury radius of ≈2400 km may have a significant impact on a simulation’s results (in our work, the scale separation between planet radius and ion gyro-radius remains of about a factor of 10).

First, scaling down the planet radius can have an impact on the global shape of the magnetosphere. To avoid this problem, the magnetic moment of the planet is scaled down proportionally to R³ in order to preserve the shape of the magnetosphere. Therefore, the global shape of the interaction region is preserved, and in particular the bow shock d_BS/R and magnetopause d_MP/R stand-off distance in units of the planet radius are preserved, as discussed in Sect. 3. Moreover, adiabatic particle acceleration processes arising in the dipolar magnetic field of the planet (such as betatron and Fermi acceleration) are unaltered by the planet rescaling as well. This is particularly important for the discussion in Sect. 4 on the loss-cone angle around the planet.

Second, scaling down the planet radius can have an impact on magnetic reconnection in the tail in the case of southward IMF. Magnetic reconnection induces the formation of a region of drifting plasma close to the neutral point, the so-called diffusion region, where particles are no longer frozen in the magnetic field. Depending on the species, this region is called the ion (electron) diffusion region and it extends for about ten ion (electron) inertial lengths around the neutral point. By scaling down the planet radius, the ion (electron) diffusion region in the tail moves closer to the planet surface. In our case, given the scaled-down planet we are using (R = 5.5 d_i, where $d_{\mathrm{i}}\equiv c/{\omega }_{\mathrm{pi}}\equiv c\sqrt{m_{\mathrm{i}}/4\pi n_{\mathrm{SW}}{e}^2}$ is the ion inertial length), the ion (electron) diffusion region is marginally (largely) separated from the planet. Thus, we expect ion dynamics in the outflow to be affected by the scaling down, at least partially, while this should not be the case for electrons. In particular, we expect the signatures of high-energy electrons observed in our simulations as a product of magnetic reconnection to hold true when using a real size planet.

Although we expect these two classes of plasma processes to be well-reproduced using the scaled-down parameters in our simulations, particular care must be taken when analyzing specific plasma quantities. Indeed, due to the rescaling of the plasma-to-cyclotron frequency ratio, electrons with energies above few tens of keV fall out of the range of validity of our model since their speed becomes ultra-relativistic (by reducing the ratio ω_pe/ω_ce, we are de facto reducing the value of the light speed c in our simulations). Moreover, by scaling down the ion-to-electron mass ratio, we reduce the ratio between the electron thermal speed and the flow speed (v_the/V_SW) from 5 to around 1.2. This means that the solar-wind electrons interacting with the planet in our simulations are injected with a slightly more anisotropic energy distribution function as compared to the real case. The relevance of these caveats to our simulation results remains unknown and should be further investigated in future works. We stress that global, fully kinetic simulations of planetary magnetospheres using realistic plasma parameters remain prohibitive using the present state-of-the-art HPC facilities.

Since this work represents a first step toward a fully kinetic global modeling of the Hermean environment, we chose to use a simple as possible realistic solar-wind configuration. Thus, we chose a purely northward or southward IMF. This choice facilitates the analysis of the reconnection sites and outflows. A more realistic configuration would require a strong B_x component in the IMF. Such an in-plane component, typically found on Mercury, tends to create a foreshock region (absent from our simulations) and drives strong north-south asymmetries in the magnetosphere. A more realistic modeling of this kind will be addressed in future works.

To conclude, the results of this work show that high-energy electrons (up to tens of keV) are generated in the magnetotail of Mercury in the case of southward IMF. Such electrons form as a consequence of magnetic reconnection in the tail neutral line located at X ≈ −2 R and move toward the planet drifting dawnward. When reaching the internal shells of the magnetosphere, the motion of these electrons (which can be treated as adiabatic) is strongly affected by loss-cone precipitation on the planet surface. We find that almost all drifting electrons are lost on the planet before completing a full longitudinal drift orbit. Therefore, the highest concentration of electrons in the range of hundreds of eV to some keV is on the nightside, especially in the post-midnight sector. This result supports electron observations by the Mariner10/PLS instrument, showing high-energy electrons around closest approach during its first Mercury flyby. Such in situ observations are in agreement with our simulation with the southward IMF. Moreover, we present synthetic electron data useful for the ongoing and future interpretation of BepiColombo/MEA first-Mercury-flyby observations. A detailed comparison between our simulations and MEA observations for this flyby will be done in the near future once the data have been fully calibrated. Finally, we envision a characterization of this energetic electron population by BepiColombo during its nominal science phase.

Acknowledgments

This work was granted access to the HPC resources at TGCC under the allocation AP010412622 and A0100412428 made by GENCI via the DARI procedure. We acknowledge the CINECA award under the ISCRA initiative, for the availability of high performance computing resources and support for the project IsC93. We acknowledge the support of CNES for the BepiColombo mission. Part of this work was inspired by discussions within International Team 525: “Modelling Mercury’s Dynamic Magnetosphere in Anticipation of BepiColombo” at the International Space Science Institute, Bern, Switzerland. We acknowledge support by ESA within the PhD project “Global modelling of Mercury’s outer environment to prepare BepiColombo”. Analysis of BepiColombo and Mariner10 data was possible thanks to the AMDA (http://amda.irap.omp.eu) and CLWeb (http://clweb.irap.omp.eu) tools.

References

All Figures

Fig. 1. Overview of large-scale structure and magnetosphere boundaries in our simulations. Panels a–f show the results of RunN. Panels a–c show dipolar YMSO = 0 cuts of the ion density ni (a), magnetic field amplitude |B| (b), and plasma current density along the y-axis Jy (c). Panels d–f show equatorial ZMSO = 0 cuts of the same quantities for RunN. Panels g–l show the results of RunS in the same format as those of RunN. All plots represent a time of t = 15.9 R/Vx. The average bow shock (dashed black line) and magnetopause (solid black line) profiles found by Winslow et al. (2013) using MESSENGER observations are added. The white region around the planet corresponds to grid cells with zero macro-particles.

In the text

Fig. 2. Cut of results of our two simulations RunN (left panels) and RunS (right panels) along the subsolar line (x-axis). Panels a–b: proton and electron densities. Panels c–d: magnetic field components and amplitude. Panels e–f: current density along YMSO of protons, electrons, and total. Panels g–h: magnetic Pmag ≡ B2/2μ0, ram $P_{\text{ram}}\equiv m_{\text{i}}{n}_{\text{i}}{V}_{x,\text{i}}^2+m_{\text{e}}{n}_{\text{e}}{V}_{x,\text{e}}^2$, and total pressures. These cuts were performed at the same time as Fig. 1 (t = 15.9 R/Vx). We highlight the bow shock and magnetopause stand-off distances obtained from our model (solid red vertical lines) and the ones obtained using the average model of Winslow et al. (2013) constrained by MESSENGER observations (dashed red vertical lines).

In the text

Fig. 3. Time evolution of electron temperature in nose and tail reconnection regions. Dipolar YMSO = 0 cut of simulations RunN (top panels) and RunS (bottom panels) at different times (time advancing from left to right). We show the color map of the electron temperature Te = (2T⊥, e + T∥, e)/3 in units of [eV] in logarithmic scale. On top of that, we superpose the in-plane magnetic field lines (gray solid lines). The white region around the planet corresponds to grid cells with a number of macro-particles too small to reliably compute the value of Te.

In the text

Fig. 4. Three-dimensional view of high-energy electron population in simulation RunS. The three panels show three different views of the same quantities. Grid cells with values of perpendicular (parallel) temperatures greater than 500 eV are shown in orange (purple). We show an equatorial cut of the electron density in grayscale. Magnetic field lines are shown around the neutral point in the tail, and these lines are colored according to the value of Bx (blue and red for sunward and antisunward directions, respectively). Both quantities are computed at a time of 15.9 R/Vx, which is the same as in Figs. 1 and 2. Mariner10 and BepiColombo first Mercury flyby trajectories are shown using gray and blue solid lines, respectively.

In the text

Fig. 5. Magnetic field amplitude along magnetic field lines |B(s)| [nT] computed at different local times (LT). LT 12 corresponds to subsolar longitude (yellow) and LT 0 corresponds to local midnight (dark blue). Panels (a)–(c) correspond to different equatorial distances from the center of the planet: 1.1 R (a), 1.3 R (b), and 1.5 R (c). Resulting values for the loss-cone angle from Eq. (3) are also shown in each panel. Plots done at a time of 15.9 R/Vx, which is the same as in Figs. 1–4.

In the text

Fig. 6. Synthetic electron energy distribution functions (eEDFs) sampled along Mariner10 (a)–(d) and BepiColombo (e)–(h) first Mercury flybys in our simulation boxes using northward IMF (left panels) and southward IMF (right panels) conditions. The plots are done using simulation results at time = 14.7R/Vx. The energy ranges and resolutions used to compute the eEDF are the same as those of the instruments PLS (onboard Mariner, top) and MEA1 (onboard BepiColombo, bottom). The horizontal axis shows the time from closest approach (CA) in minutes, and t = 0 is highlighted using vertical red dashed line. On the time axis, we highlight the bow shock and magnetopause crossings obtained by the two spacecraft in situ. These are highlighted using gray areas and vertical black dashed lines. The electron density was obtained by integrating the synthetic eEDF in energy at each given time.

In the text