© ESO 2018

1 Introduction

Mercury has a weak global magnetic field of internal origin that was first discovered by the Mariner 10 spacecraft in 1975–1976 through three flybys of Mercury (e.g. Ness et al. 1974). Later, The MErcury Surface, Space ENvironment, GEochemistry, and Ranging (MESSENGER) spacecraft (Solomon et al. 2001) provided a more accurate determination of the internal magnetic field of Mercury using two flybys in 2008, known as M1 and M2 (e.g. Anderson et al. 2008, 2010), and nearly four years of magnetic field observations around Mercury from early 2011 until 2015 (e.g. Anderson et al. 2011a, 2012; Johnson & Hauck 2016). The global planetary field of Mercury has been estimated as a single dipole with a magnetic moment of 195 ± 10 nT displaced northward by 484 ± 11 km (~0.2 R_M) from the centre of the planet, where R_M = 2440 km is the radius of Mercury (Anderson et al. 2011a, 2012). The magnetic dipole moment of Mercury is directed southward, and its axis is tilted <3° from Mercury’s spin axis (e.g. Anderson et al. 2011a). The interplanetary magnetic field (IMF) and the supersonic flow of the solar wind plasma are continuously interacting with the intrinsic magnetic field of Mercury, resulting in the formation of a “mini-magnetosphere”, that is qualitatively similar to Earth’s magnetosphere. In general, this interaction forms a collision-less bow shock that decelerates and diverts the solar wind plasma and magnetic fields around the magnetospheric obstacle of Mercury (e.g. Anderson et al. 2011a; Masters et al. 2013), forms a magnetosheath with heated plasma between the bow shock and the magnetopause (e.g. Fairfield & Behannon 1976; Raines et al. 2011), and creates an extended magnetotail with a central current sheet (e.g. Slavin et al. 2010, 2012a; Sun et al. 2015; Poh et al. 2017).

Despite the similarities in the overall structures of the magnetospheres of Mercury and the Earth, magnetic field observations by the MESSENGER magnetometer (MAG; Anderson et al. 2007) have indicated that Mercury’s magnetospheric phenomena occur on spatial and temporal scales that are very different to those at Earth and any other magnetized planet in the solar system (e.g. Slavin et al. 2009, 2010; DiBraccio et al. 2013; Raines et al. 2015). For example, the magnetic reconnection rate is estimated to be ~0.15 at Mercury, which is nearly an order of magnitude higher than that observed at Earth (DiBraccio et al. 2013). In addition, the Dungey cycle, the re-circulation of the solar wind energy between the dayside magnetopause and nightside magnetotail reconnection sites, has been estimated to be of the order of ~2 min at Mercury. This is nearly 30 times faster than the Dungey cycle at Earth and about five orders of magnitude faster than at Jupiter (Slavin et al. 2009).

These differences arise mainly because of the weak intrinsic magnetic field of Mercury, which is over two orders of magnitude weaker on the surface of Mercury than on the Earth, and because of the solar wind plasma intensity and the strength and orientation of the IMF near the orbit of Mercury. The solar wind plasma density near Mercury (on average ~30 cm⁻³) is approximately five times higher than that near the Earth. Therefore, the dynamic pressure is nearly five times stronger at Mercury. This high dynamic pressure of the solar wind and weak planetary magnetic fields move the dayside magnetospheric boundaries, i.e. the bow shock and magnetopause, approximately eight times closer to Mercury than to Earth (e.g. Raineset al. 2015, and references therein). In addition, the IMF strength near the orbit of Mercury (typically ~20 nT) is nearly four times stronger than near the Earth. Hence, the solar wind plasma β, a ratio of the solar wind thermal pressure to the magnetic pressure, and Alfvén Mach number M_A, a ratio of the solar wind to the Alfvén speed, are lower (β ≈ 0.5 and M_A ≈ 4.0) compared to those for the solar wind around the Earth (β ≈ 1.0 and M_A ≈ 7.0). The low plasma β and low M_A as well as a highly dominating component of the IMF along the solar wind flow near Mercury, as opposed to the nearly 45° IMF relative to the solar wind flow direction at Earth, have been suggested as the main reasons for the high magnetic reconnection rate at Mercury, which, in contrast to the Earth, is independent of the magnetic shear angle (DiBraccio et al. 2013).

These fundamental differences, together with the observed fast and transient phenomena in Mercury’s magnetosphere (e.g. a few seconds flux transfer events (Slavin et al. 2012b) and plasmoids (Slavin et al. 2012a; DiBraccio et al. 2015)) indicate that the magnetosphere of Mercury is highly dynamic, sensitive, and responsive to the solar wind plasma and IMF variations, much more than any other magnetized planet in the solar system (Burlaga 2001; Slavin et al. 2008, 2009). Solar wind variations contribute considerably to altering the structure of Mercury’s magnetosphere (e.g. Slavin et al. 2009; Anderson et al. 2011a; Varela et al. 2015). Therefore, to understand the structure of Mercury’s magnetosphere and its response to the solar wind plasma and IMF variations, we need to understand the interaction between the solar wind and Mercury’s magnetosphere and distinguish between the contributions from external and internal magnetic sources (e.g. Raines et al. 2015; Johnson & Hauck 2016; James et al. 2017).

Due to the lack of an upstream solar wind plasma monitor, it is difficult to estimate the solar wind parameters and their variations during a passage of a spacecraft through the dynamic magnetosphere of Mercury. This is even more pronounced for a spacecraft like MESSENGER due to the limited field of view of its plasma instrument, the Fast Imaging Plasma Spectrometer (FIPS; Zurbuchen et al. 1998; Andrews et al. 2007), especially in directions transverse to the Mercury–Sun line which limits observations of the solar wind plasma when MESSENGER is outside Mercury’s magnetosphere (e.g. Zurbuchen et al. 2008; Raines et al. 2011). Thus, there is no complete information about the solar wind plasma parameters, flow direction, and their variations, neither when MESSENGER is outside nor when it is inside Mercury’s magnetosphere (e.g. Korth et al. 2011; Baker et al. 2013; Winslow et al. 2013; Dewey et al. 2015). These challenges also hold for the future ESA/JAXA mission to Mercury, BepiColombo (Benkhoff et al. 2010), and its plasma packages including SERENA on Mercury Planetary Orbiter (MPO; Orsini et al. 2010) and MPPE on Mercury Magnetospheric Orbiter (MMO; Saito et al. 2010).

Different methods have been applied to compensate for the lack of an upstream solar wind plasma monitor and to fill in the gaps of unobserved fractions of velocity space distributions inside the magnetosphere of Mercury. For example, Korth et al. (2011, 2012) have estimated plasma pressure near the plasma sheet by maintaining pressure balance between the sheet plasma and the observed magnetic field near the equator on the nightside of Mercury’s magnetosphere. They estimated the average pressure near plasma sheet is ~1.45 nPa normalized to Mercury’s heliocentric distance of 0.39 AU (Korth et al. 2012). As another example, Winslow et al. (2013) used plasma parameters predicted by the WSA-ENLIL model for Mercury (Baker et al. 2013) and magnetic field observations by MESSENGER to estimate the response of the dayside magnetopause and bow shock to the solar wind dynamic pressure and Alfvén Mach number. They found that the average subsolar stand-off distance of the magnetopause is 1.45 R_M for a mean solar wind dynamic pressure of 14.3 nPa and the average subsolar distance of the bow shock is 1.96 R_M for a mean Alfvén Mach number of 6.6.

Here weuse a three-dimensional self-consistent hybrid model of plasma (kinetic ions and charge neutralizing fluid electrons) that runs onGraphics Processing Units (GPUs; Fatemi et al. 2017). We use our model to infer the upstream solar wind dynamic pressure and Alfvén Mach number from magnetic field observations along the trajectory of MESSENGER inside the magnetosphere of Mercury. We compare our simulation results with MESSENGER observations that have been published, modelled, and analysed in detail before. We also compare the inferred solar wind dynamic pressure from our simulations with those predicted by the WSA-ENLIL model upstream of Mercury. We show that there is a good agreement between our simulation results and MESSENGER observations. Finally, we show that our model can be used as an upstream solar wind plasma monitor and can provide estimates for plasma parameters inside Mercury’s magnetosphere. This has direct implications for observations by MESSENGER and the future ESA/JAXA mission to Mercury, BepiColombo.

2 Model

We use the AMITIS code, the first GPU-based three-dimensional self-consistent hybrid plasma model that uses a single CPU–GPU pair (Fatemi et al. 2017). This model uses only a single CPU and a single GPU, yet runs at least 10 times faster and is more energy and cost efficient than its parallel CPU-based predecessors (Fatemi et al. 2017). In this model, the ions are charged macro-particles and the electrons are a mass-less charge neutralizing fluid. The Lorentz force and the equation of motion are used to advance particle trajectories in time. The electric field, E, is directly calculated from the electron momentum equation. Faraday’s law, ∂B∕∂t = −∇ × E, is used to advance the magnetic field B in time using an implicit-explicit scheme explained in detail in Fatemi et al. (2017). The model has been successfully applied to study plasma interactions with the Moon (Fatemi et al. 2017) and with the asteroid 16 Psyche (Fatemi & Poppe 2018).

In our model we use a Mercury Solar Orbital (MSO) coordinate system centred at Mercury’s centre of mass, where the + x axis is pointing to the Sun, the + y axis is opposite to the orbital motion of Mercury and points toward dusk, and the + z axis is pointing to the north (normal to the xy-plane) and completes the right-handed coordinate system. In addition, we also use a Mercury Solar Magnetospheric (MSM) coordinate system centred at Mercury’s dipole moment (Anderson et al. 2012) to analyse and explain the location of the bow shock and magnetopause in our simulations. The difference between the MSO and MSM coordinate systems is the location of the origin of these coordinate systems.

2.1 MESSENGER orbit selection

To infer the solar wind plasma parameters upstream of Mercury and in order to compare our simulations with observations, we have selected three orbits of MESSENGER that have previously been analysed and compared with hybrid simulations. These orbits include the first Mercury flyby on 14 January 2008 known as the M1 flyby, and two regular orbits on 23 April 2011 (DOY 113 between 15:00 and 22:00, hereafter D113) and 01 July 2011 (DOY 182 between 05:00 and 11:00; hereafter D182). The M1 flyby has been studied before (e.g. Slavin et al. 2008; Anderson et al. 2010, 2011a; Raines et al. 2011) and compared with hybrid simulations (Müller et al. 2012). The D113 and D182 orbits have been compared with simulations (Richer et al. 2012; Herčík et al. 2016) and compiled into a statistical analysis of Mercury’s magnetosphere (Winslow et al. 2013).

A portion of the trajectory of every orbit is shown in Fig. 1 in a cylindrical MSO coordinate system. Only the M1 flyby passed near the equator (xy-plane), while the D113 and D182 orbits, similar to other nominal orbits of MESSENGER in 2011, passed over the poles with a closest approach of ~200 km over ~60° northern latitude of Mercury (Solomon et al. 2007). As shown in Fig. 1, the selected orbits cover different areas of the magnetosphere including the equatorial region (M1), the Sun-midnight plane (D182), and an oblique angle from the midnight meridian (D113). These assure us that our model–data comparison does not only focus on a specific magnetospheric region and its associated phenomenon, but has also been validated against different magnetospheric locations based on the available MESSENGER observations.

Fig. 1 Portion of MESSENGER’s trajectory around Mercury during the M1 flyby (red line) and two nominal orbits of MESSENGER on 23 April 2011, D113 (green line), and on 01 July 2011, D182 (blue line) in a cylindrical MSO coordinate system. The arrows show the direction of MESSENGER’s motion along each orbit. Superimposed are the approximate locations of the magnetopause (thick dashed line) and bow shock (thin dashed line) obtained from Winslow et al. (2013).

2.2 Inverse problem approach

Due to the lack of an upstream solar wind monitor and no direct observation of the solar wind plasma by MESSENGER, there is no complete information about solar wind plasma parameters, i.e. density, flux, thermal speed, and dynamic pressure, for different passages of MESSENGER throughout Mercury’s magnetosphere. Therefore, we take an inverse problem approach using the AMITIS code to infer the solar wind plasma parameters upstream of Mercury based on magnetic field observations by the MAG instrument on MESSENGER. For every orbit shown in Fig. 1, we take the average of the observed magnetic fields outside the magnetospheric disturbances of Mercury and apply it to our model as the only known parameter in our simulations (for further details, see Sect. 2.3). Then we perform ~40–50 simulations for each orbit with various solar wind dynamic pressures, and compare the magnetic fields from our simulations with those obtained from MAG instrument along the trajectory of MESSENGER to find the best agreement between our simulations and observations for the location of the magnetospheric boundaries (bow shock and magnetopause), and the overall intensity and orientation of the magnetic fields (for further details, see Sects. 2.3 and 2.4). This approach not only estimates the solar wind dynamic pressure, and consequently the Alfvén Mach number, upstream of Mercury, but also provides a general understanding of the plasma environment inside Mercury’s magnetosphere, as well as detailed information along the trajectory of MESSENGER.

2.3 Simulation parameters and assumptions

We assume that Mercury is a spherical object of radius R_M = 2440 km without an exosphere and that its surface is a perfect plasma absorber. We place a southward oriented magnetic dipole along the -z axis with a magnetic moment of 195 nT, displaced 484 km northward in the MSO coordinate system (Anderson et al. 2010, 2011a). We ignore the small tilt of the magnetic moment from Mercury’s spin axis. Since the intrinsic magnetic moment of Mercury has been estimated using observations (Anderson et al. 2010, 2011a), we consider it a constant in our model and do not change it in our simulations.

We use a simulation domain of size –7R_M ≤x ≤ +6 R_M and –10R_M ≤ (y, z) ≤ +12 R_M with a regular-spaced Cartesian grid with cubic cells of size 200 km (~0.08 R_M). We use 16 macro-particles (only protons) per cell at the inflow boundary (x = +6 R_M) where the solar wind enters the simulation box. Each macro-particle, also known as a super-particle, represents a large number of real particles to make particle simulations computationally effective; nevertheless, the charge-per-mass ratio of each macro-particle used to solve the equation of motion is equal to the charge-per-mass ratio of realparticles (Birdsall & Langdon 2005). We advance particle trajectories using a time step of Δ t = 0.001 s, which is nearly 3 × 10⁻⁴ of the solar wind proton gyroperiod away from Hermean magnetospheric disturbancesand is about 3 × 10⁻³ of a proton gyroperiod near Mercury’s magnetic poles. This small time step assures that the gyromotion of the solar wind protons is fully resolved in our simulations. While the inflow (x = +6 R_M) and the outflow (x = –7 R_M) boundaries are perfect plasma absorbers, the boundaries along the y- and z-axes are assumedto be periodic for particles (i.e. particles are transported to the opposite side of the simulation domain) and electromagnetic fields (i.e. a copy is made of the electromagnetic fields that occur at one side of the simulation domain and are injected at the other side).

For the selected orbits shown in Fig. 1, the IMF strength and orientation remains almost steady for nearly one hour before (after) MESSENGER moved into (out of) the magnetosphere with small variations. In addition, Winslow et al. (2012) used MESSENGER observations and showed that the one-hour average is suitable for the dominant component of the IMF (i.e. B_x), but that the other components may vary within shorter time periods (e.g. 30 min) (Winslow et al. 2012, 2013). However, for the selected orbits in this study, especially for D113 and D182, we did not find considerable variations between one-hour and 30-min averages of the IMF. Therefore, for the orbits analysed here, we take a one-hour average of magnetic field data observed by the MAG instrument before (after) MESSENGER moved into (out of) Mercury’s magnetosphere, except for the z-component of the IMF, B_z, during the M1 flyby. Taking a one-hour average of B_z before the inbound and after the outbound bow shock crossing during the M1 flyby gives a southward IMF with B_z ≈ –1.0 nT. However, the previous analyses of the M1 flyby (e.g. Slavin et al. 2008, 2010; Anderson et al. 2011a) considered B_z ≈ +4 nT by taking the average of the magnetic fields for a shorter period than one hour during the outbound bow shock crossing. Therefore, for consistency with previous data analyses, we also consider a northward IMF with B_z = +4 nT for the M1 flyby. The averaged values are listed in Table 1 and we apply them as inputs into our model and keep them constant at the inflow boundary of the simulation box. The IMF orientation and strength during M1 and D113 did not show considerable changes before the inbound and after the outbound bow shock crossings, during D182; however, the IMF orientation and strength changed. Therefore, as listed in Table 1, we consider two separate sets of simulations for D182: (1) a southward IMF with a considerable component along the -y axis before the inbound bow shock crossing and (2) a northward IMF after the outbound bow shock crossing. We note that the outbound B_x component is nearly 30% stronger than the inbound.

In contrast to the IMF, the upstream solar wind plasma parameters, including plasma density n_sw, velocity u_sw, and temperature T_sw, are unknowns. The solar wind dynamic pressure P_dyn = , where m_i is the solar wind proton mass, and the Alfvén Mach number M_A = are the main parameters that control the shape and structure of the magnetospheric boundaries (e.g. Slavin et al. 2008, 2009; DiBraccio et al. 2013; Winslow et al. 2013 and they are a function of variables including B, n_sw, and u_sw. Since the magnetic field B is already known from observations and listed in Table 1, we only consider n_sw and u_sw as the variables in our simulations. For simplicity and to decrease the degrees of freedom, we assume that the solar wind temperature for both ions and electrons is constant, T_i = T_e ≈ 12 eV, and is approximately equal to the average solar wind proton temperature near the orbit of Mercury (e.g. Marsch et al. 1982). Although n_sw and u_sw are unknowns and are considered as free parameters in our simulations, we only select them within the expected and/or probable ranges for the solar wind plasma near the orbit of Mercury, i.e. 16 ≲ n_sw ≲ 120 cm⁻³ and 270 ≲|u_sw|≲ 650 km s⁻¹ (e.g. Winslow et al. 2013).

Although Mercury has an eccentric orbit around the Sun (eccentricity ~0.2), and its Keplerian speed varies considerably between perihelion (~56 km s⁻¹) and aphelion (~38 km s⁻¹) (e.g. Murchie et al. 2014), for simplicity we take its average orbital speed (50 km s⁻¹) and compensate for it as a downward component of the solar wind plasma flow in all simulations presented in this study. Therefore, we assume that the solar wind flows along the -x axis with a fixed 50 km s⁻¹ component along the +y axis in the MSO coordinate system.

In addition, Mercury has a large conductive core, that induces magnetic fields from a time-varying IMF or from large dynamic pressure variations in the solar wind (e.g. Smith et al. 2012; Hiremath 2012; Hauck et al. 2013; Johnson & Hauck 2016). In our simulations we assume that Mercury has a uniform resistive interior with resistivity η = 10⁷ Ω × m. Since the solar wind plasma and IMF orientation and strength remain constant at the inflow boundary of the simulations presented here, and since Mercury’s mantle has very low conductivity, no electromagnetic induction is generated by the interior of Mercury over a constant solar wind and IMF. Therefore, the uniform resistivity assumption for the interior of Mercury, although crude, is a valid approximation in this study.

Since the AMITIS has been intentionally developed and optimized to run on a single CPU–GPU pair, every single hybrid simulation run for Mercury for the simulation domain explained here takes nearly 30–40 h to reach 200 s. This high-performance tool enables us to run simulations simultaneously up to the number of available GPUs on our work laptops, office desktops, and super-computers equipped with GPUs. With our currently available resources, we can run on average about 10 simulation runs simultaneously. With the help of this high-performance tool, we have made 25–50 simulation runs for every orbit shown in Fig. 1 and applied different solar wind plasma densities and velocities to find the best agreement between our simulations and the magnetic field observations. The simulation parameter ranges are explained and motivated in more detail in Sect. 3. For every simulation we keep the solar wind plasma and IMF constant at the inflow boundary and run the simulation to ~200 s. Then we compare the magnetic fields and the location of the different magnetospheric boundaries obtained from each simulation with those observed by the MAG instrument along the trajectory of MESSENGER. The duration of each simulation run (200 s) is equivalent to the completion of over 50 solar wind proton gyrations, and nearly two Dungey cycles inside the magnetosphere, which is long enough for the development of the entire magnetosphere. It is also equivalent to the time required for a solar wind with average velocity of 370 km s⁻¹ to sweep over our simulation box more than two times.

2.4 Determination of the magnetospheric boundaries from simulations

Winslow et al. (2013) made a survey throughout MESSENGER magnetic field observations from 23 March 2011 to 19 December 2011. For every orbit within this period, they visually inspected the time that MESSENGER crossed Mercury’s magnetospheric boundaries, which are generally the bow shock and magnetopause crossings before and after the magnetospheric transients denoted as the inbound and outbound crossings. Due to the solar wind variations, multiple crossings of each boundary is highly probable, which has been considered by recording the first and last boundary encounters (Winslow et al. 2013). The D113 and D182 orbits, shown in Fig. 1, are within the period analysed by Winslow et al. (2013). For the M1 flyby the magnetospheric boundary crossings were examined by Slavin et al. (2009), Anderson et al. (2011a), and Raines et al. (2011).

Similar to Winslow et al. (2013), we categorize the magnetospheric boundary crossings into four events: the inbound and outbound bow shock, and the inbound and outbound magnetopause crossings. Then we calculate the electric currents in our simulations from the general Ampère’s law, . We use the intensity and direction of the currents to determine the location of the magnetospheric boundaries from simulations. Then we estimate the time at which each boundary crossing has occurred along the trajectory of MESSENGER and compare them with those estimated by Winslow et al. (2013) and Slavin et al. (2009). Finally, we calculate the time offset between the estimated boundaries from our simulations and the mid-point location of the magnetospheric boundaries obtained from observations to find the best agreement between our simulations and observations.

3 Results

Here we present our hybrid simulation results for the solar wind plasma interaction with Mercury for the orbits shown in Fig. 1. For every orbit, we take the averaged IMF listed in Table 1, run a number of simulations for different solar wind plasma density and velocity, and compare magnetic fields from our simulations with those observed by MESSENGER. Here our primary intention is to estimate the upstream solar wind dynamic pressure and Alfvén Mach number that result in the best agreement between our hybrid simulations and the magnetic field observations; our main focus is the magnetospheric boundary determination. We also study the global structure of Mercury’s magnetosphere during the selected orbits.

3.1 Orbit D113: 23 April 2011

On 23 April2011, MESSENGER encountered the magnetosphere of Mercury twice (each orbit takes ~12 h). While inthe solar wind and away from magnetospheric disturbances, the MAG instrument did not observe any large variations inthe IMF direction and intensity during this day. Therefore, we present our simulation results for this orbit first and only focus on the second encounter, which occurred between 17:00 and 20:30 approximately, mainly because this period has been studied before and compared with hybrid simulations (Richer et al. 2012). As shown in Table 1, the averaged IMF for this orbit mainly lies on the equatorial plane and makes a nearly 45° angle with the solar wind plasma flow direction with a minor northward component.

In order to estimate the upstream solar wind dynamic pressure for this orbit we performed over 40 simulation runs for the upstream solar wind plasma density ranging between 16 and 26 cm⁻³ with a span of 2 cm⁻³, and for the solar wind velocity along the -x axis ranging between 270 and 340 km s⁻¹ with a span of 10 km s⁻¹. As explained in Sect. 2.3, we always assume 50 km s⁻¹ for the solar wind velocity along the +y axis to account for the Keplerian speed of Mercury in the MSO coordinate system. These parameters cover the dynamic pressure range between ~2.0 and ~5.0 nPa, Alfvén Mach number between ~2.7 and ~4.3, and plasma β between ~0.21 and ~0.33. We selected this range because our pre-analysis studies (not shown here) suggested that there is a generally good agreement between our simulations and observations for these parameter ranges of the solar wind.

3.1.1 Magnetospheric boundary determination

As explained in Sect. 2.4, we use our simulations and compare magnetospheric boundary crossings (i.e. bow shock and magnetopause) obtained from each of them along the orbit of MESSENGER with those estimated by Winslow et al. (2013) for D113 from MESSENGER magnetic field observations. This comparison enables us to estimate the upstream solar wind plasma dynamic pressure and Alfvén Mach number that resulted in the observed location of the magnetospheric boundaries. Moreover, this analysis also helps us to better understand how the bow shock and magnetopause respond to different upstream solar wind plasma parameters.

Figure 2 shows the time difference between the magnetospheric boundary crossings obtained from our hybrid simulations and those estimated by Winslow et al. (2013). A negative (positive) time shift corresponds to occurrence of the boundary crossing earlier (later) in the simulations compared to observations. For example, an earlier occurrence when the spacecraft is moving into the magnetosphere means that the magnetospheric boundary has spatially extended further out in our simulations compared to its actual location observed by MESSENGER. We also show a fitted liner regression model into our simulations using an ordinary least-squares method to calculate the Alfvén Mach number and solar wind dynamic pressure at every boundary crossing and to estimate the errors involved in our calculations. The linear assumption, although crude, is valid because of the small ranges of the Alfvén Mach number and solar wind dynamic pressure applied in our simulations.

Figure 2a suggests that the Alfvén Mach number, M_A, during the inbound bow shock crossing is 3.1 ± 0.2, which corresponds to a dynamic pressure of 2.7 ± 0.3 nPa, where P_dyn = /μ₀. On the other hand, Fig. 2b suggests that the solar wind dynamic pressure has slightly increased from ~2.7 nPa during the inbound bow shock crossing to 3.1 ± 0.5 nPa during the inbound magnetopause crossing. Figure 2c also suggests that the trend of increasing solar wind dynamic pressure has continued to the outbound magnetopause crossing where the dynamic pressure is 4.2 ± 0.5 nPa. Duringthe outbound bow shock crossing, however, Fig. 2d suggests that the Alfvén Mach number is 3.6 ± 0.2, which indicates that the dynamic pressure has slightly decreased to 3.6 ± 0.4 nPa. In general, these dynamic pressures are lower than the typical solar wind dynamic pressure near the orbit of Mercury, which is approximately 6–7 nPa.

Comparison between the bow shock (Fig. 2a and d) and magnetopause (Fig. 2b and c) crossings suggest that the location of the bow shock is more sensitive to variations in the solar wind compared to that of the magnetopause location for this orbit. This is mainly because the magnetopause forms much closer to the planet where the magnetic pressure is stronger. Therefore, the magnetopause shows smaller variations in size and location for the limited dynamic pressure ranges we have used in our simulations. Moreover, a comparison between Fig. 2a, b, c, and d indicates smaller statistical variations in the determination of the inbound bow shock crossing compared to other boundaries. As shown in Fig. 3, this is mainly related to the spacecraft geometry and the shape and structure of the magnetospheric boundary during MESSENGER passage throughout the magnetosphere.

Using the WSA-ENLIL solar wind model, Baker et al. (2013) have predicted the solar wind plasma parameters near the orbit of Mercury. Based on their model results (see their Fig. 4), n_sw ≈ 30 cm⁻³ and v_sw ≈ 370 km s⁻¹ for 23 April 2011 (DOY 113), which results in a dynamic pressure of ~7 nPa. This dynamic pressure is nearly two times that inferred from our hybrid simulations using MESSENGER magnetic field observations.

3.1.2 General structure of the magnetosphere

Figure 3 shows a snapshot of our hybrid simulation results from one of the simulations presented in Fig. 2 that showed the best agreement with observations during the outbound bow shock crossing (fordetails see Sect. 3.1.3). In this simulation, the upstream solar wind plasma density n_sw = 22 cm⁻³, the solar wind plasma velocity u_sw = [-310.0, +50.0, 0.0] km s⁻¹, the dynamic pressure P_dyn ≈ 3.6 nPa, the Alfvén Mach number M_A ≈ 3.6, and the plasma β ≈ 0.28. The trajectory of MESSENGER for the D113 orbit is shown in this figure and the UTC time during the passage is marked by arrows in Fig. 3a,b. As shown in Fig. 3a, MESSENGER entered the magnetosphere at ~17:18, passed over the north pole with a closest approach of ~1.14 R_M from the centre of the planet, and moved out of the magnetosphere at ~20:12 (Winslow et al. 2013). The orbit of MESSENGER, as shown in Fig. 3a, made a ~70° angle to the Mercury–Sun line.

Figure 3 shows that the interaction between the supersonic flow of the solar wind plasma and the weak intrinsic magnetic field of Mercury creates an Earth-like magnetosphere containing a bow shock, magnetopause, magnetotail, and funnel-shaped polar cusps. A collisionless bow shock is evident upstream with a large jump in the magnetic field strength and direction shown in Fig. 3b,d, plasma density enhancement shown in Fig. 3e,f, and solar wind velocity deceleration evident in Fig. 3g,h. As shown in Fig. 3c,d, the bow shock forms at x ≈ +2.43 R_M near the subsolar point in the MSM coordinate system, which is extended about 0.5 R_M further upstream compared to the typical distance of the bow shock at Mercury (1.95 R_M Winslow et al. 2013). This is due to the low Alfvén Mach number and low solar wind dynamic pressure used in this simulation compared to the typical Mach number of ~6.6 near the orbit of Mercury (Winslow et al. 2013). However, this distance is much closer to the planet compared to the relative distance of the bow shock at Earth, which is typically at ~15 R_E, where R_E ≈ 6370 km is the radius of the Earth (e.g. Baumjohann & Treumann 1996). The bow shock current intensity is ~70 nA m⁻² at the subsolar point, shown in Fig. 3c, which is also much weaker than that at the Earth, and it generates a weaker magnetic field with a smaller overshoot at the shock (Baumjohann & Treumann 1996; Masters et al. 2013). This is mainly due to the low plasma β and low Alfvén Mach number upstream of Mercury compared to those near the Earth (e.g. Raines et al. 2015, and references therein).

As shown in Fig. 3e–h, when the solar wind passes through the bow shock, its density increases to ~3n_sw while its velocity decreases to ~0.35 |u_sw | at the subsolar point and gets deflected around the magnetosphere, forming the magnetosheath. The magnetosheath gets broader and larger downstream, but its thickness near the subsolar point is ~0.8 R_M for the low dynamic pressure used in this simulation, and is approximately 2 times larger than the typical thickness of the Hermean magnetosheath. The magnetopause, shown in Fig. 3c,d as the innermost intense current near Mercury on the dayside, forms at x ≈ +1.70 R_M in the MSM coordinate system, whereas the average subsolar distance for the magnetopause is ~1.45 R_M (Winslow et al. 2013), which is again due to the low dynamic pressure used in this simulation. Shown in Fig. 3c,d, the magnetopause current flows primarily from dawn to dusk (along the +y axis) with a maximum intensity of ~180 nA m⁻² at the subsolar point with a northward component along the +z axis as it flows around Mercury.

The funnel-shaped polar cusps, which enable the direct access of the solar wind plasma into the surface of Mercury (e.g. Killen et al. 2001; Massetti et al. 2003; Winslow et al. 2012), can be seen in Fig. 3f over the poles near the surface. We see from Fig. 3f that the plasma density over the northern cusp (~70 cm⁻³) is larger than that over the southern cusp (~55 cm⁻³). This asymmetry in the solar wind plasma density between the southern and northern cusps is mainly associated with the northward displacement of the planetary magnetic field. The total magnetic field over the north pole of Mercury is stronger than that over the south pole. This allows an easier access of the solar wind plasma to the southern hemisphere on the dayside compared to the northern hemisphere, and results in plasma density increases over the northern cusp compared to the southern cusp. In addition to the cusps, there are two other notable features including the current sheet in the magnetotail, which can be seen in Fig. 3b and d on the nightside close to the magnetic equator (z ≈ +0.2R_M) and the signature of quasi-trapped solar wind protons on the nightside, evident in Fig. 3f close to the planet between x = −1.0 R_M and x = −1.7 R_M. Our simulations show that the current sheet has a plasma density of 10–20% of the upstream solar wind and the plasma velocity reaches over 150% of the solar wind velocity because of magnetic reconnection in the magnetotail. This reconnection causes the ions to move away or towards Mercury from the reconnection site (X-line). Those that move towards Mercury, if not moving to the dayside magnetosphere, have direct access to Mercury’s high latitudes on the nightside mainly along the open field lines. Since Mercury’s magnetosphere is highly dynamic, a time-series of our simulations (not shown here) suggests that the trapped particles at Mercury’s nightside is not a permanent feature of the magnetosphere (e.g. Luhmann et al. 1998), thus we call them the quasi-trapped plasma. Our simulations presented in Fig. 3e,f show that the density of the quasi-trapped particles on the nightside is comparable to the upstream solar wind density, and that their energy can reach ~10–20% higher than the upstream solar wind energy as they bounce between the two magnetic mirror points and drift duskward.

3.1.3 Comparison with MESSENGER observation

Figure 4 shows a comparison between the magnetic field obtained from our hybrid simulation presented in Fig. 3 (red lines), MESSENGER magnetometer observations (black lines), and the undisturbed magnetic dipole of Mercury (dashed lines) along the trajectory of MESSENGER for orbit D113. Figure 4a–d show that there is a good agreement in the overall trend of the magnetic fields between our hybrid simulation and MESSENGER observations. In addition, we see that the location of the magnetospheric boundaries have been estimated correctly in our model, especially for the outbound bow shock crossing. The location of magnetospheric boundary crossings can also be seen in the electric current density calculated from our simulations and shown in Fig. 4e. As we showed previously in Fig. 2, we estimated the dynamic pressure during the inbound crossing to be lower than that used in the simulation results presented in Figs. 3 and 4. Therefore, the bow shock and magnetopause boundaries shown in Fig. 4 have moved slightly closer to the planet during the inbound part of the orbit compared to those observed by MESSENGER. However, the general trend of the magnetic fields is similar to the observations.

In addition to the magnetic fields, the plasma parameters shown in Fig. 4f also reveal some of the general characteristics of magnetospheric plasma along the trajectory of MESSENGER. These include sharp jumps in plasma density and velocity at the bow shock (~17:18 at inbound and ~20:12 at outbound) and plasma density enhancement by nearly a factor of 2 and velocity reduction by nearly 75% in the magnetosheath between ~17:18 and ~17:35 during the inbound and between ~18:52 and ~20:12 during the outbound. Another notable feature is the plasma density enhancement near 17:50, which is an indication of a northern cusp crossing, also evident in Fig. 3f. Our model suggests that due to the geometry of the orbit, the quasi-trapped particles could not have been observed in this orbit.

In Fig. 5 we present another example of our model–data comparison for a lower dynamic pressure than that presented inFig. 4. As shown in Fig. 5, there is a better agreement between our simulations and magnetic field observations during the inbound crossing, but not at the outbound. In this simulation, the upstream solar wind density n_sw = 18 cm⁻³ and the solar wind velocity is u_sw = [–300.0, +50.0, 0.0] km s⁻¹, resulting in a dynamic pressure P_dyn ≈ 2.7 nPa and an Alfvén Mach number M_A ≈ 3.1. We see from Fig. 5d,e that the location of the bow shock and magnetopause agree fairly well with MESSENGER observations during the inbound crossing, but not at the outbound. As shown earlier in Fig. 2, our model estimated that the solar wind dynamic pressure and Alfvén Mach number were lower during the inbound than at the outbound. Thus, as Fig. 5d,e show, a solar wind dynamic pressure higher than 2.7 nPa is required to push the magnetopause and bow shock closer to the planet during the outbound magnetospheric crossing.

Richer et al. (2012) have also compared the magnetic field observations along the D113 orbit with a three-dimensional hybrid model of plasma. They examined the effects of two mathematical descriptions for the intrinsic magnetic field of Mercury including (1) a single dipole with northward displacement with similar parameters to those used in our simulations, and (2) a single dipole combined with a quadrupole fitted to the northward displaced dipole fields. They concluded that despite the similarities in the structure of the magnetic fields observed over the north pole, the topology of the fields over the south pole present considerable differences between the two models. However, their simulations for a single dipole with northward displacement presents notable differences compared to our simulations shown in Fig. 4. Richer et al. (2012) used a solar wind plasma density of 32 cm⁻³ and a bulk flow speed of 430 km s⁻¹, which results in a dynamic pressure ~10 nPa. As presented in Fig. 2, this dynamic pressure is ~3 times higher than the value we have estimated for the upstream solar wind pressure using our simulations. Nevertheless, there is a good agreement for the location of the inbound bow shock crossing between MESSENGER observations and the Richer et al. (2012) single dipole model. Conversely, as shown in Fig. 3a by Richer et al. (2012), the outbound bow shock crossing is located further upstream (time difference is nearly 15 min), which suggests that an even higher dynamic pressure is required in Richer’s model to capture the location of the bow shock in the same place as observed by MESSENGER. Moreover, the magnetic field strength from Richer et al. (2012) simulations near the closest approach is much higher than that observed by MESSENGER, which is perhaps related to the large dynamic pressure they have used in their simulations. These disagreements between the Richer et al. (2012) simulations and the MESSENGER observations may also suggest that the dynamic pressure used in their simulations, which is close to those estimated by Baker et al. (2013) from the WSA-ENLIL model for D113, is higher than the actual solar wind dynamic pressure during this orbit. This perhaps confirms the dynamic pressure estimated by our simulations (~3.5 nPa) for the D113 orbit.

Fig. 2 Time difference between our hybrid simulations and the mid-point location of the magnetospheric boundary crossings obtained from MESSENGER magnetic field observations by Winslow et al. (2013). The dashed line is a linear fit of our results using an ordinary least-squares fitting function. Panels a and b: inbound bow shock and magnetopause crossings. Panels c and d: Outbound magnetopause and bow shock crossings for the second magnetospheric transient on 23 April 2011 (D113). Pdyn = is the solar wind dynamic pressure and MA = is the Alfvén Mach number.

Fig. 3 Hybrid simulation results for the magnetospheric transit of MESSENGER on 23 April 2011 between 16:00 and 21:00 (D113) presented in the MSO coordinate system. Panels a and b: magnitude of the magnetic field in logarithmic scale. Panels c and d: electric current density calculated from the general Ampère’s law, flow normal to the presented planes. Panels e and f: plasma density in logarithmic scale and normalized to the upstream solar wind density, nsw = 22 cm−3. Panels g and h: magnitude of the solar wind velocity normalized to the upstream solar wind velocity, |usw | = 314 km s−1. The top panels are cuts in the equatorial plane (xy-plane at z = 0), viewed from Mercury’s north pole, and the bottom panels are cuts in the midnight meridian plane (xz-plane at y = 0), viewed from the orbital motion of Mercury (i.e. the -y axis). Mercury is shown by a circle, centred at the origin of the coordinate system. The direction of the solar wind and the IMF are shown by yellow and white arrows, respectively, in panels a and b. Streamlines in panels c and d show magnetic field line tracing, and in panels g and h show the plasma flow direction. A portion of MESSENGER’s orbit on 23 April 2011 (D113), also shown in Fig. 1, is shown in all panels and the UTC time during the passage is marked by arrows in panels a and b.

Fig. 4 Panels a–d: magnetic field comparison between our hybrid model simulations (red lines), MESSENGER magnetometer observations (black lines), and undisturbed intrinsic magnetic dipole of Mercury (blue dashed lines) along the trajectory of MESSENGER on 23 April 2011 (D113) between 16:00 and 21:00 UTC. Panel e: magnitude of the electric current density calculated from our simulations using the general Ampère’slaw. Panel f: solar wind bulk flow speed normalized to the upstream solar wind speed |usw | = 314 km s−1 (purple line) and solar wind plasma density normalized to the upstream plasma density nsw = 22 cm−3 obtained from our hybrid model simulations along the trajectory of MESSENGER. The mid-point location of the bow shock (BS) and magnetopause (MP) boundaries estimated by Winslow et al. (2013) as well as the closest approach (CA) to the planet are shown by the vertical lines.

Fig. 5 Hybrid simulation comparison with MESSENGER observations for the upstream solar wind dynamic pressure ~2.7 nPa, similar to our model estimated dynamic pressure during the inbound magnetospheric crossing shown in Fig. 2. The figure format is the same as that shown in Fig. 4.

Fig. 6 Time difference between our hybrid simulations and the location of the magnetospheric boundary crossings during the M1 flyby obtained from MESSENGER magnetic field observations by Slavin et al. (2008). The dashed lineis a linear fit of our results using an ordinary least-squares fitting function. Panels a and b: inbound bow shockand magnetopause crossings. Panels c and d: outbound magnetopause and bow shock crossings. Pdyn = is the solar wind dynamic pressure and MA = is the Alfvén Mach number.

3.2 M1 flyby: 14 January 2008

As shown in Fig. 1, during the M1 flyby MESSENGER entered Mercury’s magnetosphere from the magnetotail nearly at x = –4 R_M, crossed the plasma sheet, moved to the nightside of Mercury with a closest approach at ~200 km altitude from the surface, and moved out of the magnetosphere after probing the magnetosheath and bow shock on the dayside (Slavin et al. 2008; Anderson et al. 2011a). In order to estimate the solar wind dynamic pressure and Alfvén Mach number upstream of Mercury during this flyby we performed over 40 simulation runs for the upstream solar wind plasma density ranging from 26–38 cm⁻³ with a span of 2 cm⁻³, and with the velocity along the -x axis ranging between 340 and 400 km s⁻¹ with a span of 10 km s⁻¹ and a constant speed of 50 km s⁻¹ along the +y axis. These parameters cover a dynamic pressure range between ~5.0 and ~10.2 nPa, an Alfvén Mach number between ~4.3 and ~6.2, and a plasma β between ~0.34 and ~0.50.

Slavin et al. (2008), using magnetic field observations, have estimated the times that MESSENGER passed the different magnetospheric boundaries during the M1 flyby. Figure 6 shows the time difference between our simulations and those observations. Figure 6b and c suggest that the solar wind dynamic pressure did not change considerably from the inbound to the outbound magnetopause crossing, and that it was 7.2 ± 1.1 nPa. The Alfvén Mach number, however, showed a slight variation from 5.2 ± 0.4 during the inbound bow shock crossing (Fig. 6a) to 5.8 ± 0.8 during the outbound bow shock crossing (Fig. 6d). These values correspond to a solar wind dynamic pressure of ~7.3 nPa for the inbound and ~9.0 nPa for the outbound bow shock crossing.

In general, our simulations suggest that the solar wind remained relatively steady during the M1 flyby, and that its dynamic pressure was ~7 nPa, which is close to the typical solar wind dynamic pressure near the orbit of Mercury. However, as shown in Fig. 6b-d, our model was not able to readily determine the location of the inbound and outbound magnetopause and the outbound bow shock crossing, and a relatively large error is involved in our estimations compared to those shown in Fig. 2. As shown later in Fig. 7, this is mainly because of the geometry of the M1 flyby that made it difficult to estimate the location of the boundaries from our simulations.

Figure 7 shows the global structure of the solar wind plasma interaction with Mercury during the M1 flyby from one of our hybrid simulations presented in Fig. 6. In this simulation n_sw = 32 cm⁻³ and u_sw = [–360.0, +50.0, 0.0] km s⁻¹, leading to a subsolar dynamic pressure of ~7.1 nPa, an Alfvén Mach number ~5.1 and a plasma β ≈ 0.42, which are close to the typical solar wind plasma parameters near the orbit of Mercury. The general characteristics of the magnetosphere (see Fig. 3 and Sect. 3.1.2) are again evident in Fig. 7. Since the solar wind dynamic pressure is higher in this simulation compared to that shown in Fig. 3, the bow shock and magnetopause boundaries at the subsolar point have been pushed closer to Mercury by the solar wind. As shown in Fig. 7c,d, the subsolar distance of the bow shock and magnetopause in the MSM coordinate system are at 1.70 R_M and 1.45 R_M, respectively.While the magnetopause location is the same as the average magnetopause distance obtained from MESSENGER observations (Winslow et al. 2013), the bow shock stands ~0.2 R_M closer to the planet compared to the average distance of the bow shock (~1.95 R_M as estimated by Winslow et al. 2013). A comparison between Figs. 3d and 7d also shows that the magnetopause current density at the subsolar point is higher for the M1 flyby (~300 nA m⁻² in Fig. 7d, but ~180 nA m⁻² in Fig. 3d), which is an indication of the higher dynamic pressure applied in the simulation presented in Fig. 7 compared to that applied in Fig. 3.

A notable feature evident from Fig. 7e is a partial ring current at Mercury’s nightside near the equatorial plane. Consistent with MESSENGER observations (e.g. Schriver et al. 2011; Korth et al. 2014), our simulations show that this half-ring contains quasi-trapped particles only near the nightside equator that move duskward toward the dayside magnetopause. The quasi-trapped particles, also visible in Fig. 7f in the nightside near Mercury, are either lost by impacting the surface or by crossing the magnetopause into the solar wind.

Since the direction of the IMF is mainly parallel to the solar wind flow (shown by the arrows in Fig. 7a,b), a quasi-parallel shock region forms near the subsolar point, allowing the reflected solar wind ions to flow upstream and form a foreshock that disturbs the upstream solar wind plasma and fields. Fig. 7b,h show signatures of magnetic field perturbations and solar wind plasma velocity reductions in the foreshock region (+1 R_M ≤ x ≤ +4 R_M and − 3 R_M ≤ z ≤−1 R_M). However, wecould not find clear evidence of plasma density reductions there as an indication of a foreshock cavity (e.g. Schwartz et al. 2006; Blanco-Cano et al. 2009).

Moreover, as shown in Fig. 7c,d, the quasi-parallel shock near the subsolar point and the large solar wind dynamic pressure upstream perturb magnetospheric boundaries such that distinguishing the electric currents and determining the magnetospheric boundaries upstream are not an easy task. In addition, the trajectory of MESSENGER during the M1 flyby nearly lies in the equatorial plane of Mercury where it is difficult to calculate the magnetopause currents, as shown in Fig. 7c. These are the main reasons that we could not clearly determine the location of the magnetospheric boundaries, as shown in Fig. 6. However, the foreshock, and thus its associated phenomena, is a kinetic process that can only be explained by kinetic models like our hybrid plasma model used in this study, and not by MHD models.

Figure 8 compares the magnetic field and solar wind plasma density and velocity obtained from our hybrid simulation presented inFig. 7 with the magnetic fields observed by MAG and plasma counts observed by FIPS on MESSENGER. Figure 8a–d show that there is generally an agreement between our simulation results (red line) and MESSENGER magnetic field observations along the trajectory of MESSENGER during the M1 flyby. They also show that the solar wind plasma interaction with the intrinsic magnetic field of Mercury creates features that are distinct from pure dipole magnetic fields (dashed blue line). For example, there is a jump in the intensity of the magnetic fields at the inbound bow shock crossing, evident in Fig. 8a and d near 18:08, and the deflection of the magnetic field in the outbound magnetosheath, shown in Fig. 8a,b between 19:14 (magnetopause) and 19:19 (bow shock). We see from Fig. 8e that the location of the magnetospheric boundaries have been correctly estimated by our model, especially for the inbound bow shock crossing. As shown in Fig. 6, the upstream solar wind dynamic pressure remained relatively steady during the M1 flyby. Thus, the location of all magnetospheric boundaries have been estimated correctly by our model, as also shown in Fig. 8e.

During the M1 flyby, MESSENGER crossed the plasma sheet and provided estimates of the plasma environment at the sheet with proton density 1–10 cm⁻³, proton temperature ~170 eV, and a steady plasma beta ~2 (Raines et al. 2011). Figure 8f from our simulation along the trajectory of the M1 flyby shows that the solar wind plasma density was 5–30 cm⁻³ between the inbound magnetopause crossing, i.e. ~18:43, and the closest approach ~19:05. Plasma temperature in the plasma sheet estimated from the velocity space distribution of the solar wind protons in our simulation (not shown here) is ~130 ± 70 eV. These results are within the range obtained from the MESSENGER plasma observations by Raines et al. (2011). Moreover, a comparison between simulation results presented in Fig. 8f and the number of counts measured by the FIPS on MESSENGER (Fig. 8g taken from Raines et al. 2011) shows a fair agreement between simulations and plasma observations after the inbound magnetopause crossing (~18:43). Since we have not accounted for the FIPS field of view and coverage (~1.4π) to calculate plasma parameters in our simulation, and since the FIPS field of view is limited to directions transverse to the Mercury–Sun line (Raines et al. 2011), we cannot comment on disagreements between our simulations and FIPS plasma observationsbefore the inbound magnetopause crossing. However, since MESSENGER has passed the magnetosheath between ~18.08 and ~18:43, we expect the plasma flux to be higher in that region, as shown by our simulations in Fig. 8f, compared to that between the inbound magnetopause and the closest approach.

Anothernotable feature observed during the M1 flyby is the observation of the dayside double boundary layer (DBL) near the outbound magnetopause, also known as the double current sheet and double magnetopause (Slavin et al. 2008; Anderson et al. 2011b; Müller et al. 2012). Different scenarios have been suggested for the formation of the DBL at Mercury, including the plasma pressure gradient at the inner magnetosphere, perhaps because of the quasi-trapped particles on closed field lines (Müller et al. 2012) and/or the existence of exospheric plasma (Slavin et al. 2008), and plasma injection through the cusps or magnetosphere flanks (Liljeblad et al. 2015). The DBL at Mercury is mainly occurring on the dawnside (Liljeblad et al. 2015) and in contrast to that of the Earth shows either little or no dependence on the IMF direction (Anderson et al. 2011b; Liljeblad et al. 2015). The DBL during the M1 flyby is evident in the MESSENGER magnetic field observations shown in Fig. 8c,d as two consecutive large magnetic field drops, one near 19:10 and one near 19:14. As explained by Slavin et al. (2008), since the outer boundary layer (the one observed around 19:14) separates the magnetosheath and magnetosphere, we consider that to be the magnetopause boundary crossing. However, in contrast to MESSENGER observations (Slavin et al. 2008) and previous hybrid simulations for the M1 flyby (Müller et al. 2012), our model does not show evidence of the DBL, neither in the magnetic field nor in the total electric current density shown in Fig. 8e. This is perhaps due to the larger simulation cell size in our simulations (200 km) compared to those chosen by Müller et al. (2012) (adaptive grids with cell sizes smaller than 80 km near Mercury), or the lack of enough plasma to form a large pressure gradient to form the DBL.

However, it is worth noting that Müller et al. (2012) did not consider the northward displacement of Mercury’s intrinsic magnetic field or the tenuous exosphere of Mercury in their simulation. Moreover, despite the large x-component of the IMF upstream (see Table 1), Müller et al. (2012) assumed that the IMF is only along the y-axis, whereas we assumed that the y component of the IMF is zero in our simulations. However, since Müller et al. (2012), consistent with the MESSENGER observations, have observed the DBL in their simulations, and since MESSENGER magnetic field observations show a large component for the IMF along the y-axis after the outbound bow shock crossing (Fig. 8b between 19:20 and 19:45), it raises the necessity of revisiting the effects of the IMF direction in the formation of the DBL using simulations, which is outside the scope of this study.

3.3 Orbit D182: 01 July 2011

As shown in Fig. 1, the orbit of MESSENGER on 01 July 2011 (D182) is almost in the midnight meridian plane, and moves from the dayside magnetosphere into the deep midnight of Mercury after probing the northern cusp. Our analysis of the MESSENGER magnetic field observations (see Table 1) indicate that the strength and direction of the IMF one hour before (after) MESSENGER moved into (out of) the magnetosphere of Mercury have changed considerably. One hour before the inbound, the IMF had a large component along the -y axis and it was nearly southward, whereas one hour after the outbound it was mainly in the xz-plane with a northward component. In addition, the IMF was ~1.2 times stronger during the outbound compared to the inbound, which affects the solar wind Mach number and plasma β during the inbound and outbound crossings.

Similar to the previously analysed orbits, we used the inverse problem approach to infer the solar wind dynamic pressure from our hybrid simulations through comparison with MESSENGER magnetic field observations. We made two sets of simulations: one for the inbound IMF and one for the outbound, and for each set we ran nearly 15 simulation runs. Despite changes in the IMF direction and strength, our simulations (not shown here) suggest that the solar wind dynamic pressure had ~1.5 nPa variations during the D182 orbit, but onaverage it was ~7.5 nPa. Simulation results from the WSA-ENLIL model predicts that the solar wind dynamic pressure during the D182 orbit is ~10.5 nPa (see Fig. 4 in Baker et al. 2013), which is nearly 30% higher than that estimated from our model.

Figure 9 compares the global structure of the Hermean magnetosphere during the D182 orbit before the inbound (top panels) and after the outbound (bottom panels) magnetospheric crossings in the midnight meridian plane (almost the orbital plane of MESSENGER). As shown in Table 1 and by the arrows in Fig. 9a,b, the IMF direction changed from southward before the inbound to northward after the outbound bow shock crossing. The results presented inFig. 9 are obtained from two simulation runs where the solar wind plasma density n_sw = 30 cm⁻³ and velocity u_sw = [–380.0, +50.0, 0.0] km s⁻¹ are the same for both of the runs, but the IMF orientation and strength are different, as listed in Table 1. Therefore, the solar wind dynamic pressure for both of the simulations was ~7.4 nPa, but the Alfvén Mach number changed from M_A = 5.1 during the inbound to M_A = 4.3 during the outbound. Similarly, plasma β = 0.38 during the inbound changed to β = 0.27 during the outbound.

In general, Fig. 9 shows that there are considerable differences in the structure of the Hermean magnetosphere during the inbound and the outbound magnetospheric crossings. Comparison between Fig. 9a and b show that during a northward IMF, the magnetic field strength gets weaker at the south pole because of the northward displacement of Mercury’s magnetic dipole. As shown in Fig. 9b, more opposing fields form around the southern cusp during the outbound magnetosphere crossing, which may result in the formation of a quadrapole field structure and magnetic null-points around the southern cusp. In addition, the foreshock disturbances, evident in Fig. 9b, also perturb the magnetic fields in the southern magnetosphere in the midnight meridian plane during the northward IMF. However, since the IMF during the inbound has a considerable component along the y-axis, the foreshock cannot be seen in the midnight meridian cut in Fig. 9a, but instead it forms on a plane perpendicular to those shownin this figure, as discussed and shown later in Fig. 10.

Comparison between the electric currents shown in Fig. 9c and d indicates that the dayside magnetopause current during a southward IMF (290 nA m⁻² at the subsolar point) is stronger than that for a northward IMF (200 nA m⁻² at the subsolar point). This is in agreement with previous observations (e.g. Slavin et al. 2009; Anderson et al. 2011a). Moreover, we see that the magnetotail current lies further downstream of the planet (x ≤ –2 R_M) and is slightly shifted northward during the southward IMF, whereas for the northward IMF it lies closer to the planet (x ≤−1.5 R_M) and forms near the magnetic equatorial plane.

Comparison between Fig. 9e and f shows that the quasi-trapped particles at Mercury’s midnight meridian are observable during a northward IMF. It is classically known that a southward IMF is a more favourable orientation for magnetic reconnection on the dayside magnetopause of Mercury (Slavin et al. 2008, 2009). Magnetic reconnection, which is a fundamental plasma process, converts magnetic energy into particle energy and leads to the transport of energy and momentum into planetary magnetospheres (Paschmann et al. 2013). However, our simulations suggest that the quasi-trapped particles in the magnetosphere of Mercury are mainly entering the magnetosphere from magnetospheric flanks near the equatorial plane (not shown here) and the magnetic reconnection at the dayside is not the driver for their entrance into the nightside magnetosphere. This is mainly associated with the direction of the magnetic field that provides easier access of the solar wind particles into the magnetosphere of Mercury through drift motions during the northward IMF compared to the southward IMF. These particles become subject to gradient and curvature drift motions and form a partial ring structure near the equator at the nightside (Baumjohann & Treumann 1996; Fatemi & Poppe 2018). Figure 9e,f also show that the polar cusps have a broader loss cone and move more towards mid-latitudes during the northward IMF compared to that for the southward IMF, which is in agreement with previous observations (e.g. Winslow et al. 2014).

Figure 10 shows a cross section of Mercury’s magnetosphere on the dayside at x = +1.0 R_M (first two columns) and nightside at x = –1.5 R_M (last two columns) for the inbound (top panels) and outbound (bottom panels) IMF configurations during the D182 orbit. Generally speaking, there is a north–south asymmetry in the structure of the magnetosphere mainly associated with the northward displacement of the planetary dipole field and the orientation of the IMF. In addition, there is also a dawn–dusk asymmetry (in the MSO coordinate system) that is due to the dawnward component of the solar wind plasma velocity applied in our simulations. However, a comparison between the top panels and bottom panels in Fig. 10 indicates that the IMF orientation has a major contribution in the north–south and dawn–dusk asymmetries in Mercury’s magnetospheric structure. This is in agreement with previous observations (e.g. Slavin et al. 2010; DiBraccio et al. 2013; Korth et al. 2014), and shows that the IMF orientation strongly controls the dynamics of the solar wind plasma interaction with Mercury.

Figure 10e,f show that the foreshock cavity forms at close distances to Mercury near the south pole during the northward IMF and that it extends to over 6 R_M upstream where the size of the perturbed area is comparable to or larger than the size of the planet. During the southward IMF, however, since the magnetic field has a large component along the -y axis, the foreshock cavity does not seem to affect the structure of the magnetosphere at close distances to Mercury much, as shown in Fig. 10a,b. In addition to the foreshock, trailing waves associated with the bow shock that propagate outward can also be seen upstream of Mercury in Fig. 10a,e. As we explained before and is evident from Fig. 10d, there is no signature of quasi-trapped particles close to Mercury on the nightside during the southward IMF. In contrast, as shown in Fig. 10h, the quasi-trapped particles during the northward IMF has a flux of 20–40% of the upstream solar wind flux for the solar wind and IMF parameters chosen in this simulation.

Figure 11a–d compare MESSENGER magnetic fields observed during the D182 orbit (black lines) with our hybrid simulations, presented in Figs. 9 and 10, for the IMF conditions during the inbound (red lines) and outbound (green lines) crossings. To provide better visualization, we have subtracted the undisturbed magnetic dipole fields from observations and from simulations. We see that there is a good agreement between MESSENGER magnetic field observations and our hybrid simulations if we assume that there was a change in the IMF orientation and magnitude before or near the closest approach. We see from Fig. 11b that nearly 60 nT changes in the B_y on the dayside magnetosheath between ~06:55 and ~07:05 has been well reproduced by the inbound simulation (red line). This occurs because the IMF had a large y-component during the inbound that gets diverted around the dayside magnetosphere in the magnetosheath in the plane of MESSENGER’s orbit. In addition, we see from Fig. 11a,c that the inbound magnetopause crossing (~07:05) and its associated magnetic field changes are explained fairly well with the inbound simulation. Fig. 11e also shows the agreement in determination of the inbound magnetospheric boundary crossings between MESSENGER observation and our inbound simulation.

A comparison between the inbound and outbound simulations in Fig. 11e also suggests that the magnetopause current layer has moved closer to the planet during the northward IMF (outbound) compared to the southward IMF (inbound). We can also see this by comparing the location of the magnetopause currents in Fig. 9c,d and the size and structure of the dayside magnetospheric cavity in Fig. 9e,f. This is in contrast to previous works that have suggested the dayside magnetopause at Mercury lies further upstream during a northward IMF than during a southward IMF (e.g. Slavin & Holzer 1979; Ip & Kopp 2002; Kidder et al. 2008; Trávníček et al. 2010). There are at least three possible scenarios to explain this disagreement: (1) there are some variations in the IMF orientation and strength between the inboundand outbound simulations we have presented here. There is no change in the solar wind dynamic pressure between the two simulation runs, but we have shown in Fig. 10 that the dynamic magnetosphere of Mercury is highly controlled by the IMF, which may affect the size and structure of the dayside magnetosphere; (2) as the orientation of the IMF changes, the magnetic reconnection sites also change correspondingly, which highly affects the transport of the solar wind plasma into the magnetosphere and alters the dynamics and structure of the magnetosphere; and (3) in all the above-mentioned studies that discuss the location of the dayside magnetopause and its response to the IMF direction, the northward displacement of the dipole moment has not been considered. These three plausible reasons need to be investigated and examined in detail in a separate study.

We showed in Fig. 9f that there are quasi-trapped particles in Mercury’s nightside when the IMF has a northward component. The signature of these particles has been investigated before using simulations (e.g. Trávníček et al. 2010; Herčík et al. 2016). They have also been observed by MESSENGER FIPS instrument at the nightside of Mercury during the D182 orbit (Herčík et al. 2016). In agreement with previous simulations and observations, Fig. 11f also shows the detection of the trapped particles along the D182 orbit near the plasma sheet between ~07:30 and ~07:40 only for the outbound simulation where the IMF is northward. The density of the trapped particles is comparable to the upstream solar wind density, while their velocities are ~15% of the upstream solar wind velocity, shown in Fig. 11f,g. Possibly, a slight reduction in the magnitude of the magnetic fields observed between ~07:30 and ~07:40, shown in Fig. 11d, is either related to the existence of trapped particles or is an indication of a current sheet crossing, or both, observed only in the outbound hybrid simulation.

As shown in Fig. 11a–d, during the outbound magnetosheath crossing (between ~08:26 and ~09:08), we see high-frequency large amplitude waves in MESSENGER magnetic field observations, which makes the identification of the boundary crossings difficult. We do not observe these waves in our simulations, mainly because our simulation cell sizes are too large to capture such wave activities. However, the boundary crossings are more evident in our simulations than in the observations. The electric current and a large change in plasma density around 08:50 in our outbound simulation, shown in Fig. 11e,f, suggest that the outbound bow shock crossing occurred nearly 15 min earlier than suggested by Winslow et al. (2013) from MESSENGER magnetic field observations. This shows that our model not only estimates the solar wind plasma parameters upstream of Mercury, but may also help to better determine magnetospheric boundary crossings.

Herčík et al. (2016) have also compared the results of their hybrid simulations with MESSENGER observations for the D182 orbit. The solar wind dynamic pressure in their simulation is ~5 nPa, which is slightly smaller than the dynamic pressure we have estimated from our simulations. However, the location of the magnetospheric boundaries from their simulations, showed in Fig. 10b by Herčík et al. (2016), is in good agreement with those estimated from MESSENGER magnetic field observations, but the magnitude of the magnetic field from their simulations not only does not match MESSENGER observations during the closest approach, but also shows a non-dipolar structure close to Mercury. Asstated by Herčík et al. (2016), the actual size of Mercury in their simulations is half of its real size (~1200 km), which may affect the physics of the interaction and the shape and structure of the magnetosphere.

Fig. 7 Global structure of the Hermean magnetosphere during the M1 flyby obtained from our hybrid simulations in the same formatas Fig. 3.

Fig. 8 Comparison between our hybrid simulations and MESSENGER observations along the trajectory of MESSENGER on 14 January 2008 during the M1 flyby. The format of the presented results is the same as that shown in Fig. 4. Panel f: upstream solar wind speed |usw | =365 km s−1 and solar wind plasma density nsw = 32 cm−3. Panel g: plasma counts observed during the M1 flyby by the FIPS instrument on MESSENGER, digitized from Fig. 3 in Raines et al. (2011).

Fig. 9 Hybrid simulation results for the magnetospheric transit of MESSENGER on 01 July 2011 (D182) for the IMF conditions (top panels) before the inbound and (bottom panels) after the outbound bow shock crossings. Panels a and b: magnitude of the magnetic field in logarithmic scale. Panels c and d: electric current densities calculated from the general Ampère’s law moving along the y-axis, Jy, perpendicular to the plane shown here. Panels e and f: plasma density in logarithmic scale and normalized to the upstream solar wind density, nsw = 30 cm−3. Panels g and h: magnitude of the solar wind velocity normalized to the upstream solar wind velocity, |usw | ≈ 385 km s−1. All the panels are cuts in the midnight meridian plane (xz-plane at y = 0), viewed from the orbital motion of Mercury (i.e. the -y axis). The rest of the figure format is the same as that shown in Fig. 3.

Fig. 10 Hybrid simulation results for the D182 orbit in the planes perpendicular to those shown in Fig. 9. Panels a–d: for the inbound and panels e–h: for the outbound IMF conditions. The cuts are in the yz-plane panels a,b,e and f: at x = +1.0RM, and panels c,d,g and h: at x = –1.5 RM. Panels a,c,d and e: magnitude of the magnetic field in logarithmic scale. Panels b,d,f and h: plasma flux normalized to the upstream solar wind flux Γsw = 1.15 × 1013 m−2⋅s. The solar wind flows into the planes shown here with 50 km s−1 along the +y axis and the direction of the IMF in the yz-plane are shown by black arrows in panels a and e.

3.4 Response of the magnetopause to the solar wind dynamic pressure

Winslow et al. (2013) used MESSENGER magnetic field observations to estimate the stand-off distance of the magnetopause as a function of the solar wind dynamic pressure predicted by the WSA-ENLIL model by Baker et al. (2013). They found that the shape of the magnetopause is well described by the Shue et al. (1997) model, given by , where R is the distance of the magnetopause from the dipole centre, R_ss is the subsolar stand-off distance of the magnetopause from the dipole centre, θ is the angle between R and the dipole-Sun line, and α is a flaring parameter (Winslow et al. 2013). They found that R_ss = 1.45 R_M and α = 0.5 provides the best fit to the midpoint magnetopause crossings observed by MESSENGER, yields to a power law function with a sixth-root relationship between the stand-off distance, R_ss, and the solar wind dynamic pressure, P_dyn, given by Winslow et al. (2013): (1)

In Fig. 12 we show R_ss as a function of P_dyn obtained from Eq. (1) (blue line). We also show the lower (light blue) and upper (dark blue) bounds for that function, and the acceptance area (in yellow). Also shown is the response of the magnetopause to large solar wind dynamic pressure changes during strong coronal mass ejections (CMEs) and high-speed stream (HSS) events observed by MESSENGER (green asterisks; Slavin et al. 2014), and the green dashed line shows the function best fitted to the extreme events, taken from Fig. 19 by Slavin et al. (2014). The solar wind dynamic pressure for the extreme solar events was estimated from the pressure balance between the solar wind plasma and the magnetic pressure at the location of the magnetopause, assuming that the plasma pressure inside the magnetopause is negligible (Slavin et al. 2014).

Since for each of the presented orbits we made tens of simulation runs (in total over 100 simulations) for different solar wind dynamic pressures, we compiled all of these simulations to estimate the response of the magnetopause to the solar wind dynamic pressure from our model results. The red circles in Fig. 12 show our simulation results where we have estimated the R_ss at the subsolar point in the MSM coordinate system (i.e. at y = 0 and z = 0) as a function of the upstream solar wind dynamic pressure, P_dyn = . Similar to analyses by Winslow et al. (2013) and Slavin et al. (2014), we have also provided a sixth-root least-squares fit to our simulation results, shown by the solid red line, which yields a function given by (2)

A comparison between our model-predicted R_ss as a function of P_dyn and those estimated by Winslow et al. (2013) from MESSENGER observations combined with the WSA-ENLIL solar wind prediction model, we find that our simulation results agree fairly well with the lower limit of the Winslow’s empirical model for dynamic pressures below ~10 nPa and agrees fairly well with the Slavin et al. (2014) observations during the extreme solar events. However, since we do not have any observational points from our simulations for solar wind dynamic pressures larger than 10 nPa in this study, our fitted curve should be interpreted with care.

Fig. 11 Comparison between our hybrid simulations and MESSENGER observations along the trajectory of MESSENGER on 01 July 2011 between 06:00 and 10:00 UTC (D182). Panels a-d: magnetic field comparisonbetween MESSENGER observations (black lines) and our hybrid simulations for the inbound (red lines) and outbound (green lines) IMF conditions listed in Table 1 along the trajectory of MESSENGER. To provide better visualization, we have subtracted the intrinsic magnetic fields of Mercury, BD, from observations and from simulations. Panel e: electric current density. Panel f: plasma density normalized to the upstream solar wind plasma density nsw = 30 cm−3. Panel g: plasma velocity normalized to the upstream solar wind plasma speed |usw | ≈385 km s−1, obtained from our simulations along the trajectory of MESSENGER. The mid-point location of the bow shock (BS) and magnetopause (MP) boundaries estimated by Winslow et al. (2013) as well as the closest approach (CA) to the planet are shown by the vertical lines.

Fig. 12 Stand-off distance of the magnetopause, Rss, from the centre of Mercury at the subsolar point in the MSM coordinate system as a function of the solar wind dynamic pressure, Pdyn. The solid blue line is from a MESSENGER empirical model by Winslow et al. (2013) shown in Eq. (1), where the lower limit and upper limit of the model are marked by the light blue and dark blue dashed line, respectively. The acceptance area between the lower and upper limits is in yellow. Taken from Fig. 19 by Slavin et al. (2014), the green asterisks are MESSENGER observations obtained during three extreme solar wind dynamic pressure events (i.e. CMEs and HSS) and the green dashed line is the function fitted to the observations by Slavin et al. (2014). The red circles are from our hybrid simulations and the red solid line is the curve fitted to our simulations, shown in Eq. (2), using a least-squares curve fitting function.

4 Discussion

The lack of an upstream solar wind plasma monitor when a spacecraft is inside the highly dynamic magnetosphere of Mercury limits the interpretations of observed magnetospheric phenomena and their correlations with upstream solar wind variations. An accurate understanding of the solar wind plasma and its variations as it interacts with Mercury is crucial to better understand the morphology of the interaction, the structure of the magnetosphere of Mercury, and its associated phenomena and their correlations with solar wind variations.

We used the AMITIS code, the first three-dimensional electromagnetic hybrid model of plasma that runs on a single CPU–GPU pair (Fatemi et al. 2017), to infer the solar wind plasma parameters upstream of Mercury from magnetic field observationsinside the magnetosphere. We selected three orbits of MESSENGER, shown in Fig. 1, that have been analysed and compared with hybrid simulations before. With the help of our high-performance model and by comparing our simulation results with magnetic field observations along the trajectory of MESSENGER inside the magnetosphere of Mercury, we estimated the upstream solar wind dynamic pressure, Alfvén Mach number, and plasma β during each of the selected orbits (Figs. 2 and 6). We showed in Figs. 4, 8, and 11 that our model provides good agreement with MESSENGER observations. Furthermore, our simulations, shown in Figs. 3, 7, 9, and 10, provide global aspects of the solar wind plasma interaction with Mercury during the selected orbits of MESSENGER including the bow shock, magnetopause, magnetotail, polar cusps, quasi-trapped particles, and partial ring current-like structures.

We compared the solar wind dynamic pressure inferred from our simulations with those predicted by the WSA-ENLIL model (Baker et al. 2013). For the D113 and the D182 orbits, the dynamic pressure inferred from our simulations is different to those predicted by the WSA-ENLIL model. However, as shown in Figs. 4 and 11, there is a good agreement between our simulations and the MESSENGER magnetic field observations during those orbits. This may confirm that the dynamic pressure inferred from our simulations is closer to the actual solar wind dynamic pressure, compared to that predicted by the WSA-ENLIL model.

In all the simulation results presented here, we have used regularly spaced Cartesian cubic grids of size Δ L = 200 km (~0.08 R_M), which is ~4δ_i, where δ_i = c/ω_pi is the ion inertial length. Generally, in a hybrid model of plasma ΔL ≫ δ_e, where δ_e = c/ω_pe is the electron inertial length. The typical δ_e near the orbit of Mercury is nearly one kilometre, which makes Mercury a suitable object for hybrid modelling. However, making Δ L much smaller than 200 km (e.g. 50 or even 100 km) when all simulation parameters are realistic, like ours, currently makes simulations computationally very expensive or not possible, even with the help of GPUs, when tens or hundreds of simulation runs are required, like in this study. This relatively large Δ L constrains our model–data comparison, and particularly affects the determination of the magnetospheric boundary crossings, even with the help of interpolation between neighbouring cells, which is the approach we have taken in this study. The finite cell size effect is evident in our simulation results presented in Fig. 12, where each cluster of our simulation points (red circles) are separated horizontally by the size of Δ L. Thus, we expect that the location of the magnetospheric boundaries presented here has an error comparable in size to Δ L.

Moreover, since the orbital speed of MESSENGER changes along its elliptical orbit around Mercury, we could not present our simulation results shown in Figs. 2 and 6 as a function of distance for the magnetospheric boundaries. Thus, we have presented them as a function of time, which allowed us to directly compare our simulations with observations.

In addition to the finite simulation cell size effect, we assumed a constant Keplerian speed for Mercury, compensated in our simulations by a 50 km s⁻¹ dawnward component for the solar wind velocity in the MSO coordinate system. Realistically, this speed would vary between ~30 to ~70 km s⁻¹ due to the eccentric orbit of Mercury, which results in <4% variations in the solar wind dynamic pressure. Therefore, a <4% error is involved in the solar wind dynamic pressures inferred from our simulations. In addition, this results in ≤6° variations in the solar wind aberration angle, which affects the location of the magnetospheric boundaries. We have examined the effect of the solar wind aberration angle in our simulation results (not shown here), and we found that it mainly affects the location of the bow shock downstream. Our simulations showed negligible variations in the upstream stand-off distance of the magnetopause and bow shock when the nearly 6° variations inthe solar wind aberration is taken into account.

Moreover, for simplicity and to reduce one level of freedom, we assumed a constant solar wind temperature in our simulations. We have also examined the effects of different solar wind plasma temperature in our simulation results (not shown here) and found that the temperature affects the global structure of the solar wind interaction with Mercury, mainly because the plasma β changes. Wealso found that the solar wind plasma temperature affects the location of the magnetospheric boundaries downstream and its effects become more noticeable at further distances, like the inbound magnetospheric crossing during the M1 flyby. However, due to the lack of the solar wind plasma information upstream, we could not provide any estimation for the upstream solar wind temperature, and thus assumed it to be constant in our simulations and chose it to be the typical plasma temperature near the orbit of Mercury.

In this study we inferred solar wind plasma dynamic pressure and Alfvén Mach number upstream of Mercury from magnetic field observations, mainly based on the magnetospheric boundary crossings. Although there is a good agreement between our simulations and observations shown in Figs. 4, 8, and 11, our estimations for the solar wind plasma parameters are highly dependent on the accuracy in the determination of the magnetospheric boundaries from observations. As an example, we show in Fig. 11 that our model suggests a different location for the outbound bow shock crossing compared to that estimated from observations.

In all of our simulations presented here, we assumed that the intrinsic magnetic field of Mercury is a single dipole field with a northward offset. However, MESSENGER observations have indicated signatures of remanent magnetization in Mercury’s crust, known as magnetic anomalies (Johnson et al. 2015; Hood 2016). The observed crustal magnetic anomalies are mainly near the Caloris impact basin near the north pole with a strength of a few nano-Tesla at 40 km altitude above the surface, and their radial component opposes the intrinsic magnetic dipole field of Mercury (Hood 2016). The crustal magnetic fields have not been considered in our simulations and they may affect the solar wind interaction with Mercury at close distances to the surface. Understanding the effects of the crustal fields on the solar wind interaction with Mercury is not within the scope of this paper.

The recent observations of Mercury’s gravity field and Earth-based radar observations of Mercury’s librations have suggested that Mercury has a large electrically conductive iron-rich core with a radius of ~2020 km (~0.8 R_M) encompassed by a ~420 km solid silicate mantle and crust (Smith et al. 2012; Hiremath 2012; Hauck et al. 2013). Large increases in the solar wind dynamic pressure considerably move the magnetopause and magnetotail current sheets closer to Mercury and intensify their current densities (Slavin et al. 2014; Jia et al. 2015). These time-varying changes in electric currents and their associated magnetic fields generate induced currents at the boundaries between the conductive core and the relatively resistive mantle of Mercury to cancel out magnetic field perturbations inside the core. Outside the conductive core, however, the induced currents temporarily increase magnetic fields in Mercury’s magnetosphere, from a few nano-Tesla to hundreds of nano-Tesla, depending on the solar wind dynamic pressure variations (Jia et al. 2015). It has been theoretically suggested that the induced magnetic fields keep the magnetopause at ~1.1–1.2 R_M (Hood & Schubert 1979; Suess & Goldstein 1979; Jia et al. 2015), and thus the solar wind would not easily access the surface of Mercury even during large dynamic pressure changes.

Jia et al. (2015), using a global MHD simulation, have studied the electromagnetic response of Mercury’s interior for large solar wind dynamic pressure changes. According to their simulations, a solar wind dynamic pressure change from ~11 nPa to ~66 nPa can induce ~200 nT fields at the surface of Mercury, which is comparable to or even larger than the strength of the surface field generated by the intrinsic magnetic dipole moment. As numerically shown and discussed by Jia et al. (2015), the induced field influences the structure of the magnetosphere by moving the magnetopause ~0.1 R_M further upstream compared to the conditions when the interior magnetic induction is not taken into account (e.g. Siscoe & Christopher 1975; Slavin & Holzer 1979; Slavin et al. 2009). Thus, the surface of Mercury is more shielded from direct impact of the solar wind during high dynamic pressure changes when the interior conductivity is considered. However, as discussed in detail by Slavin et al. (2014), the stand-off distance of the magnetopause during extreme solar events observed by MESSENGER (shown by the green line and asterisks in Fig. 12) is closer to the surface than the values predicted before (e.g. Hood & Schubert 1979). Slavin et al. (2014) suggest that the intense magnetic reconnection at Mercury may oppose the shielding effects of the induced currents, which moves the subsolar magnetopause much closer to Mercury.

However, in our simulations, we did not account for the conductive interior of Mercury. Instead, we assumed that Mercury is a uniform resistive obstacle to the solar wind. This is a valid assumption for this study, and there are two particular reasons why the conductive interior of Mercury does not affect our simulation results presented here: (1) we assumed that the solar wind plasma and IMF remain constant in time upstream of Mercury during our simulations; therefore, there is no time-varying magnetic field and/or plasma dynamic pressure to generate interior induced magnetic fields, and (2) MESSENGER magnetic field observations during the three selected orbits analysed here (i.e. the M1 flyby, D113, and D182) do not show any large changes in the IMF before (after) MESSENGER moved into (out of) the magnetosphere. Therefore, we expect that there is either no contribution or very minimal contribution from the conductive core of Mercury in the MESSENGER orbits analysed here. Perhaps there are still some remnants from a previously induced magnetic field from solar wind and IMF changes which are diffusing into the interior over time. We do not completely reject the possibility that the solar wind pressure or IMF changes during the long period (a few hours) that MESSENGER was passing through the magnetosphere for the orbits analysed here, but there is no convincing evidence of such changes in the observed magnetic fields presented in Figs. 4, 8, and 11.

We conclude that during quiet times for the solar wind, i.e. no large changes in the solar wind plasma or IMF, our model-estimated fitted function (Eq. (2)) shown by the red line in Fig. 12 provides a valid determination of Mercury’s magnetopause stand-off distance as a function of the solar wind dynamic pressure. Intriguingly, this function agrees fairly well not only with the lower limit of the empirical function obtained from MESSENGER for low solar wind dynamic pressures (light blue dashed line in Fig. 12), but also with MESSENGER observations during extreme solar events (green asterisks in Fig. 12 from Slavin et al. 2014). Since our model is a kinetic model (solar wind ions are particles), magnetic reconnection is self-consistently present in our simulations, which may result in erosion of the magnetosphere and may change the magnetopause stand-off distance, as suggested previously by Slavin et al. (2014). However, studying the magnetic reconnection and its effects on the magnetosphere of Mercury is very important, and is left for future investigations.

In addition to the magnetopause response to the solar wind plasma variations, Winslow et al. (2013) have analysed the response of the bow shock to the solar wind. Since the IMF orientation for the majority of the simulations we have presented here is nearly quasi-parallel to the solar wind flow, the foreshock region formed near the subsolar point (e.g. evident in Figs. 7b,9b) does not allow us to study the bow shock response to solar wind variations from our simulations.

Mercury has a very thin, surface-bound, and collision-less neutral exosphere observed by Mariner 10, MESSENGER, and ground-based observations from the Earth (e.g. Broadfoot et al. 1974; Potter & Morgan 1985; Zurbuchen et al. 2008; Doressoundiram et al. 2009). Ionization of exospheric neutrals through various processes including charge exchange with solar wind ions, photoionization, and electron impact ionization produces planetary ions (e.g. Killen et al. 2008). These ions, mainly Na⁺ and O⁺ groups, have been observed to be substantially enhanced at high latitudes over the north pole, with higher densities on the dayside than on the nightside (~0.7 cm⁻³ for Na⁺-group and ~0.05 cm⁻³ for O⁺ -group on the dayside, compared to20–30 cm⁻³ solar wind density) (e.g. Killen et al. 2001; Zurbuchen et al. 2011). In our model we only included a single solar wind plasma species and did not account for the tenuous planetary ions. Even so, the simulation results presented here are in good agreement with the MESSENGER observations, mainly because the planetary ions have negligible density compared to the solar wind plasma density at high altitudes, i.e. along the orbit of MESSENGER. Thus, they cannot significantly perturb the magnetic fields of Mercury. However, they make critical contributions to processes at close distances to Mercury and may affect the local electric and magnetic fields, which remain open for future investigation.

5 Summary and conclusions

With the aid of our high-performance hybrid model of plasma (the AMITIS code), we inferred the upstream solar wind plasma parameters (i.e. dynamic pressure and Alfvén Mach number) and their variations over time from magnetic field observations inside the magnetosphere of Mercury. We showed that for the inferred upstream parameters, there is a good agreement between our simulations and magnetic field observations along the trajectory of MESSENGER. We have built a large library of different solar wind plasma parameters and IMF configurations to study the solar wind interaction with the magnetosphere of Mercury. This library may also help us to estimate solar wind plasma parameters upstream of Mercury for other orbits of MESSENGER that have not been analysed here. Knowledge of the solar wind plasma and its variations is essential in order to understand a highly dynamic magnetosphere like Mercury’s. It is also crucial in order to interpret observationsfrom the missions to Mercury that cannot directly observe the solar wind plasma, like MESSENGER and BepiColombo.

Acknowledgements

All MESSENGER data used here are publicly available at NASA’s Planetary Data System (PDS). From there we have used the MESSENGER MAG reduced data archive (MESS-E/V/H/SW-MAG-4-SUMM-CALIBRATED-V1.0) in the MSO coordinate system, file names MAGMSOSCIAVGXXXXX_10_V08. Data from the modelling results are available upon request to the corresponding author, Shahab Fatemi (shahab@irf.se or shahabfatemi@gmail.com). This research was partly conducted using resources provided by the Swedish National Infrastructure for Computing (SNIC), project SNIC 2017/11-15 at the High Performance Computing Center North (HPC2N), Umeå University, Sweden.

References

All Tables

All Figures

Fig. 1 Portion of MESSENGER’s trajectory around Mercury during the M1 flyby (red line) and two nominal orbits of MESSENGER on 23 April 2011, D113 (green line), and on 01 July 2011, D182 (blue line) in a cylindrical MSO coordinate system. The arrows show the direction of MESSENGER’s motion along each orbit. Superimposed are the approximate locations of the magnetopause (thick dashed line) and bow shock (thin dashed line) obtained from Winslow et al. (2013).

In the text

Fig. 2 Time difference between our hybrid simulations and the mid-point location of the magnetospheric boundary crossings obtained from MESSENGER magnetic field observations by Winslow et al. (2013). The dashed line is a linear fit of our results using an ordinary least-squares fitting function. Panels a and b: inbound bow shock and magnetopause crossings. Panels c and d: Outbound magnetopause and bow shock crossings for the second magnetospheric transient on 23 April 2011 (D113). Pdyn = is the solar wind dynamic pressure and MA = is the Alfvén Mach number.

In the text

Fig. 3 Hybrid simulation results for the magnetospheric transit of MESSENGER on 23 April 2011 between 16:00 and 21:00 (D113) presented in the MSO coordinate system. Panels a and b: magnitude of the magnetic field in logarithmic scale. Panels c and d: electric current density calculated from the general Ampère’s law, flow normal to the presented planes. Panels e and f: plasma density in logarithmic scale and normalized to the upstream solar wind density, nsw = 22 cm−3. Panels g and h: magnitude of the solar wind velocity normalized to the upstream solar wind velocity, |usw | = 314 km s−1. The top panels are cuts in the equatorial plane (xy-plane at z = 0), viewed from Mercury’s north pole, and the bottom panels are cuts in the midnight meridian plane (xz-plane at y = 0), viewed from the orbital motion of Mercury (i.e. the -y axis). Mercury is shown by a circle, centred at the origin of the coordinate system. The direction of the solar wind and the IMF are shown by yellow and white arrows, respectively, in panels a and b. Streamlines in panels c and d show magnetic field line tracing, and in panels g and h show the plasma flow direction. A portion of MESSENGER’s orbit on 23 April 2011 (D113), also shown in Fig. 1, is shown in all panels and the UTC time during the passage is marked by arrows in panels a and b.

In the text

Fig. 4 Panels a–d: magnetic field comparison between our hybrid model simulations (red lines), MESSENGER magnetometer observations (black lines), and undisturbed intrinsic magnetic dipole of Mercury (blue dashed lines) along the trajectory of MESSENGER on 23 April 2011 (D113) between 16:00 and 21:00 UTC. Panel e: magnitude of the electric current density calculated from our simulations using the general Ampère’slaw. Panel f: solar wind bulk flow speed normalized to the upstream solar wind speed |usw | = 314 km s−1 (purple line) and solar wind plasma density normalized to the upstream plasma density nsw = 22 cm−3 obtained from our hybrid model simulations along the trajectory of MESSENGER. The mid-point location of the bow shock (BS) and magnetopause (MP) boundaries estimated by Winslow et al. (2013) as well as the closest approach (CA) to the planet are shown by the vertical lines.

In the text

	Fig. 5 Hybrid simulation comparison with MESSENGER observations for the upstream solar wind dynamic pressure ~2.7 nPa, similar to our model estimated dynamic pressure during the inbound magnetospheric crossing shown in Fig. 2. The figure format is the same as that shown in Fig. 4.
In the text

Fig. 6 Time difference between our hybrid simulations and the location of the magnetospheric boundary crossings during the M1 flyby obtained from MESSENGER magnetic field observations by Slavin et al. (2008). The dashed lineis a linear fit of our results using an ordinary least-squares fitting function. Panels a and b: inbound bow shockand magnetopause crossings. Panels c and d: outbound magnetopause and bow shock crossings. Pdyn = is the solar wind dynamic pressure and MA = is the Alfvén Mach number.

In the text

	Fig. 7 Global structure of the Hermean magnetosphere during the M1 flyby obtained from our hybrid simulations in the same formatas Fig. 3.
In the text

Fig. 8 Comparison between our hybrid simulations and MESSENGER observations along the trajectory of MESSENGER on 14 January 2008 during the M1 flyby. The format of the presented results is the same as that shown in Fig. 4. Panel f: upstream solar wind speed |usw | =365 km s−1 and solar wind plasma density nsw = 32 cm−3. Panel g: plasma counts observed during the M1 flyby by the FIPS instrument on MESSENGER, digitized from Fig. 3 in Raines et al. (2011).

In the text

Fig. 9 Hybrid simulation results for the magnetospheric transit of MESSENGER on 01 July 2011 (D182) for the IMF conditions (top panels) before the inbound and (bottom panels) after the outbound bow shock crossings. Panels a and b: magnitude of the magnetic field in logarithmic scale. Panels c and d: electric current densities calculated from the general Ampère’s law moving along the y-axis, Jy, perpendicular to the plane shown here. Panels e and f: plasma density in logarithmic scale and normalized to the upstream solar wind density, nsw = 30 cm−3. Panels g and h: magnitude of the solar wind velocity normalized to the upstream solar wind velocity, |usw | ≈ 385 km s−1. All the panels are cuts in the midnight meridian plane (xz-plane at y = 0), viewed from the orbital motion of Mercury (i.e. the -y axis). The rest of the figure format is the same as that shown in Fig. 3.

In the text

Fig. 10 Hybrid simulation results for the D182 orbit in the planes perpendicular to those shown in Fig. 9. Panels a–d: for the inbound and panels e–h: for the outbound IMF conditions. The cuts are in the yz-plane panels a,b,e and f: at x = +1.0RM, and panels c,d,g and h: at x = –1.5 RM. Panels a,c,d and e: magnitude of the magnetic field in logarithmic scale. Panels b,d,f and h: plasma flux normalized to the upstream solar wind flux Γsw = 1.15 × 1013 m−2⋅s. The solar wind flows into the planes shown here with 50 km s−1 along the +y axis and the direction of the IMF in the yz-plane are shown by black arrows in panels a and e.

In the text

Fig. 11 Comparison between our hybrid simulations and MESSENGER observations along the trajectory of MESSENGER on 01 July 2011 between 06:00 and 10:00 UTC (D182). Panels a-d: magnetic field comparisonbetween MESSENGER observations (black lines) and our hybrid simulations for the inbound (red lines) and outbound (green lines) IMF conditions listed in Table 1 along the trajectory of MESSENGER. To provide better visualization, we have subtracted the intrinsic magnetic fields of Mercury, BD, from observations and from simulations. Panel e: electric current density. Panel f: plasma density normalized to the upstream solar wind plasma density nsw = 30 cm−3. Panel g: plasma velocity normalized to the upstream solar wind plasma speed |usw | ≈385 km s−1, obtained from our simulations along the trajectory of MESSENGER. The mid-point location of the bow shock (BS) and magnetopause (MP) boundaries estimated by Winslow et al. (2013) as well as the closest approach (CA) to the planet are shown by the vertical lines.

In the text

Fig. 12 Stand-off distance of the magnetopause, Rss, from the centre of Mercury at the subsolar point in the MSM coordinate system as a function of the solar wind dynamic pressure, Pdyn. The solid blue line is from a MESSENGER empirical model by Winslow et al. (2013) shown in Eq. (1), where the lower limit and upper limit of the model are marked by the light blue and dark blue dashed line, respectively. The acceptance area between the lower and upper limits is in yellow. Taken from Fig. 19 by Slavin et al. (2014), the green asterisks are MESSENGER observations obtained during three extreme solar wind dynamic pressure events (i.e. CMEs and HSS) and the green dashed line is the function fitted to the observations by Slavin et al. (2014). The red circles are from our hybrid simulations and the red solid line is the curve fitted to our simulations, shown in Eq. (2), using a least-squares curve fitting function.

In the text