Influence of Mercury's Exosphere on the Structure of the Magnetosphere

Mercury is embedded in a tenuous and highly anisotropic sodium exosphere, generated mainly by plasma‐surface interactions. The absolute values of the sodium ion density are still under debate. Observations by MESSENGER's Fast Imaging Plasma Spectrometer (FIPS) instrument suggest the density of exospheric ions to be several orders of magnitude lower than the upstream solar wind density, indicating that the sodium exosphere has no substantial influence on the magnetospheric current systems. However, MESSENGER magnetic field observations of field line resonances revealed sodium ion densities comparable to the upstream solar wind density. To investigate how a dense exosphere would affect the current systems within Mercury's magnetosphere, we apply an established hybrid (kinetic ions, fluid electrons) model and conduct multiple model runs with gradually increasing exospheric density, ranging from no sodium ions at all to comet‐like configurations. We demonstrate how a sufficiently dense exosphere leads to self‐shielding of the sodium ion population from the ambient electric field and a significant inflation and symmetrization of Mercury's magnetosphere, which is decreasingly affected by the dipole offset. Once the sodium ion density is sufficiently high, Region 2 field‐aligned currents emerge close to the planet. The modeled Region 2 currents are located below the orbit of MESSENGER, thereby providing a possible explanation for the absence of these currents in observations. The sodium exosphere also closes a significant fraction of the Region 1 currents through Pedersen and Hall currents before the “guiding” magnetic field lines even reach the planetary surface. The modeled sodium ion and solar wind densities agree well with observations.


Introduction
Mercury, the first planet in our solar system, is the smallest in size with a radius of R M = 2,440 km (Smith et al., 2012) and the densest with its iron core taking up about 80% of its radius (Johnson et al., 2016). Magnetic field data obtained from the Mariner 10 and MESSENGER missions showed that the molten core generates a weak intrinsic magnetic field Glassmeier, 1997;Slavin et al., 2008). Mercury's magnetic field at the equator is about 190 nT. It can either be described by a dipole that is antiparallel to its rotation axis, shifted northward by 0.2 R M  or equivalently by a centered dipole superimposed with a significant quadrupole (Wicht & Heyner, 2017). The weak dipole strength can possibly be explained by a feedback coupling between Mercury's conducting interior and its magnetospheric environment (Glassmeier et al., 2007;Heyner et al., 2011).
The interaction of the weak Hermean magnetic field with the impinging solar wind results in a small magnetosphere that qualitatively resembles the Earth's magnetosphere. Compared to Earth's Dungey timescale of about an hour, Mercury's Dungey timescale is a factor of 30 shorter (Slavin et al., 2009. The high reconnection rate at Mercury is independent of the upstream interplanetary magnetic field (IMF) direction, and thus, the Hermean magnetosphere is highly dynamic and susceptible to sudden variations of the upstream solar wind conditions (DiBraccio et al., 2013). Instead, reconnection at the Hermean dayside magnetopause rather seems to be controlled by the difference in plasma β across the magnetopause (Poh et al., 2016).
Earth's magnetotail consists of the northern and southern tail lobes, which are nearly devoid of plasma. The magnetic field of the lobes originates from the polar regions of the planet. At low magnetic latitudes, that is, near the equatorial plane, these lobes are divided through the central plasma sheet, while the northern and southern plasma mantle envelopes the respective lobes at the high latitudes (Frank, 1985). The plasma mantle is populated by planetary ions originating from the high latitudes of the ionosphere and magnetosheath plasma that has entered the magnetosphere through the polar cusp regions. How much solar wind plasma is able to enter through the cusps is heavily dependent on the B z -component of the IMF (Rosenbauer et al., 1975). The ions within the plasma mantle flow downstream, parallel to the enveloping magnetopause. However, as the ions experience an E×B drift toward the plasma sheet, they become stratified depending on their velocity parallel to the magnetic field. Thus, slow-moving particles enter the plasma sheet close to the planet, while fast-moving particles enter the plasma sheet at farther distances or are lost down the tail (Frank, 1985;Pilipp & Morfill, 1978;Rosenbauer et al., 1975). At Mercury, first observations of the southern plasma mantle have been presented by DiBraccio et al. (2015). These authors showed how particles from the southern plasma mantle are a significant contributor of solar wind plasma for the plasma sheet. In a subsequent study of 1,051 MESSENGER transits through the southern magnetosphere, Jasinski et al. (2017) showed that the southern plasma mantle was encountered in only 94 cases. Their analysis showed how the presence of the southern plasma mantle is facilitated by an antisunward B x -component of the IMF and when the magnitude of the IMF is below 25 nT , which is about the average IMF magnitude at Mercury (Winslow et al., 2013). Furthermore, the sign of the IMF's B y -component is linked to the plasma mantle being predominantly present in the duskward or dawnward side of the magnetosphere (Jasinski et al., 2017). This is due to the magnetospheric tail exhibiting a tilting behavior depending on the B y -component of the IMF (Sibeck et al., 1985). As observations of the northern plasma mantle are not feasible with MESSENGER's orbit configuration, future missions such as BepiColombo will enable the analysis of these regions (Benkhoff et al., 2010).
Mercury's surface is embedded in a dilute gas envelope, consisting mainly of sodium (Cheng et al., 1987;Milillo et al., 2005). This sodium exosphere is generated by four major processes: thermal desorption (TD), photon-stimulated desorption (PSD), micrometeoroid impact vaporization (MIV), and surface sputtering (SP). Through Earth-bound telescope observations, Mangano et al. (2015) found that sodium doublet emission peaks appear at high latitudes in both the northern and southern hemispheres, that is, in regions where SP and PSD processes are most intense or close to the equator where TD dominates. Orsini et al. (2018) showed that the latitudes with the most sodium emissions are controlled by the space weather at Mercury, and the low-latitude emissions occur when the planet is hit by interplanetary coronal mass ejections (ICMEs). Recently, analysis of MESSENGER data showed how ICMEs with the highest dynamic pressure ever recorded and intense southward magnetic field components were able to push and erode the magnetopause below Mercury's surface . In such cases, the dayside surface is exposed to direct solar wind bombardment, which locally increases the sputtering rates (Winslow et al., 2020).
At Earth, the ionosphere is an important part of the global magnetospheric current system. The midlalitude to high-latitude current structure within Earth's magnetosphere connects the magnetopause currents to the planet through Region 1 field-aligned currents (R1-FACs or Birkeland currents) and Region 2 field-aligned currents (R2-FACs). The large height-integrated Pedersen conductivity, that is, Pedersen conductance, within Earth's dense ionosphere reaches values of about 1-10 S. This conductance prevents currents from reaching the terrestrial surface; that is, they are closed through perpendicular currents inside the ionosphere (Kivelson & Russell, 1995). For the case of Mercury, Cheng et al. (1987) used Mariner 10 observations to estimate the Pedersen conductance (resulting from pickup processes of sodium ions) to be about 0.1-0.3 S, that is, multiple orders of magnitude lower than the values found at Earth. Leblanc et al. (2013) showed that the velocity of sodium neutrals is about three orders of magnitude smaller than the bulk velocity of the solar wind. The tenuous sodium exosphere is partially photoionized, and sodium ions with an average density in the range of 5.1·10 −3 -2 cm −3 have been observed by the Fast Imaging Plasma Spectrometer (FIPS Andrews et al., 2007) throughout Mercury's magnetosphere (Dewey et al., 2018;DiBraccio et al., 2015;Gershman et al., 2014;Korth et al., 2014;Raines et al., 2014Raines et al., , 2015. These in situ observations suggest the density of the exospheric ions to be multiple orders of magnitude lower than the average upstream solar wind proton density of about 60 cm −3 (Winslow et al., 2013), thereby suggesting that the sodium exo-ionosphere has no dynamic influence on the current systems within the magnetosphere. However, FIPS observations rather provide a lower limit on ion densities, due to nonideal spacecraft attitudes and ions outside of the observable energy range of 100 eV/e to 13.3 keV/e , which would be especially important for recently generated, cold exospheric ions yet to be picked up by the ambient electric field.
The notion of a much higher sodium ion density inside Mercury's magnetosphere was put forward by a recent study by James et al. (2019). These authors used MESSENGER magnetic field data and analyzed 566 instances where the frequencies of the local magnetic field lines resonate with the local plasma frequency. By using a power-law model with the exponent as a free parameter for the plasma mass density profile along the field lines and the KT17 model of Mercury's average magnetospheric field (Korth et al., 2015, James et al. (2019) obtained an estimate of the plasma mass density close to the surface of about 500 amu/cm 3 . Furthermore, James et al. (2019) argued that the plasma close to the surface consists mostly of sodium ions. Therefore, their derived plasma mass density would correspond to a sodium ion number density of about 22 cm −3 , which is much closer to the average upstream solar wind density than the sodium ion densities observed by FIPS. However, James et al. (2019) had to restrict their analysis to the closed planetary field lines of the KT17 model in order to be able to establish boundary conditions that lead to standing waves. Thus, James et al. (2019) were not able to obtain any information on the plasma densities for the cusp and remote nightside regions. Taking into account both FIPS observations and the study of James et al. (2019), the values of the sodium ion density in Mercury's magnetosphere are a subject of ongoing discussion and may differ by an entire order of magnitude between different approaches. Due to the large gyroradii of the sodium ions, a more dense sodium exo-ionosphere may make nonnegligible contributions to magnetospheric current systems. In particular, increased ion density will severely alter the near-surface conductivity and thus have a critical influence on the closure of FAC systems. Anderson et al. (2014) subtracted Mercury's dipole field  and average external field disturbances (from magnetopause and tail currents) from MESSENGER data and calculated the FACs from the magnetic field residuals in the northern hemisphere. Their analysis reveals two regions on the high-latitude dawn and dusk side with equal amplitude currents flowing radially toward and away from the planet, respectively. These currents were identified as the Hermean R1-FAC system. The associated closure currents have to flow from dawn to dusk either through the exosphere or inside the planet. Using a shell model for Mercury's mantle with a radially symmetric conductivity profile and assuming the conductance of the exosphere to be negligible, Anderson et al. (2014) proposed that these internal currents have a radial path through the upper crust and flow in lateral direction at the core-mantle boundary. The mantle's height-integrated conductivity, derived from the magnitude of the R1 currents, has been estimated to 1 S, which is about a factor of 3-10 larger than the exospheric Pedersen conductance derived by Cheng et al. (1987). However, with regard to the uncertainties in the sodium ion densities close to the surface, it is reasonable to assume that the exospheric ions increase the low-altitude conductivity outside of the planet and partially close the FACs before they ever penetrate Mercury's surface, thereby drastically reducing the amount of current connecting to the planetary interior. Anderson et al. (2014) also argued that the apparent absence of R2-FACs at the altitudes covered by MESSENGER (between 400 and 1,000 km) implies that R2-FACs exist either at lower altitudes or within latitudes of ±30°from the magnetic equator.
The goal of our study is to assess to what extent an increased sodium ion density affects the overall structure of Mercury's magnetosphere. We apply an established hybrid model (ions kinetic, electron fluid), including a sophisticated latitude-dependent exosphere model that takes the four processes TD, MIV, PSD, and SP into account. With a hybrid model, we can directly assess the effects of the sodium ion pickup and the associated modification of the electric and magnetic fields. The total sodium density of the exosphere will be used as a free parameter, considering surface densities derived from an exosphere model by Gamborino et al. (2019), the results of James et al. (2019) and even larger values that would correspond to Venus/Mars-like interaction scenarios. By comparing the runs with different exospheric densities, we will investigate what densities are needed to actually affect the magnetospheric current systems. Furthermore, we will assess whether the effects of a dense exosphere could be observed in the altitudes covered by the Mercury Magnetospheric Orbiter (MPO), which is part of the upcoming BepiColombo mission (Benkhoff et al., 2010;Glassmeier et al., 2010).
In section 2, we introduce the model for the interaction of Mercury's magnetosphere with the upstream solar wind. The latitude-dependent exosphere model is presented in section 3. We investigate how a variable sodium exosphere alters the global structure of the solar wind ions, sodium ions, and the current system within the Hermean magnetosphere in section 4. In section 5, we summarize the results of this study.

Hybrid Model AIKEF
To model the feedback of sodium ions on Mercury's magnetosphere, we use the three-dimensional hybrid model AIKEF (Müller et al., 2011) which treats electrons as a massless, charge-neutralizing fluid, and the ions kinetically. The AIKEF model has already been applied to Mercury by Müller et al. (2012) and Exner et al. (2018). The former study focused on the "double magnetopause" signature that was identified in magnetic field data of the first two MESSENGER flybys (Slavin et al., 2008). Müller et al. (2012) showed that this current sheet feature developed as solar wind protons were diverted into the dayside magnetosphere by nightside reconnection events. Exner et al. (2018) investigated the dynamic interaction of Mercury's magnetosphere with an ICME and found that at least two different sets of solar wind upstream conditions are required to explain magnetic field observations from the dayside and nightside magnetosphere, indicating that the upstream conditions changed severely on timescales shorter than a MESSENGER orbit.
Under nominal solar wind conditions, the gyroradii of solar wind protons at Mercury are on the order of 0.08 R M (Winslow et al., 2013), and the motion of the protons can be described in the guiding center picture. The use of a hybrid model to capture proton dynamics is therefore not imperative. However, when considering planetary heavy ions with their gyroradii comparable to Mercury's radius, individual ion motion needs to be taken into account. Both hybrid (Fatemi et al., 2018;Herčík et al., 2013;Müller et al., 2012;Parunakian et al., 2017;Richer et al., 2012;Trávníček et al., 2010) and magnetohydrodynamics (MHD) models Pantellini et al., 2015;Varela et al., 2016;Yagi et al., 2010) have been applied to describe solar wind dynamics near Mercury. As we analyze the feedback of exospheric sodium ions onto Mercury's magnetosphere, which have a pickup gyroradius of about 2.4 R M , the application of a hybrid model is indeed mandatory.
In the AIKEF model, the static, Cartesian three-dimensional grid is defined with respect to the Mercury Anti Solar Orbital (MASO) coordinate system, where the origin is at the center of the planet, and the x-axis is aligned with the upstream solar wind velocity. The z-axis is antiparallel to Mercury's dipole moment, while the y-axis completes the right-handed system, thereby pointing (roughly) in the direction of Mercury's orbital motion. All model runs presented in this study use the same box dimensions of 20 R M in x and z and 12 R M in y direction with Mercury located a distance of 7 R M from the inflow plane. Each cubic cell has an edge length of 0.056 R M which is more than sufficient to resolve the exospheric sodium ion as well as solar wind proton gyroradii. At the beginning of the simulations, each cell is filled with 20 macroparticles per species, which have the same charge-to-mass ratio as the ions they represent. Therefore, our model is applicable to assess how the unusual high sodium ion densities found by James et al. (2019) would affect the magnetospheric dynamics.
As a representative setup for the upstream conditions at Mercury, we use average values for the IMF magnitude (Winslow et al., 2013) and electron and ion temperatures (Pierrard et al., 2016). We use the same IMF direction as employed in Exner et al. (2018). This IMF is mostly parallel to the y-axis, which is a typical IMF direction at Mercury (Winslow et al., 2013). The corresponding electric field therefore mostly points northward. Therefore, pickup ions will be mostly confined to the xz-plane, thereby facilitating a more straightforward analysis. The results of Exner et al. (2018), which did not include an exosphere, may also be considered as an additional baseline to assess the influence of the sodium exosphere in the simulations presented here. We employ a low dynamic pressure of 2.1 nPa of the impinging solar wind, which is about a factor of 7 smaller than the observed average value (Winslow et al., 2013). Consequently, Mercury's magnetosphere will be inflated compared to the nominal case and the influence of sodium ions on the magnetosphere, and current systems will be more accessible. The corresponding solar wind density and velocity are n 0 = 10 cm −3 and 500 km/s, respectively. These conditions correspond to an Alfvénic Mach number of M A = 3.25 and a plasma β of 0.18. As discussed in Müller et al. (2011), AIKEF includes a self-consistent Note. Physical parameters of Mercury are obtained from Anderson et al. (2012) and Winslow et al. (2013).
coupling between Mercury's space environment and the diffusion of the magnetic field in its interior. The inner resistivity profile is modeled in agreement with geophysical models (Anderson et al., 2018;Hauck II et al., 2013). This approach has been successfully applied in Jia et al. (2015Jia et al. ( , 2019 and Exner et al. (2018). All parameters have been summarized in Table 1.

Exosphere Model
The sodium exosphere of Mercury is generated by four major processes: TD, MIV, PSD, and SP; see Killen et al. (2007); Milillo et al. (2005); Raines et al. (2015); Gamborino and Wurz (2018); Gamborino et al. (2019). Since they are caused by the solar radiation, TD and PSD are confined to the dayside, while bombardment with micro-meteoroids yields, on average, a radially symmetric release of sodium particles (Wurz et al., 2010). Proton SP mainly happens at the footpoints of the cusp regions. The cusp locations are controlled by the upstream solar wind conditions (He et al., 2017;Winslow et al., 2012). Under nominal solar wind conditions, four distinct surface regions with high solar wind influx have been identified: two regions in the dayside hemisphere at high northern and southern latitudes (below the polar cusps) and two regions associated with the open-closed field line boundary in the nightside magnetosphere Rong et al., 2018;Wang & Ip, 2011). In times of increased solar wind ram pressure, for example, during the passage of ICMEs, the cusps are being pushed to lower latitudes and experience an enhanced particle influx (Winslow et al., 2013(Winslow et al., , 2017, or the magnetopause is even pushed below the surface Winslow et al., 2020). In such instances, the dayside surface is under direct bombardment of the solar wind, increasing the local SP rates.
Still, with observations showing that sodium ion density is smaller than the solar wind density by 2-4 orders of magnitude (Milillo et al., 2005;Raines et al., 2015), sodium ions were treated as test particles in previous simulation studies. Using an MHD approach for Mercury's magnetosphere, Seki et al. (2013) found that with high surface conductivity, near-surface plasma convection is inhibited; ions are not able to leave into the magnetosphere and are confined near the surface.
Employing an MHD model with low solar wind ram pressure, Yagi et al. (2017) showed that sodium ions were able to gyrate around Mercury to form a "sodium ring," while, at higher upstream ram pressure, only a partial sodium ring was formed at the nightside. Paral et al. (2010) combined a hybrid model of the electromagnetic fields with a test particle model consisting of PSD and SP. These authors showed that depending on the B z -component of the upstream magnetic field and resulting reconnection patterns, planetary ions are confined to low surface altitudes or are able to travel to higher altitudes and subsequently downstream for positive and negative B z -components, respectively. Gamborino et al. (2019) conducted a Monte-Carlo simulation for the exosphere which included gravitational escape, radiation pressure, and surface adsorption as loss processes. They were able to explain the sodium column density profile derived from MASCS/UVVS measurements on 23 April 2012 . Gamborino et al. (2019) obtained values for the scale height, surface density, and column density for each process (see Table 2). In agreement with Wang and Ip (2011), PSD is associated with the highest column density of 8.0·10 15 m −2 .
For our exosphere model, we adopt the parameters from Gamborino et al. (2019). As Mercury's gas envelope is collision less , a barometric law is not suitable. Therefore, our model representation of Mercury's exosphere is based on an empirical analytical model that was initially developed by Saur et al. (2008) and Roth et al. (2014) to describe local exospheric inhomogeneities at Enceladus and Europa (see also Arnold et al., 2019;Kriegel et al., 2011). The neutral density profile is described by a superposition of one radially symmetric profile and seven Gaussian bells, with their peaks located at different latitudes. The individual contributions represent MIV, TD, four SP, and two PSD source regions. The combined neutral gas profile is given by (1) with r ! being the radial distance from the origin of the MASO system, n MIV,0 and n i,0 , H r,MIV , and H r,i representing the surface density and the scale height of each source process, respectively. Δθ i is the angle between the local zenith direction at the maximum source density (the peak of the respective bell). H Θ,i is representing the angular width of each bell, and therefore, each profile dominates in a particular surface area. The latitudes of the maxima and angular widths were empirically estimated from Fig. 6 in Gamborino et al. (2019). The profile is assumed to be hemispherically symmetric with respect to the xzplane, that is, the planetary equator. All parameters are are given in Table 2.
We choose our surface densities n i,0 such that our modeled column densities are consistent with the column densities given by Gamborino et al. (2019). Furthermore, considering our hybrid model's peak resolution of 136 km near the surface, we need to increase the scale height term associated with TD process from 57 (Gamborino & Wurz, 2018) to 100 km to adequately resolve that contribution.
The resulting neutral profile is shown in Figure 1 and shows highly anisotropic features. Close to the dayside surface, PSD and TD make the strongest contribution to the local neutral density. However, due to their larger scale heights, MIV and SP dominate the neutral density at altitudes above 2 R M above the surface at the dayside. At high latitudes on the nightside, MIV is the only process for generating planetary ions. Huebner and Mukherjee (2015) showed that the photoionization frequency f Na,Earth of sodium atoms is in the range of 7.26-7.91·10 −6 s −1 at Earth's average orbital position for low and high solar activity, respectively. The ionization frequency varies along Mercury's highly elliptic orbit. To obtain the value at Mercury's average orbital position, we rescale this frequency f Na,Mercury =f Na,Earth (1/0.387 AU) 2 = 4.85-5.28·10 −5 s −1 . Within Mercury's geometric shadow, neutrals are ionized by charge exchange and electron impacts. We assume the ionization frequency within the shadow to be one fifth of the dayside value.
Even though Mercury's exosphere is collision less, the ionization of neutral particles and the subsequent pickup of the sodium ions can be represented by an effective collision frequency, associated with a Pedersen conductivity (Goertz, 1980;Neubauer, 1998). With sodium column densities obtained by Mariner 10, Cheng et al. (1987) found values of about 0.1-0.3 S for the Pedersen conductance of Mercury's sodium exosphere. To estimate the Pedersen conductance of each process from our exosphere model, we use the respective column densities, the updated values of the planetary dipole strength obtained by MESSENGER, and a photoionization frequency of 5.28 · 10 −5 s −1 , as the MESSENGER mission took place around solar maximum. In this way, we find that the PSD process in the high-latitude dayside regions corresponds to the maximum Pedersen conductance of 0.982 S, which is an order of magnitude larger than the value obtained by Cheng et al. (1987). The other processes of TD, MIV and SP contribute, in decreasing order, a Pedersen conductance of 0.103, 3.06·10 −4 , and 8.80·10 −5 S in the latitude regions where they are predominant. In other words, within the high-latitude regions, the Pedersen conductance, contributed by the PSD process, is of similar magnitude as the conductance of the mantle of 1 S derived by Anderson et al. (2014). Additionally, when the solar wind interaction with the magnetosphere is considered, the actual heavy ion densities may even be higher than assumed here (Vernisse et al., 2017; see also section 4); that is, the actual conductance may even exceed our estimates. It is therefore reasonable to assume that the inclusion of an exosphere might lead to a significant modification of current closure in these regions.
To study the effect of different planetary ion densities on the magnetic fields and currents in Mercury's magnetosphere, we carry out multiple model runs, treating the neutral density of the sodium exosphere as a free parameter. In the five model scenarios considered, the "baseline" neutral gas density from the Gamborino model is multiplied with factors of 0, 1, 5, 50, and 500, respectively.

Simulation Results
In order to illustrate the influence of sodium ions on the structure of Mercury's magnetosphere, we show hybrid model results in Figures 2 and 3 that are organized as follows: From left to right, the columns of the subfigures display the results from the model runs for sodium neutral profiles that are factors of 0, 1, 5, 50, and 500 of the profile defined by Equation 1, respectively. As the x-axis is aligned with the upstream solar wind velocity and the IMF is mostly aligned with the y-axis, the convective electric field points northward, and thus, the pick-up ions are initially accelerated mainly in the positive z direction. Therefore, in the following figures, we show the model output within the xz-plane. The model output for the xy-plane is similar to the results discussed in Exner et al. (2018).
From top to bottom row, Figure 2 shows the sodium ion density, the sodium ion bulk velocity and its streamlines, as well as the electric field magnitude with its field lines, respectively. Figure 3 has the same layout as Figure 2, but it shows the density and bulk velocity for the solar wind protons and the bottom row presents the magnetic field.
In all plots shown, the magnetopause boundary calculated from the Shue model (Shue et al., 1997) is represented by a dashed line. The two parameters characterizing the Shue magnetopause boundary are the subsolar standoff distance R SS and the flaring parameter α. To maximize the size of the region that can be populated by exospheric ions, we apply a low upstream solar wind dynamic pressure of 2.1 nPa, corresponding to a stand-off distance of R SS = (1.7 ± 0.1) R M (Winslow et al., 2013). The associated flaring parameter is α = 0.55, representing an open magnetopause. In our model run without an exosphere, we identify the subsolar location of the magnetopause to R SS = 1.6 R M (where u x = 0 km/s, not shown), which is within the range of values derived by Winslow et al. (2013). The subsolar point is located at the same position in all five cases, since the Shue model does not consider the counter pressure associated with an ionosphere. Within the period of a single Dungey cycle (about a minute, see Slavin et al., 2009Slavin et al., , 2010, the standard exosphere model given by Equation 1 inserts sodium ions with densities of 1 and 10 2 cm −3 near the magnetopause boundary at the subsolar point and above the cusps, respectively. These number densities are comparable to the upstream solar wind number density of 10 cm −3 . Therefore, the presence of the sodium ions influences the solar wind-magnetosphere interaction in ways we present in the next paragraphs.

Global Sodium Ion Structure
In Figures 2b-2d, the sodium ion density is fragmented into three ray-like structures, two of which develop north of the planet. The uppermost density enhancement with a maximum density of 0.1-0.5 cm −3 is associated with pickup cycloids of sodium that has been ionized outside of the magnetosphere. In perfectly uniform fields, pickup sodium ions would achieve a maximum velocity of 2u 0 = 1,000 km/s and a gyroradius of r g = 2.5 R M (Simon et al., 2007). However, as the majority of sodium ions within the northernmost ray are picked up within the magnetosheath where the solar wind protons have been decelerated, the local convective field is decreased. This, in turn, accelerates the sodium pickup ions to lower maximum velocities of only 550 km/s; see Figure 2g. Therefore, the width and height of the cycloidal arcs in the northernmost ray are 14 and 4 R M , respectively, and are slightly smaller than the respective values of 2πr g = 15.7 R M and 2r g = 5 R M in the uniform upstream fields (Simon et al., 2007). The more extended exospheres in Figures 2c and 2d inject sodium ions into stronger electric fields at the "flanks" of the magnetosheath (see Figures 2i-2n) so that the maximum velocity inside the northernmost arc reaches 750 km/s.
The other two rays in Figure 2b originate from the high-density cusp regions where sodium ions follow the draped northern and southern polar field lines. Consequently, the sodium ions accumulate inside the magnetosphere along the magnetopause boundary and envelop the tail lobes in the high latitudes. Similar to Earth, we identify these regions as the northern and southern plasma mantle regions (Frank, 1985;Fuselier et al., 2019;Pilipp & Morfill, 1978;Rosenbauer et al., 1975). These two rays reach their largest |z| values in narrow bands between z = 3.5-4 R M and z = −2.5-3 R M , respectively; that is, their thickness amounts to 0.5 R M , consistent with observations . The difference in the z offset of the rays is mainly due to the dipole offset and the large-scale asymmetries in the magnetosphere caused by the IMF inclination against the planetary dipole.
Before the MESSENGER era, knowledge about Mercury's dipole had been obtained only from Mariner 10 observations. However, the dipole strength had been overestimated by a few hundred nT as well as no offset of the dipole could be obtained. Under these assumptions, Delcourt et al. (2003) conducted simulations of test particle trajectories of sodium ions from a spherical exosphere. The simulation domain encompassed a simplified magnetosphere, consisting of the superposition of the dipole field and a Harris sheet model for the magnetotail. These sodium ion trajectories are focused into symmetric northern and southern rays, similar to our results in Figure 2b. However, the rays in their model converge into the equatorial plasma sheet at wakeside distances of x = 4 R M , while our modeled rays at this x value are nearly aligned with the magnetopause boundary and are just beginning to converge toward the plasma sheet.
Even though the exospheric profile in this work is symmetric between Mercury's northern and southern hemispheres, the modeled footpoint region of the southern ray appears to be more focused than its northern counterpart (in Figures 2b-2d). This is due to different trajectories of the ion guiding centers in both hemispheres. Most of the northern cusp ions are able to leave the planet's vicinity along the northern dipole field lines. In contrast, the pickup of sodium ions generated in the −x and −z quadrants is interrupted by the planet itself.
To illustrate the guiding center motion of the sodium ions in the southern hemisphere, we highlight three streamlines that are each 100 km apart, colored violet, red, and yellow, respectively. These sodium ions are generated in the −x and −z quadrant of Figure 2g, in the outermost region of our exosphere model (see Figure 1). At first, these streamlines are parallel and close together as the pickup direction of the associated sodium ions is the same inside the undisturbed solar wind. Important to note is that every ion generated along the same streamline follows the same path. However, after crossing the bow shock, the paths are deflected in different directions. The violet streamline turns around the planet toward the north into the magnetosheath. The sodium ions following the red and yellow streamlines are deflected toward Mercury's south pole. The red streamline shows how sodium ions precipitate onto the south-polar surface and may contribute to local surface erosion . The sodium ions that do not move toward the surface, denoted by the yellow streamline, are decelerated by the gradient of the plasma pressure within the southern cusp (Winslow et al., 2017). Therefore, the local sodium ion density is increased, thus making the southern cusp population slightly more dense than the northern counterpart. Subsequently, these sodium ions leave the southern cusp along the draped southern field lines. In summary, the asymmetry in the sodium ion densities near the cusps is partially generated by the influence of the solar wind on the magnetosphere and therefore could not be obtained by Delcourt et al. (2003). Furthermore, because the electric field of the upstream solar wind points northward, only the southern cusp is open for precipitation of pickup sodium ions to the surface (see red streamline in Panel 2g).
Ionospheric ions that propagate against the direction of the ambient electric field, that is, the southern rays in our results, have been observed during multiple Cassini flybys of Titan (Edberg et al., 2011;Modolo et al., 2007;Szego et al., 2007). Edberg et al. (2011) argued that this motion could stem from, for example, ambipolar electric fields or ion motion parallel to the local magnetic field lines. Indeed, in the case of Mercury, ion motion along the draped field lines in the southern plasma mantle is the result of local magnetic and particle pressure.
Following the sodium ions within the northern and southern plasma mantle into the downstream tail lobes in Figures 2b-2d, it appears that the rays cease at x values of about 6-8R M ; that is, the sodium ion density decreases to values below 1 cm −3 , indicated by the yellow regions. However, these rays are actually slightly deflected out of the xz-plane and continue to propagate mostly parallel to it at distances below values |y| < 2 R M . The motion of sodium ions out of the xz-plane in Figures 2b-2d is therefore a result of a nonvanishing ycomponent of the ambient electric field.
Within the 50-fold exosphere in Figure 2d, the maximum |z| values of the rays equalize to about z = ±2 R M , losing the northward offset that the rays displayed within less dense exospheres. Thus, the influence of the dipole offset is diminished with a thicker exosphere. This is the result of sodium densities that are comparable to the solar wind extending farther into the upstream solar wind than the region dominated by the planet's dipole. The distance where the dipole strength is comparable to the IMF strength is about 2 R M from the center of Mercury. In contrast, the 50-fold exosphere extends to a distance of 2.1 and 2.4 R M at the subsolar point and above the cusps, respectively, a farther distance of about 0.6 R M from the magnetopause given by the Shue model. Therefore, with a dense enough exosphere, Mercury would interact with the solar wind mainly through its ionosphere rather than through its magnetic field, that is, the interaction is more comet/Venus/Mars like in nature. This comet-like behavior is most visible in the extreme case of a 500-fold more dense exosphere in the rightmost column of Figure 2 where the exosphere pushes the bow shock upstream to x = −4.5 R M .
In a recent hybrid study, Egan et al. (2019) investigated how the solar wind interacts with a Mars-sized planet, its ionosphere and the resulting ion escape. This planet has been assigned a weak dipole field. The strength of which was successively increased between consecutive simulation runs. At small dipole strengths, the planet's intrinsic magnetic field is shielded by the farther extending ionosphere, consistent with our 500-fold exosphere run for Mercury. By subsequent increase of the dipole strength, the resulting magnetopause position exceeds the standoff distance caused by the interaction with the ionosphere alone. When an ICME has pushed the magnetopause below the surface, FIPS observed sodium ion densities of 3-4 cm −3 along the equatorial dayside segment of MESSENGER's orbit (Winslow et al., 2020). These authors argued that the increased impact of the solar wind particles onto Mercury's dayside surface leads to increased sputtering rates. Such an increase in the exospheric base density would be in agreement with the modeled sodium ion densities in the 50-fold exosphere run, where the values in the dayside equatorial region amount to 1-10 cm −3 .
Measurements of the nightside plasma sheet (x < 3 R M ) obtained sodium ion densities of about 0.1-1 cm −3 , with an average value of 0.66 cm −3 . The modeled sodium ion densities in the plasma sheet for the fivefold, 50-fold, and 500-fold increased exosphere runs are in agreement with these measurements as they yield maximum sodium ion densities in the sheet of 0.1, 0.5, and 1 cm −3 , as can be seen in Figures 2c-2e, respectively. For the standard exosphere case in Figure 2b, however, the densities found in the plasma sheet of about 0.01 cm −3 are about an order of magnitude lower than observations. This could be due to an underestimation of the nightside exosphere in the model of Gamborino et al. (2019) due to lack of nightside observations by MESSENGER or the underestimation of the nightside ionization frequency in our model. Jasinski et al. (2017) analyzed MESSENGER transits through the southern plasma mantle in the 1 < x < 4 R M segment. These authors showed that the average sodium ion densities amount to 0.004 cm −3 , with singular observations of up to 0.05 cm −3 . In the 1 < x < 4 R M segments of the southern plasma mantle of Panels 2b-2d, sodium ion densities reach values of 0.001-1, 0.005-5, and 0.01-20 cm −3 in the standard, fivefold, and 50-fold increased exosphere runs, respectively. The sodium ion densities in the northern plasma mantle in Panels 2b-2d do not differ significantly from their southern counterparts. Therefore, the standard and fivefold exosphere agree with the average values of the sodium ion density observed in the southern plasma mantle.

Journal of Geophysical Research: Space Physics
The E×B drift inside of the downstream magnetopause boundary is directed toward the nightside plasma sheet. Depending on the magnitude of the sodium ion's velocity along the magnetic field lines, slow particles are diverted into the plasma sheet at distances closer to the planet than faster particles, a significant portion of which is lost downtail . Jasinski et al. (2017) showed that sodium ions with low velocity divert into the plasma sheet around x = 3 R M . The corresponding flux of sodium ions into the plasma sheet has been derived to a value of about 6·10 7 cm −2 s −1 .
Multiplying the density distributions and velocity fields of each column in Figure 2, we can estimate the flux j Na =nu z of sodium ions that divert into the plasma sheet between 4 < x < 8 R M . Originating from the southern plasma mantle, we find sodium ion fluxes into the plasma sheet of about 1·10 7 , 3·10 7 , and 8·10 7 cm −2 s −1 for the standard, fivefold, and 50-fold exosphere, respectively. Therefore, the fivefold and 50-fold exospheres lead to sodium ion fluxes comparable to observations. Our results indicate that the fluxes originating from the northern plasma mantle do not differ from the southern fluxes by more than a factor of 2. Poh et al. (2017) showed that protons only contribute a third of the necessary pressure to the stress balance within the plasma sheet. These authors proposed that the remaining portion might be contributed by heavy ions. Indeed, our results indicate that an exosphere with a sufficient density might provide the needed amount of sodium ions in the plasma sheet through sodium ion fluxes from the plasma mantle.
In total, the model runs using the fivefold and 50-fold exosphere agree with most MESSENGER observations of sodium ion densities. The orbits of the upcoming BepiColombo mission will provide equal coverage of both hemispheres of Mercury at low altitudes (Benkhoff et al., 2010). This will enable a more accurate estimate which of our exosphere models represents the actual exosphere of Mercury.
Sodium ions are inserted with an initial velocity of v i = 0  and diminish the local con- The nearly stationary sodium ions shield themselves from the electric field of the IMF, analogous to a Faraday cage (Bagdonat, 2005). This leads to a weakened pickup process and thus to an even further accumulation of sodium ions. This feedback successively increases the shielding of the sodium ions from the electric field of the solar wind. The pickup process is limited to the erosion of the outer edges of the high-density regions . Within the northern cusp region in Figure 2g, pickup mainly takes place at the outer edges of sodium densities below values of 1 cm −3 .
As the sodium ion density is larger in the southern cusp, its electric field is always significantly weaker than in the northern counterpart, see Figures 2l-2n. However, the electric field within the low-latitude magnetosphere in Figure 2k is locally increased to 30 mV/m in three distinct regions. This corresponds to about three times the value in the upstream solar wind and is depicted by red colors in Figure 2k. Within the nightside, the high electric field is the result of the ring current and the neutral current sheet within the magnetotail. The crescent-shaped region of enhanced electric field structure on the dayside is located along the outermost closed magnetic field lines and thus indicates the magnetopause location. The inner structure on the dayside, just 0.2 R M above the surface, stems from the ring current. Generated mainly by the gradient and curvature drifts within the dipole field, drifting solar wind ions generate their own electric field. With the introduction of denser exospheres in Figures 2l-2n, the shielding effect of the sodium ions and associated electric field decrease the ring current's electric field. Furthermore, the protons within the ring current are gradually dominated with sodium ions. This is because the ram pressure of the upstream solar wind is low and therefore the magnetosphere is able to encompass the larger gyroradii of sodium ions for the formation of an equatorial "sodium ring" around the planet (Yagi et al., 2017). Therefore, when sodium ion densities increase, the associated decrease of the electric field will ultimately overcompensate the field associated with the ring current. Indeed, this is visible in Figures 2m and 2n where sodium ion densities larger than 10 cm −3 , that is, comparable to the upstream solar wind density, reach into the ring current region. This is consistent with the results by Paral et al. (2010), where sodium ions, derived from a test particle model, populate the equatorial ring current at densities comparable to the applied upstream solar wind density.

Solar Wind
In previous simulation works about the solar wind interaction with Mercury, a possible influence of the exosphere has not been considered, as the sodium ion densities measured by FIPS were low enough to be regarded as test particles (Dong et

Journal of Geophysical Research: Space Physics
near-surface sodium ion densities that are up to two orders of magnitude higher than those derived from FIPS measurements, which are comparable to ion densities of the solar wind and, thus, likely to affect the global magnetosphere. To compare our model output to these preceding studies, we show Mercury's magnetosphere without exospheric pickup in the leftmost column of Figure 3. The solar wind in Figure  3a crosses the bow shock at a subsolar surface altitude of 0.9 R M . There, the solar wind density is increased to values of 35 cm −3 , which corresponds to a factor of 3.5 of the upstream solar wind density. A corresponding jump by a factor of 3.5 is seen in the solar wind velocity in Figure 3f and magnetic field in Figure 3k. However, due to the low solar wind ram pressure in our model, the jump is slightly weaker than the maximum factor of 4 expected from the Rankine-Hugoniot relations. For the upstream parameters considered here, Mercury's solar wind interaction is therefore "weaker" compared to a magnetosphere resulting from average upstream conditions. Through FIPS measurements within the wakeside plasma sheet of selected MESSENGER orbits (in 2011-2012) when Mercury was near its perihelion, that is, under influence of high solar wind ram pressure, Gershman et al. (2014) obtained solar wind densities in a range of about 1-10 cm −3 with an average solar wind density of 7.66 cm −3 . Within the modeled plasma sheet (see Figure 3a), the solar wind density is in the range of about 0.05-0.5 cm −3 , a value smaller than the observed value by Gershman et al. (2014) by one order of magnitude. However, as we employ lower than average solar wind ram pressure, fewer solar wind ions are able to penetrate into Mercury's magnetosphere and populate the plasma sheet. Indeed, for quiet solar wind conditions when MESSENGER crossed the nightside plasma sheet (in 2013-2015), Dewey et al. (2018) showed that the average solar wind density roughly halved to around 3.1 cm −3 , to which our modeled solar wind densities are in a reasonable range of values. It is therefore likely that the low solar wind ion density in the plasma sheet of our model is indeed a result of undisturbed, low ram pressure conditions. Two major features within the polar cusps can be discussed with our results, that is, the solar wind ion density in the immediate vicinity of the surface and the maximum solar wind density at higher altitudes covered by MESSENGER. For the latter, the solar wind density inside the cusps consists of ions traveling into the cusp and ions reemerging due to the magnetic mirror force. In the case of a planet-centered dipole, the cusps should be of the same size (Delcourt et al., 2003) and populated by equal solar wind densities. To illustrate the ion dynamics near the cusps in Mercury's offset dipole field, we highlight two streamlines of the solar wind bulk velocity at 0.02 R M above and below the magnetic equator in Figure 3f. Both streamlines are deflected at the magnetopause boundary into the northern and southern hemisphere, respectively. After reaching about halfway into the northern cusp, the northern, yellow streamline is sharply diverted northward, while the southern, red streamline is only marginally deflected southward. This is due to the offset dipole, in which the mirror points for a larger fraction of the solar wind ions are beneath the southern surface, that is, a significant fraction of the solar wind ions in the southern cusp to impact onto the surface than in the northern counterpart, where ions are able to reemerge instead. The resulting bulk velocity in the northern cusps is therefore stronger affected than in the southern cusp; that is, the yellow streamline is stronger deflected than the red streamline.
In other words, it is expected that the maximum solar wind density is of lower value in the southern cusp than in the northern counterpart (Exner et al., 2018;Fatemi et al., 2018;Trávníček et al., 2010). Indeed, in Figure 3a, the maximum solar wind density within the northern cusp ranges between 60 and 70 cm −3 at surface altitudes of 0.5 R M , while the maximum solar wind density within the southern cusp is in the range of about 50-55 cm −3 , which corresponds to values of 10-30% smaller than the northern maximum. Within the northern cusp, the high-density region is located at surface altitudes of about 0.4-0.6 R M and forms an elongated structure that approximately extends over 1.5 and 0.5 R M in length and width, respectively. These northern cusp dimensions agree with reported cusp crossings of up to 0.6 R M , as MESSENGER's terminator orbits had been nearly perpendicular to this structure. For the average solar wind density observed in the northern cusp, these authors reported a lower limit of at least 10 cm −3 . Indeed, when moving from left to right in the first row of Figure 3, the maximum solar wind density within the northern cusp is always larger than 40 cm −3 and therefore agrees with the lower limit given by Raines et al. (2014). The 500-fold exosphere case constitutes the only exception with a value of 5 cm −3 due to the dense exosphere limiting the access into the cusp.
Regarding the solar wind densities at the polar surfaces, Winslow et al. (2012) argued that the solar wind density at the surface in the southern cusp should be a factor of 4 higher than the solar wind density in the northern counterpart. Indeed, our modeled solar wind densities in these two surface regions are about 5 and 25 cm −3 in the northern and southern cusps of Figure 3a, corresponding to an increase by factor of 5.
By investigating multiple dipolarization events in Mercury's wakeside, Dewey et al. (2018) found fast plasma flows in Mercury's vicinity that move duskward with an average velocity magnitude of 300 km/s. Within surface altitudes analyzed by these authors, our modeled solar wind has a bulk velocity in the range of 200-700 km/s in Figure 3f, encompassing the observed value. With increasing distance from the planet and outside the analysis region of Dewey et al. (2018), however, the solar wind protons inside the nightside plasma sheet are accelerated to larger velocities of about 2u 0 = 1,000 km/s. Furthermore, from x = 1.2 R M , the plasma sheet in Figure 3f exhibits a warped shape; that is, its z position gradually increases from z = 0.2 R M to a z value of 1.5 R M which is an effect of the twisting of the tail lobes (the decreasing influence of the planetary dipole with distance).
The only discernible effect of the standard sodium exosphere on the magnetosphere is seen within the magnetotail in the second column of Figure 3. The plasma sheet has become parallel to the x-axis for additional 0.8 R M to x = 2 R M until it exhibits the previously warped shape again. This is the result of dominating sodium ions against the reduced magnetic field in the neutral sheet which locally symmetrizes the plasma sheet. Indeed, with the fivefold and 50-fold denser exospheres and thus larger sodium ion densities within the tail, the plasma sheet remains parallel to the x-axis to larger distances of up to x = 4 R M in Figures 3h and 3i (Simon & Motschmann, 2009).
In the cusps between the model runs with the standard and 50-fold denser exosphere, the maximum solar wind densities within the northern cusps of Figures 3b and 3d have roughly halved from 70 to about 40 cm −3 . The sodium ion density in the fivefold dense exosphere run (see Figure 3c) amounts to competing values to the solar wind within the dayside magnetosphere. This leads to a nearly disappearing solar wind 10.1029/2019JA027691

Journal of Geophysical Research: Space Physics
population, while the sodium ions dominate the equatorial dayside magnetosphere. Therefore, not only do the larger sodium ion densities in the cusp regions shield themselves from the pickup by the ambient electric field, but also the pressure gradients at the outer flanks of these sodium populations prevent the solar wind protons from entering the cusps while the sodium ion population takes over as the major carrier of the ring current. As previously stated, from the fivefold dense exosphere on, the solar wind mainly interacts with the inflated ionosphere rather than through the planetary magnetic field (see Figures 2c and 2d). This is also evident by the subsolar point of the magnetopause gradually moving toward the planetary equator (i.e., z = 0 R M ) in Figures 3c and 3d, as denoted by a white cross. The symmetrization is a result of the exosphere profile in our model that is symmetric with respect to the planetary equator. Therefore, combining all previous effects, in regions where sodium ion densities compete and dominate against the . Modeled magnetic field components and magnitude along a simplified MPO trajectory from the model runs without a sodium exosphere, and with the standard, fivefold, and 50-fold sodium exosphere included. The bottom panel shows the spacecrafts position along the trajectory, given in planetary radii in MASO coordinates; that is, the trajectory is confined to the xz-plane and starts from its apogee in the nightside at the geographic equator. The closest approach is indicated by the dashed black vertical line. The orange colored boxes in the background represent the regions when the trajectory crosses the plasma sheet, northern cusp and southern cusp from the model run without the exosphere, respectively.

10.1029/2019JA027691
Journal of Geophysical Research: Space Physics solar wind ions, a test-particle approach for the sodium ions should not be valid anymore as the global magnetosphere is affected.
In the previous section, we found that the northern and southern plasma mantle are filled by sodium ions originating from the polar regions. Panel 3a allows to correlate the downtail magnetopause boundary of the Shue model with the regions where the solar wind density is decreased to values of about 5 cm −3 , depicted in orange. Within the 1 < x < 4 R M segment of the southern mantle, the modeled solar wind density amounts to values between 0.01 and 4 cm −3 . The observed values of 0.5-1.5 cm −3 (Jasinski et al., 2017) are therefore well within the range of our results.
To assess the influence of the exosphere on the magnetic field more quantitatively, we present the magnetic field of the model runs without an exosphere, with the standard, fivefold, and 50-fold exosphere along a simplified trajectory of the Mercury Planetary Orbiter (MPO, the more closely orbiting spacecraft of the BepiColombo mission) in the upper four panels of Figure 4. The orbit of the spacecraft is assumed to be confined to the meridional xz-plane with its apogee in the nightside region, where the spacecraft crosses the geographic equator plane in northward direction (see bottom panel of Figure 4 and white orbits in Panels 3k-3n). For simplicity, we present the fields and the spacecraft's position in terms of percentage of the trajectory, starting from its apogee. Thus, the closest approach is located at 50% of the orbit in the dayside magnetosphere, indicated by a dashed black vertical line in Figure 4. The orange boxes represent the segments of the plasma sheet, northern and southern cusp crossings, as estimated from the model run without an exosphere from Figure 3a.
In the previous section, we investigated how the increasing sodium densities in the tail are able to straighten out the plasma sheet. The same effect is seen within the plasma sheet region of Figure 4, where the neutral sheet crossing is indicated by the change of the sign of the B x -component. In the absence of an exosphere, this crossing occurs at about 4% of the orbit, at a z location of about z = 0.4 R M , while the crossing has moved to the magnetic equator at about 2% of the orbit, that is, z = 0.2 R M , in the model run with the 50-fold exosphere. The sign of the B x -component changes from positive to negative and the minimum total magnetic field of the different model runs is about 45-50 nT. Consequently, the location of the X-line must be farther downtail than the apogee of the MPO orbit. Indeed, the X-line location in Panels 3k-3n is located between 3 < x < 7 R M . This is in agreement with Poh et al. (2017), who reported an average downtail location for the X-line at a x = 3 R M , which is about 1.5 R M further downtail than the apogee of the MPO trajectory that we consider here. Zhang et al. (2020) investigated the different z locations of all MESSENGER crossings of the neutral sheet and compared these to the average location of the plasma sheet. These authors showed that Mercury's tail exhibits a flapping motion that depends on the plasma density within the plasma sheet (in which only solar wind particles were considered) and the variability of the upstream solar wind conditions. However, our results (considering steady upstream conditions) indicate that the location of the neutral sheet may also depend on the exospheric density.
At closest approach, the magnetic field magnitude in the model runs with the standard, fivefold, and 50-fold exospheres are about 10, 17, and 20 nT weaker than without an exosphere. The compression of the planetary magnetic field has therefore weakened by up to 15% of the local magnetic field magnitude from the model run without an exosphere. This indicates an upstream displacement of the subsolar magnetopause location, as seen in Panels 3k-3n.
As the magnetosphere is increasingly inflated, the magnetospheric cusps are affected by the increasing exospheric density as well. Along this orbit of MPO, both cusp positions are affected, albeit the southern cusp shows more significant changes. The northern cusp crossings have moved poleward from around 35-41% of the orbit without an exosphere to about 34-40% and 32-38%, when the respective standard and 50-fold exospheres are included. The respective magnetic field magnitudes have decreased by 10 and 20 nT, which corresponds to about 3% and 7% of the local magnetic field magnitude when no exosphere is included. The southern cusp region, however, initially present at about 52-65% of the orbit when no exosphere is considered, has moved farther toward the southern pole to about 53-66% and 57-73% of the orbit with the respective standard and 50-fold exospheres included. The corresponding differences in the magnetic field magnitude are about 7 and 50 nT, that is, about 7% and 35% of the local magnetic field magnitude of the 10.1029/2019JA027691 Journal of Geophysical Research: Space Physics model run without an exosphere. Thus, both cusps have moved slightly poleward with a progressively dense sodium exosphere. Winslow et al. (2017) showed how the magnetospheric cusps shift poleward with decreasing upstream dynamic pressure. Our results reveal that by including a dense sodium exosphere (which increases the counter pressure inside the magnetosphere), the poleward motion can also be attributed to the existence of a significant sodium exosphere. Therefore, the variability of the subsolar magnetopause location of 0.2 R M by Winslow et al. (2013) may partly be contributed by the influence of a dense exosphere at low upstream dynamic pressure conditions. However, we would like to emphasize that an increase in the exospheric sodium ion density by a factor of 50-500 (as considered in the last two model runs) would likely only be a transient event caused by extreme SP, that is, during an ICME impact. Therefore, when considering multiple spacecraft orbits covering an extended period of time, such a singular event is not expected to change the average positions of the magnetospheric boundary layers beyond the variability already found by Winslow et al. (2013).

Influence of Sodium Ions on the Currents in Mercury's Magnetosphere
To investigate how the sodium exosphere affects the current structure within Mercury's magnetosphere, we show the modeled currents in Figure 5. The first and second columns of Figure 5 depict the current density for the runs without any sodium ions and with the 50-fold dense exosphere. The last column shows the difference between the current densities in these two scenarios; that is, quantities of the model run without the exosphere are subtracted from the quantities of the 50-fold dense exosphere run. The current density and the components of the current perpendicular to the local magnetic field in the xz-plane are shown in the first and second rows of Figure 5, respectively. We show the FACs within the yz-plane in the bottom row of Figure 5. To facilitate a direct comparison to the results of Anderson et al. (2014Anderson et al. ( , 2018 and Dong et al. (2019), we flip the sign of the FACs as follows: Thus, red areas/positive values in Panels 5g-5i depict FACs that flow antiparallel to the local magnetic field, while blue areas/negative values represent FACs in the same direction as the magnetic field. Thin black lines in the first and second columns of Figure 5 are used to show the magnetic field lines.
Three regions of high current density within the magnetosphere in the upper and middle rows of Figure 5are readily recognized, that is, the bow shock and magnetopause boundaries, and the nightside current sheet. While the magnitudes of the first two currents do not change by more than 10% when the sodium exosphere is included, the additional pressure of the sodium ions in Figure 5b pushes their subsolar positions farther away from the planet. In Figure 5c, the displacements are clearly discernible, that is, the two upstream current layers have moved from subsolar surface distances of about 0.9-1.4 R M and 0.5-0.8 R M for the bow shock and magnetopause positions, respectively. The thickness of these two current layers seem to be unaffected by the displacement. Winslow et al. (2013) determined that the variability of the subsolar magnetopause location can vary by about 0.2 R M at low upstream ram pressures of 2 nPa. In their investigation of induction processes due to time-varying solar wind ram pressure, Jia et al. (2015) found that the subsolar magnetopause position moves by about 0.1 R M when induction is considered. Heyner et al. (2016) showed that the magnetopause position can move toward the planet by a maximum of 0.055 R M due to magnetopause erosion by magnetic reconnection. Our results indicate that the presence of significant amounts of sodium ions changes the magnetopause position by up to 0.3 R M . The analysis of the magnetopause locations by Winslow et al. (2013) did not take into account a possible influence of a sodium exosphere. Thus, it may be likely that in low upstream pressure conditions, the presence of a dense sodium exosphere may contribute to the variability of the magnetopause location. Simultaneous observations of the upstream solar wind conditions and displacement of the magnetopause from its expected position by spacecrafts of BepiColombo could, therefore, be used as an additional proxy to determine the column density of the sodium exosphere.
The dayside ring current above the equatorial surface in Figure 5a (z = 0.2 R M , x = −1.1 R M ) exhibits a magnitude of about 60 nA/m 2 , which is depicted in blue. This value has increased by a factor of 2.5 to 150 nA/m 2 , depicted in red in Figure 5b, due to local sodium ion densities being larger than the solar wind density by at least four orders of magnitude. Furthermore, the bulk velocity of solar wind ions within the ring current 10.1029/2019JA027691

Journal of Geophysical Research: Space Physics
region is nearly zero in Figure 3i, compared a nonvanishing bulk velocity of the sodium ions of about 200 km/s in Figure 2i, and therefore, sodium ions have taken over as the major carrier of the ring current.
The inner edge of the nightside current sheet in Panel 5a is located close to the nightside surface. However, when the 50-fold exosphere in Panel 5b is considered, the inner edge of the current sheet has moved to higher altitudes at x = 0.3 R M . The latter is in agreement with observations of the inner edge of Mercury's current sheet located at 0.22 R M (Poh et al., 2017) which indicates that a significant sodium ion density is needed in the nightside current sheet to establish local stress balance. By analyzing all MESSENGER magnetometer observations for deviations from the average magnetospheric model in the polar regions , Anderson et al. (2018) investigated the dependency of the FACs on the disturbance index of Mercury's magnetosphere . These authors showed that a low disturbance index (0-20%) is associated with a small magnitude of R1-FACs in the range of about 80-120 nA/m 2 , compared to a range of about140-300 nA/m 2 when the disturbance index is high. In a multifluid model of Mercury's magnetosphere during the MESSENGER M2 flyby (where the upstream dynamic pressure was low, i.e., about 11 nPa), Dong et al. (2019) found that the R1-FACs in the northern hemisphere amount to magnitudes of 150 and 115 nA/m 2 for the descending/planetward and ascending/antiplanetward current directions, respectively.
As we employ even lower upstream dynamic pressure than Dong et al. (2019) (about 2 nPa), we expect a low magnitude for the R1-FACs. Indeed, the magnitude of ascending R1-FACs in the duskside of the northern polar cap (red arc in Figure 5g) amounts to 60-80 nA/m 2 , which is in reasonable agreement with the previous values of Anderson et al. (2018) and Dong et al. (2019). At large, the R1-FAC system in Figure 5g, indicated by black arrows, follows the shape of the R1-FAC system sketched by Anderson et al. (2014) but depicts asymmetries in detail in between the duskside and dawnside hemispheres. This is because the IMF exhibits a large B y -component and therefore the R1-FAC system is slightly rotated around the x-axis. The R1-FAC latitudes in the analysis of Anderson et al. (2014) showed no significant asymmetry between the dawnside and duskside hemispheres as temporary features of the IMF's B y -component average out; that is, the net value of the IMF's B y component vanishes over the time span covered by the MESSENGER mission. A strong current of 60-80 nA/m 2 is emerging from the equatorial surface in duskside (y = −1 R M , z = 0 R M , depicted in blue) and belongs to the southern R1-FAC system as the magnetic equator is offset to the north. No counterpart to this current is present at the dawnside surface region in the northern current system (indicated by the green ellipse). This is due to the southward directed current, in the dawnside magnetopause layer (light blue arc), reaches z valuesof −0.5 R M , and merges with the red arc of the opposing, northward directed current of the southern R1-FAC system. The combined currents then turn planetward and connect to the dawnside surface below the equator. In agreement with Anderson et al. (2018) and Dong et al. (2019), no current is noticeable below the latitudes of the R1-FACs, that is, no R2-FAC system seems to be present in the case of no exospherein Panel 5g.
When including the sodium exosphere, however, the R1-FAC system in Figure 5h appears to symmetrize between the duskside and dawnside hemispheres. The northern R1-FAC system at the dawnside does not merge with the southern system (as it has in Panel 5g) but turns northward to connect to the northern surface at medium latitudes, as depicted by the blue arc. Comparing the extent of the closed dipole field lines between Panels 5 g and 5h, the surface altitudes of the outermost field lines have increased from about 0.4-0.6 R M and 0.2-0.6 R M at the duskside and dawnsides, respectively. This is due to the increased magnetopause standoff distance when sodium ions are included, leading to a larger cross section of the magnetosphere in the yz-plane. Due to this inflation of the dawnside magnetosphere, the blue arc depicting the northern R1-FAC system is able to return to the northern surface. It seems that the descending R1 currents at the dawnside of the northern hemisphere and the ascending R1 currents at the duskside of the southern hemisphere (denoted by the blue arcs in Figure 5h) are restricted to a small range of medium latitudes and thereby increased to magnitudes of about 100 nA/m 2 . In contrast, the complementary R1 currents in the respective hemispheres (denoted by the red areas at the surface in both hemispheres) span much larger latitude ranges and therefore show lower magnitudes of about 40-60 nA/m 2 . This is consistent with to the findings of Dong et al. (2019), who showed that the maximum planetward R1 current is about 35 nA/m 2 stronger than the emerging current.
Due to the inflated magnetosphere when the sodium is included, new currents (depicted by green arrows) develop at low altitudes of up to 0.3 R M in the equatorial latitudes of both sides of Mercury in Panel 5h.
Flowing from low latitudes of the northern hemisphere to low latitudes of the southern hemisphere, these currents have a magnitude of about 20-30 nA/m 2 , which is about a factor of 4 weaker than the adjacent R1 currents. These weak currents are likely to be associated with the yet unobserved R2-FACs.
It has been shown that R2-FACs in Earth's inner magnetosphere are mostly driven by pressure gradients within the plasma sphere which extends about 3-4 R E , where R E is the radius of the Earth (Ganushkina et al., 2015;Tsyganenko, 2000). Siscoe et al. (1975) showed that a scaling factor of 10 with M i and P i being the magnetic moments and upstream solar wind dynamic pressures at Earth and Mercury, respectively, is suitable to scale the Earth's subsolar magnetopause standoff distance (≈11 R E ) down to the magnetopause standoff distance at Mercury, as obtained from Mariner 10 observations (≈1.6 R M ). In other words, these authors suggest that regarding the size of the magnetosphere, a characteristic length scale of 1 R M at Mercury would correspond to 6.9 R E at Earth. This implies that Mercury occupies such a large volume within its own magnetosphere that the Hermean plasma sphere would be located at 0.57 R M , that is, within the planet itself (Russell et al., 1988). However, MESSENGER observations showed a diminished magnetic moment of Mercury by a third of the value obtained by Mariner 10. Additionally, we employ an upstream dynamic pressure that is about a factor of 7 weaker than the average dynamic pressure observed at Mercury. Our results of Panel 5b show that the 50-fold exosphere leads to a displacement of the subsolar magnetopause by 0.3 R M farther upstream. Consequently, in our model run of the 50-fold exosphere, 1 R M at Mercury would correspond to 5.4 R E at Earth; that is, the relative size of the Hermean magnetosphere has increased by about 17%. The plasma sphere of Mercury would, therefore, still lie beneath the surface with 0.75 R M . However, Anderson et al. (2018) argued that the appearance of the ring current system is an indirect indicator for the presence of R2-FACs. For that matter, results of Yagi et al. (2017) and our results of the model run with the standard exosphere show that a partial ring current consisting mainly of sodium ions is present within the nightside, while the ring current fully encloses the planet in the 50-fold exosphere model run. Therefore, the modification of the magnetospheric structure by the sodium ions close to the surface cannot be captured by the scaling approach of Siscoe et al. (1975), and the dynamics of the sodium ions and their large gyroradii lead to the formation of a ring current close to Mercury's surface Yagi et al., 2017).
Furthermore, R2-FACs originating from the plasma sphere do not necessarily need to connect with the exo-ionosphere near the terminator plane as proposed by Anderson et al. (2018) but rather connect to the nightside region of the plasma sheet as a recent study of Liu et al. (2016) shows. These authors analyzed THEMIS data acquired between the plasma sphere and plasma sheet in Earth's magnetotail. For the first time, these authors demonstrated that R2-FACs can exist in downtail distances of 8-12 R E in the plasma sheet and then connect to the nightside ionosphere of Earth. Anderson et al. (2018) argued that MESSENGER did not observe the R2-FAC system as the spacecraft's trajectory (between 400 and 1,000 km) did not intercept these low-altitude regions. This implies that either R2-FACs exist at lower altitudes or within latitudes of ±30°from the magnetic equator . Indeed, MESSENGER descended to altitudes of about 0.3 R M above the equator and would barely graze the R2-FACs in our model. Even in the case of MESSENGER crossing the R2-FAC system, the magnetic field disturbances generated by the adjacent R1-FACs would likely obscure the magnetic field disturbances generated by the R2-FAC system which has a roughly 4 times smaller amplitude.
As an additional reason for the not observed R2-FACs by MESSENGER, we found in Panel 5g that even at very low solar wind ram pressure, when the magnetosphere is already inflated compared to the average case (Winslow et al., 2013), no R2-FACs exist within the yz-plane of our model. It seems that a significant amount of sodium ions is needed to exert an additional internal counter pressure to the solar wind to further inflate the magnetosphere. With the increased subsolar magnetopause position, the resulting magnetosphere is inflated enough to "host" R2-FACs. Therefore, in the MESSENGER era, the absence of observable R2-FACs may likewise be due to the sodium density likely being smaller than our 50-fold dense exosphere case.
Due to the expected absence of a significant ionosphere, Anderson et al. (2014) proposed a R1 closure current that flows through Mercury's mantle from the dawnside to the duskside (see Figure 4a in their paper). In our model run without an exosphere in Figure 5d, large magnitude perpendicular currents are located at the core-mantle boundary. The y-component (not shown) of these currents is negative in the polar regions, which are indicated by orange ellipses. These duskward currents flow in a similar manner as presented by Dong et al. (2019). While these authors obtain a magnitude of about 500 nA/m 2 , the magnitude resulting from our model is about 400 nA/m 2 , due to lower solar wind ram pressure in our model.
However, Anderson et al. (2014) acknowledged that a dense enough ionosphere would lead to a significant portion of the R1 closure currents be carried by ionospheric Pedersen and Hall currents instead. As can be seen in Panel 5f, when deriving the difference of the current magnitudes of the 50-fold exosphere run to the run without sodium included, the R1 closure current at the polar core-mantle boundary has decreased by about 100 and 40 nA/m 2 in the northern and southern hemisphere, respectively. With the sodium exosphere included in Panel 5e, strong perpendicular currents of magnitudes around 80-100 nA/m 2 have emerged at medium latitudes on the dayside, as indicated by yellow ellipses. These regions coincide with the highest sodium ion densities at near-surface altitudes in the dayside, as seen in Figure 2d. Comparing these regions with our exosphere model in Figure 1, the associated ion populations are predominantly due to PSD processes. Therefore, a sufficiently dense cloud of planetary ions is able to "intercept" a significant amount of the R1 closure currents before they penetrate into the planetary interior. Due to the longitude-dependent exosphere, however, these closure currents are not above the poles as Anderson et al. (2014) proposed (as that is the shortest path for them to reach from the duskside to dawnside R1-FAC) but tied to regions of sufficiently high sodium ion densities on the dayside. Furthermore, a strong perpendicular current of about 160 nA/m 2 is identified at the south pole in Figure 2d (see red area in the southern ellipse) that has no counterpart in the northern hemisphere. This is due to the dipole offset and the associated northward shift of the southern Chapman-Ferraro current system closer to the surface.

Summary
In this study, we investigate how a dense sodium exosphere would affect the magnetic fields and current systems in the small magnetosphere of Mercury. Important open questions are the contributions of the exosphere to the closure of the R1-FACs and the conditions for the existence of R2-FACs.
Direct FIPS measurements of sodium ions show densities of 2 cm −3 in the northern magnetospheric cusp and reduced values of 5.1 · 10 −3 cm −3 in the dayside magnetosphere . However, in a recent study of field-line resonances, James et al. (2019) indirectly derived a much higher sodium ion density of about 22 cm −3 above the dayside surface, that is, multiple orders of magnitude above the value from FIPS.
As these values are comparable to the upstream solar wind density, it is conceivable that Mercury's sodium exosphere can affect the global magnetospheric structure. Especially the closure of R1-FACs in Mercury's magnetosphere is not fully understood, and sodium ions might play a role to provide necessary Pedersen and Hall conductivities for closure outside of Mercury's interior (Anderson et al., , 2018. As sodium ion gyroradii are large compared to the radius of Mercury, we employ the AIKEF hybrid model (kinetic ions, massless electron fluid) that has successfully been applied to the Hermean magnetosphere by Müller et al. (2012) and Exner et al. (2018) and is able to reproduce key features of the magnetic field perturbations observed by MESSENGER. In this study, we expand our model to include a sodium exosphere. The sodium ions are photoionized and, in contrast to any preceding model, may generate currents that affect the electromagnetic fields. We represent our sodium exosphere by a superposition of Gaussian bells, taking into account the four major processes responsible for exospheric genesis at Mercury: TD, PSD, SP, and micro-meteoroid impacts. Scale heights, scale angles, and latitude dependencies of these processes have been empirically fitted to the results of a Monte-Carlo model of the sodium exosphere by Gamborino et al. (2019). Our exosphere model is symmetric with respect to the planetary equator and exhibits the largest sodium densities in the medium and equatorial latitudes of the dayside due to TD and PSD processes, respectively. To understand the effect of a dense exosphere on Mercury's magnetosphere, we treat the density of the sodium exosphere as a free parameter and conduct model runs with multiples of the surface density derived from Gamborino et al. (2019) of 0, 1, 5, 50, and 500, respectively.
In the model run without an exosphere, the subsolar positions of the magnetospheric boundary layers as well as the solar wind densities of the northern cusp and nightside plasma sheet agree well with average values measured by MESSENGER. The magnitude and direction of the modeled R1-FACs are consistent with the currents obtained from MESSENGER MAG observations, while current closure is governed mainly by currents flowing along the core-mantle boundary (Anderson et al., 2018). However, due to the large y-component of the upstream magnetic field, the R1-FACs exhibit a dawn-dusk asymmetry within the terminator plane. Furthermore, no R2-FACs are discernible in the model results.
After including the sodium exosphere model of Gamborino et al. (2019), the magnetospheric boundary layers are not yet affected. Downstream of Mercury, the sodium ions are confined to three ray-like structures, two of which have their respective footpoint located in the cusp regions. The cusps also coincide with the regions of highest sodium densities from our exosphere model. Due to their large gyroradii, the pickup ions initially move northward. Therefore, a fraction of sodium ions that have been picked up upstream of the bow shock in the southern hemisphere travel into the southern cusp and increase the sodium density, resulting in an asymmetry between the sodium densities modeled for the cusps. Farther downstream, these rays reach the inner side of the magnetopause boundary and contribute to the enhanced densities within the northern and southern plasma mantle. The third ray consists of sodium ions that have been picked up in the northern high to medium latitudes within the magnetosheath. Sodium densities within the northern cusp and plasma sheet are in reasonable agreement with FIPS measurements.
The sodium ions leave the test particle regime when the exospheric density exceeds the values from Gamborino et al. (2019) by a factor of 5. In this case, sodium ion densities amount to values comparable to the solar wind density. Consequently, the sodium ions exert an additional counter pressure onto the magnetopause from inside the magnetosphere, which results in a slight displacement of the magnetospheric boundary layers away from the planet by about 0.1 R M . The sodium ion densities at the footpoints of the two polar rays have become so large that they start to shield themselves from the ambient electric field. This leads to a less effective pickup process and therefore, to a further buildup of the sodium ion density (Bagdonat, 2005). Therefore, at the altitudes covered by MESSENGER within the northern cusp, sodium ion densities are still in good agreement with MESSENGER observations, although the neutral gas density has been increased by a factor of 5. The fivefold denser exosphere leads to sodium ion densities comparable to values obtained by James et al. (2019) above the dayside surface. This also leads to the ring current being mainly carried by the sodium ions.
An exosphere model that is increased by a factor of 50 (compared to Gamborino et al., 2019) results in a significantly inflated magnetosphere that reveals several new features, compared to the two preceding scenarios. The exosphere reaches farther into the upstream solar wind than the planetary magnetic field. Therefore, Mercury's interaction with the solar wind is mainly governed by its ionosphere rather than by its magnetic field, that is, the interaction is more comet/Venus/Mars like in nature. Consequently, the magnetopause is significantly displaced by 0.3 R M away from the planet, a displacement of the same order of magnitude as caused by induction and erosion by reconnection processes Jia et al., 2015). Due to the exospheric model being symmetric with respect to the planetary equator, the subsolar point of the magnetopause is shifted southward, no longer affected by the dipole offset. The further decreased electric field leads to a stronger shielding effect of the sodium ions at low altitudes from the pickup. Only the outer edges of the high-density regions are eroded by the pickup process, which leads to sodium ion densities that are still consistent with observations in the altitudes covered by MESSENGER. Therefore, even a high neutral gas density may lead to unexpectedly low ion densities at higher altitudes, as the pickup process nearly ceases at low altitudes.
Furthermore, the closure current for the R1-FAC system in the 50-fold exosphere model run along the core-mantle boundary has decreased in magnitude. In turn, strong perpendicular currents form within the high-latitude regions of the dayside, coinciding with regions of increased sodium ion densities. Therefore, the R1-FACs are partially closed by the Pedersen and Hall currents of the sodium exosphere before they penetrate into Mercury's interior  as the Pedersen conductance is about a factor of 50 larger than the conductance of the mantle in the 50-fold exosphere run. In this setup, additional currents develop at equatorial latitudes at surface altitudes of about 0.3 R M . These currents have been identified as the yet unobserved R2-FACs. As MESSENGER's altitude ranged from 0.3 to 1.2 R M in this latitude range, the spacecraft would have barely grazed the R2-FAC system (Anderson et al., 2018). Moreover, as the magnitude of the modeled R2-FACs is a factor of 4 smaller than that of the adjacent R1-FACs, it is likely that the magnetic field disturbances generated by R2-FACs are obscured by the magnetic field disturbances generated by the R1-FACs.
Our results indicate that even at very low solar wind ram pressure, Mercury's magnetosphere does not exhibit a strong R2-FAC system. It seems that exospheric sodium ions of sufficient density are needed to exert a counter pressure against the upstream solar wind which inflates the Hermean magnetosphere significantly.
This enables R2-FACs at low altitudes which connect to the equatorial ring current system. Therefore, it is also possible that MESSENGER did not observe R2-FACs due to the sodium density not reaching values comparable to our dense exosphere.
Overall, our model results show that within a "reasonable" range of sodium ion densities, the counter pressure exerted by these ions does not significantly change the position of the magnetopause beyond the variability already documented in the literature. Therefore, over a wide range of exospheric parameters, the magnetopause exhibits a certain "stiffness" against the addition of sodium ions to the system. This is mainly caused by the shielding effect of recently ionized sodium particles which prevents the ambient electric field from immediately accelerating the bulk of the newly generated ion population near the dayside surface.