Field Line Resonance in the Hermean Magnetosphere: Structure and Implications for Plasma Distribution

The first statistical survey of field line resonance (FLR) events is presented using magnetometer data from the entire MErcury Surface, Space ENvironment, GEochemistry and Ranging mission. Ultralow‐frequency waves are an important tool for the magnetoseismology of the Hermean magnetosphere; this study provides a completely new window onto the resonance structures and plasma density distribution in the Hermean magnetosphere. Here we assess resonance events from two categories—toroidal resonances characteristic of the classical picture of FLRs in the terrestrial magnetosphere driven by the Kelvin‐Helmholtz instability and a more comprehensive approach including all observed transverse resonances with more relaxed polarization criteria. Two hundred twenty‐three toroidal FLRs with characteristics consistent with Kelvin‐Helmholtz‐driven FLRs are found in the dayside Hermean magnetosphere. The fundamental frequencies of these waves are used to provide estimates of plasma mass density in the range of ∼ 1–650 amu/cm3. A further 343 transverse resonances are found which provide very similar density estimates to the Earth‐like FLR population. Fundamental and harmonic frequencies from all 566 events are used to fit a power law to plasma mass density along the field lines. The equatorial plasma mass density is predicted to vary approximately with R−7.5. The offset of the Hermean dipole into the northern hemisphere causes significant asymmetries in the standing wave structure. Due to the extreme warping (away from a dipolar configuration) of Mercury's magnetosphere by the solar wind, the fundamental toroidal mode is predicted to oscillate with a notably lower frequency than the fundamental poloidal mode, contrary to relative toroidal and poloidal frequencies modeled for Earth's magnetosphere.


Introduction
Field line resonance (FLR) events have long been studied at Earth, where they provide a unique way of remotely probing the closed field line region of the terrestrial magnetosphere using both ground and spaced-based instruments. Toroidal FLRs oscillate azimuthally, approximately as shear Alfvén waves, where the frequency of the fundamental and its harmonics are defined by the plasma mass density along the field line and the field morphology. Magnetoseismology uses the known structure of the magnetic field and the frequency of observed FLRs to estimate the plasma mass density distribution of the resonant field line. In this study, we survey the MErcury Surface, Space ENvironment, GEochemistry and Ranging (MESSENGER) data set for FLRs which could be used to analyze the plasma distribution at Mercury.
The first suggestion that FLRs could exist in Mercury's magnetosphere was made by Russell (1989), when Mariner 10 magnetometer data revealed the presence of ∼2s period, linearly polarized oscillations. It was initially proposed that this wave was a fourth-harmonic FLR, corresponding to an ∼ 8s fundamental. This was later contested by Southwood (1997) who suggested that the wave could not have been a pure shear Alfvén wave as it had a significant compressional component-FLRs in the terrestrial magnetosphere are transverse waves, without a significant compressional component.
Previous studies of ultralow-frequency (ULF) wave activity within the Hermean magnetosphere have shown that waves often manifest with frequencies very similar to the local H + and He + gyrofrequencies, which are ∼ 1 Hz (e.g., Boardsen et al., 2012;Slavin et al., 2008). Because the frequencies observed are so close to the local ion gyrofrequencies, it has been suggested that the magnetohydrodynamic (MHD) approach to describing ULF waves at Earth used by Dungey (1954) In the terrestrial magnetosphere, ULF waves are often associated with driving by an energy source either internal to the magnetosphere (wave-particle interactions), or external to the magnetosphere (e.g., by solar wind buffeting or KH waves forming on the magnetopause). At Earth, wave-particle interactions include the drift and drift-bounce resonance mechanisms (Southwood et al., 1969), which allow the transfer of energy from gradient-curvature drifting energetic particle populations, such as those injected by substorms, to localized poloidal MHD waves (e.g., James et al., 2013;Yeoman et al., 2008Yeoman et al., , 2010. At Mercury, this instability is unlikely to develop because the magnetosphere may be too small to trap the drifting energetic particles. The bounce period is also comparable to the local proton cyclotron periods (Blomberg et al., 2007), so MHD treatment of wave-particle interactions is unlikely to be appropriate in the Hermean magnetosphere.
Interactions with the solar wind are likely to be a major source of ULF waves at Mercury. Buffeting of the solar wind has been suggested to produce a "ringing" of the Hermean magnetosphere (Baumjohann et al., 2006;Glassmeier et al., 2004), driving waves and thus energizing plasma. The KH instability has also been suggested as a likely source of ULF wave generation through FLR (James et al., 2016;Liljeblad et al., 2016;Liljeblad & Karlsson, 2017). The mechanism for KH waves driving FLRs at Earth is described in detail by Southwood (1974) using a box model magnetosphere and Chen and Hasegawa (1974) using a dipolar field geometry. In both cases, a KH wave on the Earth's magnetopause transmits pulsations planetward in the form of evanescent fast magnetosonic waves. The fast mode wave exhibits circular polarization with handedness and azimuthal wave number, m, which is dependent upon the local time: on the duskside of the magnetosphere, the azimuthal wave number m > 0 and the wave has a right-handed (RH) polarization; on the dawnside, where m < 0, the wave has a left-handed (LH) polarization. The azimuthal wave number is an inverse measure of the scale size of the wave, where a small m corresponds to a large-scale wave event, and the sign represents the direction of propagation, where m > 0 corresponds to eastward propagation and m < 0 is westward. At some location planetward of the magnetopause, the wave reaches a minimum in amplitude, where the polarization handedness reverses. Farther planetward of the minimum, if the 10.1029/2018JA025920 frequency of the evanescent wave matches the local toroidal eigenfrequency, the fast mode wave excites the shear Alfvén mode, forming a resonance. At the resonance, the wave peaks in amplitude, oscillating linearly in the azimuthal direction, where the handedness of the wave reverses once more. For the remainder of this paper, a resonant event which follows the description above will be referred to as an "Earth-like FLR." A number of studies have recorded observations of KH waves on the magnetopause at Mercury Liljeblad et al., 2014;Slavin et al., 2008;Sundberg et al., 2010Sundberg et al., , 2012, many of which were observed on the dusk side of the magnetosphere. This dawn-dusk asymmetry in KH waves has been reproduced using global kinetic hybrid simulations (Paral & Rankin, 2013) and is attributed to the finite gyroradius of ions broadening the shear layer at dawn and thus reducing the growth rate compared to dusk. The work by Sundberg et al. (2012) showed examples where ULF waves in Mercury's inner magnetosphere were driven by KH waves on the magnetopause, typically with periods ranging from 10 to 40 s. More recently, a study by Liljeblad and Karlsson (2017) on ULF waves driven by KH waves found that there were significant numbers of ULF waves with frequencies in the range of 20-40 mHz near dawn and stated that KH waves at dawn may be more important than previously thought. Further evidence in support of the solar wind driving of ULF waves is provided by James et al. (2016), where the transverse polarization (of ∼ 400,000 spectra) is consistent with the antisunward flow of the solar wind imparting opposite handedness polarizations upon the waves in both the dawn (∼ 230,000 spectra) and dusk (∼ 170,000 spectra) sides of the magnetosphere.
The sodium in Mercury's magnetosphere is the heaviest of the ions which make up a significant proportion of the plasma at Mercury and would therefore provide the lowest ion cyclotron frequency in the Hermean plasma. The sodium ion cyclotron frequency typically ranges from 40 to 140 mHz along the Y axis and 117-180 mHz along the X axis of the Mercury Solar Magnetospheric (MSM) coordinate system (where X originates at the center of Mercury's dipole and points toward the Sun). These values are estimated using the KT17 magnetic field model (Korth et al., 2015(Korth et al., , 2017; see section 2.1) from just within the magnetopause and at the planetary surface, where the magnetic field varies from 60 to 210 nT along Y and from 175 to 270 nT along X. In this paper, we extend the work from James et al. (2016) by searching for discrete Earth-like FLRs using MESSENGER MAG data. The resonant waves described here all exhibit frequencies far below those of the lowest local sodium cyclotron frequencies present at Mercury and are therefore unlikely to be related to the ion-Bernstein waves of Boardsen et al. (2012Boardsen et al. ( , 2015 or the IIH mode as suggested by Othmer et al. (1999) and Glassmeier et al. (2004). With frequencies well below the domains of IIH and ion-Bernstein waves, and wave power predominantly in the transverse directions (i.e., almost no compressional wave power), we assume that the resonances studied here are Alfvén waves that can be described using the MHD approximation. Using the resonance events found here, we construct an average plasma mass density model for the dayside magnetosphere of Mercury.

Resonant Wave Detection
MAG data from MESSENGER were collected from regions both inside and outside of the Hermean magnetosphere due to the elliptical orbit of the spacecraft. This study only utilized data which were collected from within the magnetosphere, so data from the magnetosheath and solar wind were discarded. Crossings of MESSENGER through the magnetopause could be detected using the MAG data, characterized by a sudden rotation in the magnetic field direction, accompanied by a change in the character of the fluctuations in the field (Winslow et al., 2013). Each of the remaining field vectors were then rotated from the MSM coordinate system into a coordinate system based on the KT17 model field (Korth et al., 2015(Korth et al., , 2017 at MESSENGER's location within the magnetosphere at the time corresponding to the time of measurement. The MSM coordinate system is defined such that X is directed from the center of Mercury's dipole to the Sun, Z is aligned with the magnetic dipole (and rotation axis), and Y points toward dusk. Each individual vector is rotated such that there is one component parallel to the model field (B || ), one toroidal component (B t ) that is positive eastward, and one poloidal component (B p ) that is in the direction of the local outward pointing normal of the field line. This is done in two steps: the first step rotates the X and Y axes around the Z axis, such that the Y direction is now pointing azimuthally around the planet (the toroidal direction), the second step rotates the Z and new X axes about the azimuthal (new Y ) axis, such that Z is aligned with the local model field direction.
Using a similar method to that of James et al. (2016), wave activity is to be detected by Fourier analyzing each component of the magnetic field data using a 120-s sliding window which allows the detection of frequencies 10.1029/2018JA025920 as low as ∼ 8 mHz, while minimizing spatial ambiguity. At each step, the sliding window progresses by 30 s, and the data within are linearly detrended. As the aim of this study is to investigate MHD wave activity at Mercury, only spectral peaks with frequencies lower than the local sodium ion gyrofrequency were used, to avoid contamination from cyclotron and IIH modes. In a given component of the field, each spectral peak has a corresponding Fourier power and frequency, which can be compared directly to the Fourier power of the other two components at the same frequency, allowing the determination of the predominant wave polarization. The Fourier powers and phases of the transverse (toroidal and poloidal) components of the magnetic field data can be used to determine the eccentricity, e, of the transverse polarization ellipse (Born & Wolf, 1980), where each frequency bin has its own polarization ellipse. The handedness of each of the transverse polarization ellipses is defined using the value of k i · B (Means, 1972), where k i is the wave propagation vector for fast Fourier transform frequency bin i, calculated using the complex result of the fast Fourier transform, and B is the local ambient magnetic field. The ambient field is taken to be the KT17 model field, and by the definition of the coordinate system used, B = (0, 0, B || ), so k i · B can be simplified to k i|| B || . When k i · B > 0, the wave at frequency i exhibits a RH polarization, while k i · B < 0 corresponds to a LH polarization.
In this study we use the KT17 model by Korth et al. (2015Korth et al. ( , 2017, which is based on earlier modular field models for Earth (e.g., Tsyganenko, 2013), where each field source (e.g., internal dipole or tail current) has its own module and a corresponding shielding field which contains the source field within the magnetopause boundary. At the midpoint of each time series used to create the Fourier spectra, the KT17 model is used to trace the position of MESSENGER to four sets of footprints where possible: 1. surface footprints in the northern and southern hemispheres of Mercury, 2. footprints on the surface of Mercury's iron core at a radius of 2,020 km (Hauck et al., 2013; the trace between these footprints will be used to model wave propagation), 3. footprints on a surface of "invariant latitude"-a Mercury-sized surface centered on the planets offset magnetic dipole, similar to the method used in Korth et al., 2014 (these footprints are convenient for creating plots because they put field line mappings to the northern and southern hemisphere on the same scale), 4. a footprint in the magnetic equatorial plane (this footprint is defined by the point along the trace where Z = 0 in the MSM coordinate system, coincidentally also where |B| along the trace is at a minimum, as the field model is symmetric about Z = 0).
In order to search for Earth-like FLRs-spectral peaks adjacent in time and close in frequency had to be stitched together to create a time series of wave power, frequency, and polarization as MESSENGER traversed the magnetosphere. Each time series was then scanned, three consecutive spectra at a time (a time series based on three time windows, each overlapping the next by 90 s), for characteristics that would suggest the existence of FLR: 1. the sum of the transverse components of power must be greater than the compressional (parallel) component, 2. there must be a peak in toroidal wave power, 3. a reversal in handedness of the wave must occur in the appropriate direction for a FLR.
The direction of the flip in polarization depends on both the local time and the direction of travel of MES-SENGER: when MESSENGER is traveling planetward, the reversal would be RH to LH in the dawnside magnetosphere and LH to RH on the duskside, where the reverse is true for when MESSENGER is traveling away from Mercury.

Estimating Plasma Mass Density
The Singer et al. (1981) wave equation, can be used to accurately model wave propagation in an arbitrary magnetic field configuration. In this equation, s is the distance along the field line, is the displacement of the field line, V A is the Alfvén velocity, B 0 is the ambient magnetic field magnitude, and is the angular frequency of the wave. The parameter h is related to the separation of two adjacent field lines in the direction of the field displacement by where is the normal separation of both field lines at some point on the field line, and d(s) describes the normal separation of the two field lines at any other point along the field line. The wave equation can be solved numerically using a Runge-Kutta algorithm (e.g., Waters et al., 1995) given a model magnetic field and a plasma model. The frequencies of the wave harmonics can be found using the shooting method, where the frequency is varied until the boundary conditions of the wave are met.
To model the field displacement due to the ULF wave activity in the Hermean magnetosphere, a trace of the KT17 model provided the morphology of the field line, with the magnetic field strength, B 0 , along the field line. The assumption has been made that the wave is reflected at the surface of Mercury's iron core in both hemispheres, so the boundary conditions for the resonances require the existence of a wave node (i.e., no field displacement) at the footprints of the field lines on the core of the planet. At Earth, the boundary conditions for the waves would be such that they reflected at the ionosphere at each end of the field line, where the conductivity is high; at Mercury there is no ionosphere, and the height-integrated Pedersen conductivity above the surface is very low (Lammer & Bauer, 1997), and instead the iron core may be used to anchor the field line and close the currents associated with FLRs. For toroidal waves, h was calculated using two field lines which mapped to the same radial distance in the magnetic equatorial plane, with a separation in magnetic local time (MLT) equivalent to 0.1 R M at the equatorial footprint. To model poloidal waves, the equatorial footprints of the waves had the same MLT but were separated by 0.1 R M in the radial direction. For simplicity, the plasma mass density along the field line was modeled using a power law, where defines the shape of the distribution along the field line, eq is the equatorial plasma mass density, R is the radial distance from the center of the Hermean dipole, and R max is the farthest point along the field line. The region of space between the crust and the core of Mercury was assumed to be a continuation of the magnetosphere, using equation (3) to describe the plasma mass density. Another model was created where this region of space was assumed to be transparent and that the wave propagated at the speed of light with a refractive index based on that of SiO 2 , but this made almost no difference to the predicted plasma mass densities. This is because there is relatively little field displacement near the magnetic footprints (i.e., the nodes of the wave) compared to that near the magnetic equator.
To determine the plasma mass density along a field line, given a fundamental wave frequency, f , and a fixed power law, , the downhill simplex method (Nelder & Mead, 1965) was used to vary eq in order to minimize the difference between the fundamental solution to equation (1) and f . This was performed for each event which qualified as a resonance in section 2.1, using integer steps of in the range 0-6. Figure 1 shows an example of a possible Earth-like FLR detected using the method described in section 2, shortly after MESSENGER had crossed through the dawnside magnetopause at ∼ 10:28 UT on 27 May 2014. Figure 1a shows the poloidal (red), toroidal (green), and compressional (blue) components of the magnetic field, where wave activity is exhibited predominantly in the toroidal component from ∼ 10:29 to 10:34 UT. Figure 1b shows the Fourier power for each component, using the same color scheme as in Figure 1a, with the total power in black. Figure 1c shows the frequency of the wave as MESSENGER moves planetward through the region of wave activity. Figure 1d shows the wave polarization ellipses calculated with the Fourier coefficients from each pair of transverse spectra using the method described in section 2.1; there is a single ellipse plotted for each time window, where the frequency bin used each time corresponds to that of the peak transverse power (toroidal + poloidal) in that window. The toroidal amplitude is in the vertical axis, poloidal amplitude is along the horizontal axis, RH polarized ellipses are in red, and LH polarized ellipses are green. Figure 1e shows the radial distance from the center of Mercury's dipole to the magnetic equatorial footprint of MESSENGER in orange and the local time of the equatorial footprint in blue. The pink shaded area in panels a-e depicts the 120 s time window which is identified here as containing the wave resonance, where Figure 1f shows the power spectra for each of the three components of the magnetic field during this time window (the polarization ellipse which corresponds to the transverse components of these spectra lies at the center of this region in Figure 1e, intersected by the vertical dashed line). The gray shaded area represents the time window which corresponds to a reversal in polarization near the minimum in wave amplitude. At both reversals, the polarization becomes linear as expected. Figure 1g demonstrates the orientation of the coordinate system used for Figures 1a and 1b.

Results
The method described in section 2.2 was used to estimate the plasma mass density along the resonant field line for the example wave in Figure 1 at the fundamental frequency of ∼ 25 mHz. Figure 2a shows the solutions of the equatorial plasma mass density, eq , for power laws in the range 0 ≤ ≤ 6, where eq is between 95 and 115 amu/cm 3 . Figure 2b demonstrates how the plasma mass density would be distributed along the length of the field line, x, for each integer step of the power law. The southern portion of the field line (starting at x = 0) is on the left of the figure, the northern portion of the field line is on the right of the figure, and the magnetic equator is represented by a vertical dashed line.  Figure 3 shows the field line displacement, , along the length of the field line, x, for the first three poloidal (panel a) and toroidal (panel b) harmonics of the resonant field line from Figure 1. The second harmonic of the toroidal wave has a frequency of ∼ 58 mHz, which is close to the observed second peak in toroidal wave power in Figure 1f. The functional form of h in equation (1) differs greatly between the toroidal and poloidal components of the wave, resulting in a ∼ 10 mHz difference in the two fundamental frequencies and a considerable difference between the morphologies of the perturbed field lines. The displacement of the field line in Figure 3 is asymmetric about the magnetic equator because of the offset of the Hermean dipole relative to the planet, which shifts the largest displacement of the field line northward of the magnetic equator for the odd modes.
The method of detecting Earth-like FLRs outlined in section 2.1 was able to find 223 candidate resonance events. Figures 4a and 4b show the magnetic mapping of these events to the equatorial plane (a) and on the surface of invariant latitude (b), where both panels are oriented such that noon is at the top and dawn is to the right. Most events (212 of 223) are found in the dayside magnetosphere, particularly between dawn and noon, mapping to midlatitudes (∼ 30-60 • ) on the surface of invariant latitude.
The modal frequencies for the events in each spatial bin of the magnetic equatorial plane and the invariant-latitude surface are shown in Figures 4c and 4d, respectively. Most of the events appear to exhibit resonant frequencies between ∼ 15 and 50 mHz, and while it is likely that there is some relationship between the frequency of the resonance and the location of the resonant field line, it is difficult to draw any conclusions from these data, as there are so few events in each spatial bin. Figure 4e shows the mean predicted equatorial plasma mass density (assuming = 3) mapped to the magnetic equatorial plane. Even though there are very few events, there is a definite trend in the predicted equatorial plasma mass density, where the highest values (close to 500 amu/cm 3 ) map to locations in the equatorial plane closest to the planet, and the lowest values of plasma mass density (∼ 1 amu/cm 3 ) map to  the equatorial plane much farther from the planet. Figure 4f shows mapping of the densities in Figure 4e to the invariant latitude surface, where the highest densities are predicted on low-latitude field lines.
The spatial coverage of the small number of classical toroidal FLR events observed by MESSENGER is sparse, particularly near the flanks. At Earth FLRs can often be driven by cavity and waveguide modes (Allan & Poulter, 1992), where there is the added possibility of extra turning points, evanescent regions, and resonances. Also, if Earth occupied the same proportion of its magnetosphere as Mercury does, much of the resonant regions that are currently observed in the Earth's magnetosphere would be below the planet's surface. This leaves a relatively much smaller distance between the magnetopause and the planet to form a classical FLR, so it could be that in some circumstances, the fast mode wave might behave differently to the classical picture at Earth.
For the reasons stated above, FLRs could exist with polarization reversals which do not match the classical "Earth-like" picture, where there are exactly two polarization reversals of the incoming wave. Due to the warped nature of the Hermean field (distorted away from a dipolar configuration), a truly toroidally oscillating field line may be unrealistic-there could also be a significant poloidal component. To investigate this, the data analysis above has been modified and relaxed below to include anything that appears to be a transverse resonance. The required criteria for this new event selection condition are the following: 1. Transverse wave power must be dominant of the compressional wave power. 2. There must be a peak in either toroidal or poloidal wave power. 3. There must be a reversal in polarization across the peak in amplitude, but the direction in which the handedness switches is not constrained.
This relaxed set of criteria produced 566 events in total, a small number of which were poloidal (47 of 566), but most were toroidal. In order to compare with the Earth-like events, Figure 5 shows the predicted equatorial plasma mass density assuming a power law index of = 3 against the radial distance of the magnetic equatorial footprint, where classical Earth-like events are red and the more relaxed transverse events (including the Earth-like FLRs) are gray. The orange line presented in Figure 5 shows a least squares fit to the densities from all events. Figure 6 shows where the full set of 566 events mapped to in the magnetosphere (a and b), their frequencies (c and d), and the predicted equatorial plasma mass density assuming a power law with = 3 (e and f) in the same format as Figure 4. The coverage of these events is much more complete over the dayside magnetosphere. The modal frequency appears to be slightly higher closer to the surface of the planet, but there are still relatively few events in most of the spatial bins. The estimated eq depicts a very similar trend to that in Figure 4, albeit with a more complete spatial coverage.
The trend line in Figure 5 showing the overall fit to the radial dropoff of eq (calculated assuming = 3) with the radial distance of the equatorial footprint is linear in log space, which corresponds to a power law in linear space of the form where the plasma mass density appears to drop off with ∼R −7.5 . Table 1 contains the values a and k with uncertainties calculated for all values of in the range 0-6, all of which have a very similar profile.

Example Event
The example event presented in section 3 and Figure 1 exhibits some very similar features to what one might expect to observe in an Earth-like FLR as described by Chen and Hasegawa (1974) and Southwood (1974). Shortly after the crossing through the magnetopause, the eccentricity of the transverse wave polarization is at its most circular, where the handedness (LH) is consistent with driving by the antisunward flow of solar wind on the dawn magnetopause. At ∼ 10:30 UT while MESSENGER traverses planetward, the polarization reverses, while the amplitude of the wave is reducing-possibly indicating that the spacecraft had traversed through the turning point region, where the wave becomes evanescent. In the evanescent region, the wave is highly elliptical, as predicted by Chen and Hasegawa (1974). Figure 1b shows that the toroidal wave power is dominant through the whole time series and that there is a large peak in the toroidal wave power at ∼ 10:32-10:33 UT; the transverse polarization becomes linear and reverses in handedness across a region where the toroidal wave power reaches a maximum intensity-indicating that MESSENGER had passed through a resonant region. At the resonant region, Figure 1f shows that the wave exhibited a frequency of ∼ Note. The "⋇" symbol represents the multiplicative equivalent to "±". 25 mHz, with a small secondary peak in wave power at ∼ 58 mHz. Soon after, the wave amplitude detected by the spacecraft dropped off quickly as it moved farther planetward, away from the resonant region.
Using MESSENGER's position as it crossed the resonant region, a Runge-Kutta routine was used to trace the spacecraft position along the KT17 model magnetic field in order to determine the morphology of the local field line. The field trace and the fundamental frequency of the observed toroidal wave (25 mHz) was then used to solve equation (1) using the method described in section 2.2 to determine the equatorial plasma mass density, eq , for power laws in integer steps from = 0 to = 6. Figure 2a shows that eq takes values within the range 95-115 amu/cm 3 , where the highest values correspond to the lowest power index, = 0. James et al. (2016) estimated a much higher density of ∼ 240 amu/cm 3 for the same wave event, but the method used to make that estimate assumed that the plasma mass density was constant along the field line and did not include the portion of the field line which existed between the surface and the core of Mercury (i.e., the wave was reflected at the planet's surface). James et al. (2016) calculated this value by inverting a simple time-of-flight calculation, a method which is invalid when the wavelength is comparable to the size of the system as it provides a poor estimate of the fundamental wave frequency (Singer et al., 1981). Figure 2b shows the plasma mass density profile along the field line calculated for all seven power indices; there is relatively little difference in plasma mass density near the equator, but the difference at the footprints of the field lines are significant. Figure 2 indicates that, while the power law index, , vastly alters the inferred plasma mass density at the extrema of the field line, it makes relatively little difference to the inferred plasma mass density near the equator. The relatively small difference in plasma profile near the equator arises because the outer portion of the field line where the Alfvén speed, V A , is smallest is the most significant region for defining the wave period, (Denton & Gallagher, 2000). For the fundamental wave period, the footprints are less important for the propagation of the wave, as they are expected to contain the wave nodes, as opposed to the wave antinode in the outer portion of the field line, where the displacement of the wave electric field is slowed down by the mass of the plasma. Figure 3 shows the wave displacement along the field line, obtained by solving equation (1) assuming a power law index = 3 for the first three poloidal (panel a) and toroidal (panel b) harmonics. In both panels, it is apparent that approximately 60% of the field line exists north of the magnetic equator, which is due to the magnetic dipole being offset toward the northern hemisphere of the planet. This offset results in the northern magnetic footprint mapping to a much higher latitude than that of the southern footprint. The offset also produces a further northward shifting of the fundamental antinode; at Earth the waves antinode might be expected to exist near the equatorial plane (e.g., Denton et al., 2004;Takahashi et al., 2004), but our solution for the Hermean magnetosphere shows that the antinode is displaced significantly into the northern portion of the field line, particularly for the toroidal mode.
Due to the warped geometry of the magnetic field lines within the magnetospheres of Mercury and Earth, the solutions of equation (1) for poloidal and toroidal modes often produce different frequencies, particularly for the fundamental mode. Cummings et al. (1969) showed, for a range of power law indices in the range = 0 to = 6, that the poloidal frequency is expected to be lower than the toroidal frequency of the same field line in a dipolar field configuration. Figure 3 shows that, in this particular region of Mercury's magnetosphere, the poloidal mode has a significantly higher frequency than that of the toroidal mode. The reason for this may be the more extreme warping of Mercury's magnetic field, where the outer parts of the field lines are swept tailward of their planetary footprints (according to the KT17 model), compared to the Earth's relatively more dipolar magnetic field. This warping appears to make the antinode of the toroidal wave much wider than that of the poloidal wave, possibly indicating that the region of lowest Alfvén speed is spread further along the field line for the toroidal mode when compared to that of the poloidal mode. The second and third harmonics of both modes have much more similar frequencies, as predicted by Cummings et al. (1969).

Standard Earth-Like FLR Events
In total there were 223 Earth-like resonance events found. Figure 4 shows that almost all of the events are found on the dayside magnetosphere, mapping to midinvariant latitudes. This is the same region of the 10.1029/2018JA025920 Hermean magnetosphere which was found to contain a large amount of toroidally polarized wave activity by James et al. (2016). Approximately 60% of these events map to the dawnside of the magnetosphere and only 40% on the duskside, contrary to the results obtained by Liljeblad et al. (2016), where more waves were observed toward dusk. The frequencies of most of these events were within the range 15-50 mHz-directly comparable to those wave studied by Liljeblad et al. (2016) and Liljeblad and Karlsson (2017), which were associated with KH as their driving energy source, and the toroidal waves in Figure 8 of James et al. (2016). The modal frequencies of the spatial bins in Figure 4 do not show much of an obvious pattern, but this is likely due to the low numbers of events found in each bin. The larger number of dawnside events supports the work of Liljeblad and Karlsson (2017), where they concluded that the dawnside KH may be more important than previous research would suggest.
The equatorial plasma mass density shown in Figure 4 appears to increase monotonically toward the planet, where plasma mass densities are predicted to be up to a few 100 amu/cm 3 . Locations farther out, which map to magnetic latitudes greater than ∼50 • , have average mass densities of less than 100 amu/cm 3 . The density estimates produced here are higher than the observations made by MESSENGER's Fast Imaging Plasma Spectrometer (FIPS) instrument. Studies using FIPS Raines et al., 2011Raines et al., , 2013Raines et al., , 2014 have found the average density of H + in the region of 1-20 cm −3 , Na + around 5.1 × 10 −3 cm −3 (or 2 cm −3 in the cusps), He 2+ at 3.9 × 10 −2 cm −3 , He + at 3.4 × 10 −4 , and O + at 8.0 × 10 −4 cm −3 , which would add up to a mass density orders of magnitude lower then those predicted closest to the planet. However, the measurements made by FIPS could be in agreement with our predicted plasma mass densities at larger distances from Mercury. The FIPS instrument detects ions within the energy range of 50 eV to 20 keV (Andrews et al., 2007) and may be unable to detect the presence of cold ion populations with energies which are beyond the scope of this instrument. Also, approximately 30% of the FIPS field of view was obstructed by MESSEN-GER's sunshade (Raines et al., 2011) which could contribute toward the difference in density measurements, especially when the center of the plasma velocity distribution is not within the instrument's field of view.
Modeling of the Hermean magnetosphere and sodium ion production have predicted densities closer to those predicted in this study. Benna et al. (2010) modeled the Hermean magnetosphere during the first MESSENGER flyby of Mercury, where proton densities ranged from 8 to 10 cm −3 within a small drift belt, > 10 cm −3 around the morning sector, and 10-100 cm −3 in the cusps. It is unlikely that after being injected into the nightside of the magnetosphere, there would be enough plasma capable of drifting around to the dayside to account for the plasma mass densities predicted by this study, but the source may be through the photoionization of exospheric sodium on the dayside. Attempts to model the sodium ions in the magnetosphere have predicted number densities in the range of ∼ 10-100 cm −3 Leblanc et al., 2003;Yagi et al., 2010), corresponding to 230-2,300 amu/cm 3 , where the highest concentrations of sodium are expected in the dayside magnetosphere. The main sources of the sodium exosphere are through photon-stimulated desorption, thermal desorption, and solar wind sputtering, with some contribution from micrometeoroid vaporization and magnetospheric ion sputtering (Leblanc & Johnson, 2010); the atoms near the dayside of the magnetosphere are then photoionized to form part of the Hermean plasma. Another problem is that Delcourt et al. (2003) showed that sodium ion with 1 keV energy do not behave adiabatically in most of the Hermean magnetosphere-only within a thin shell directly above the planet's surface. This would lead to many sodium ions being pitch angle scattered into the loss cone; but for colder ions, if they exist, this region could extend to greater distances. When BepiColombo reaches Mercury in late 2025, the Plasma Wave Investigation (Kasaba et al., 2010) instrument onboard the Mercury Magnetospheric Orbiter will be able to sample the upper hybrid or plasma frequency to elucidate the total plasma number density.

Other Transverse Resonances
The resonance events shown in Figure 5 which do not follow the classical picture of an Earth-like resonance all provide very similar inferred plasma mass densities when compared to the original 223 Earth-like events. Both sets of predicted plasma mass densities follow the same trend, with an overall monotonic decrease in plasma mass density with distance from Mercury. The densities predicted here are also in good agreement with those predicted by James et al. (2016) based on the average toroidal frequencies observed at different L shells. The majority of this population of waves exhibit predominantly toroidal polarizations, but a small number of the new events are poloidally polarized. At Earth, poloidal oscillations are considered to be unreliable for the purposes of estimating the plasma mass density as their frequency can be reduced by the presence of radial pressure gradients (Denton et al., 2003;Vetoulis & Chen, 1996). If this is the case at Mercury, it should make little difference when combined with the toroidal events because there are so few poloidal events. Figures 6a and 6b show that the overall distribution of the resonance events is very similar to those in Figures 4a and 4b, where there is an abundance of resonance events in the dayside magnetosphere. In total there are an extra 343 resonances that do not fit the Earth-like description, where ∼ 60% of the observed events occurred toward the dawn side of the magnetosphere, much like the Earth-like events described in section 4.2. The similarities are also obvious when comparing the frequencies presented in Figures 6c and 6d, which show that the frequencies that are typically observed are the same as the Earth-like FLRs. Given that the wave characteristics of the new population of resonance events are so similar to the Earth-like events, it is logical that the inferred densities that are presented in Figures 6e and 6f are also very similar to those in Figures 4e and 4f. In both cases, the highest densities are predicted to exist on field lines very close to the planet, which map to low invariant latitudes.
The fits that are produced in Table 1 are global fits to all events from all local times at all times during the Hermean year. Similar fits could be found for individual MLT sectors, but there would be relatively few events for the fit. The uncertainty in the global fit is likely to be large as there is a lot of spread in the density estimates for a small range in R eq . The spread in density estimates probably arises due to several factors including changes in mass loading of the magnetic field and changes in field morphology. Changes in the mass loading of the magnetosphere can arise when magnetic activity at Mercury changes; extreme solar wind conditions and large magnetic fluctuations have been shown to increase sputtering rates and therefore also plasma mass densities (Poh et al., 2016;Raines et al., 2014). Extreme coronal mass ejection events such as those of Slavin et al. (2014) can also drastically affect the magnetospheric size and shape. The elliptical orbit of Mercury also means that the magnetosphere experiences a wide range of solar activity levels, where the interplanetary magnetic field strength doubles (James et al., 2017) and solar wind dynamic pressure goes from ∼ 11.0 to ∼ 26.5 nPa (Fujimoto et al., 2007) from aphelion to perihelion. These differences in interplanetary magnetic field and dynamic pressure would both change the magnetospheric activity levels and change the size of the magnetosphere, thus changing the natural eigenfrequencies of the closed field lines.  (5) for a wide range of power law indices in the range −12 ≤ ≤ 12 on the vertical axis and varying ranges of equatorial plasma mass density on the horizontal axis. Each panel corresponds to equatorial footprints within different ranges of R eq labeled in the top left corner of each panel. The contours emphasize the gradients in the frequency misfit, and a green circle shows the optimal values for and eq .
In the terrestrial magnetosphere, the typically accepted rates of radial plasma number density dropoff in the equatorial plane at Earth are usually in the range R −3 -R −4 (e.g., Poulter et al., 1984;Warner & Orr, 1979). A similar number density profile which varies with R −4 could exist at Mercury if the average ion mass varied with ∼ R −3.5 . This would correspond to an abundance of heavy ions (e.g., sodium) at low altitudes and low mass ions (protons) forming most of the plasma at larger radial distances.
Using the average density dropoff with radial distance provided in Figure 5, and assuming a density profile along the field lines ( ), an average model for plasma mass density can be created for Mercury's dayside magnetosphere. Figure 7 shows this average density model applied to the dawn (a), noon (b), and dusk (c) sectors of the magnetosphere, using = 3. The green dashed lines in each panel represent the altitude range within which the MESSENGER MAG data was sampled over the duration of the mission. This model suggests that the majority of the magnetosphere near the flanks consists of low (< 10 cm −3 ) densities, while high densities (> 100 cm −3 ) are prevalent close to the surface of Mercury, particularly in the subsolar region.

Harmonics and Plasma Mass Density Profiles
Previous studies of toroidal ULF waves in the terrestrial magnetosphere have made use of power laws in order to describe the mass density profile along the field lines (e.g., Cummings et al., 1969;Denton & Gallagher, 2000;Schulz & McNab, 1996;Waters et al., 1995). The power law index has been shown to alter the relative frequencies of the wave fundamental and its harmonics, so the ratio of an observed waves fundamental frequency and its harmonics could be used to infer the variation along the field line. Wharton et al., 2018).
Many of the spectra of the 566 transverse resonances from section 4.3 exhibited multiple spectral peaks representing possible harmonic frequencies, such as that observed in Figure 1f. Figure 8 shows the use of these harmonics to determine the plasma mass density profile along the field lines for eight equally spaced regions between R eq = 1.1 and R eq = 1.9. Each panel has the equatorial plasma mass density (log( eq )) along the horizontal axis and the power law along the vertical axis. The contours and color profile show the frequency misfit (log(Δ RMS )) between the modeled harmonic frequency f mi obtained using the Runge-Kutta In each of the panels of Figure 8, a green circle marks the minimum in the frequency misfit for each region, but as the contours show, the optimum in each plot is very shallow and is part of a much larger band of low f RMS which spreads across a large range in . The optimal values of and eq for each of the panels in Figure 8 are given in Table 2. The optima in the regions closest to Mercury in the range 1.1 ≤ R eq ≤ 1.7 suggest that density profiles which vary with ≈ 6 − 10 would fit with the observed waves and their harmonics-indicating that there is significantly more plasma near the field line footprints than there is near the equator, but these minima are very poorly constrained. For R eq ≥ 1.7, there is a sudden change, where becomes very small, ≈ − 1 − 0, implying that the plasma mass density may have a peak near the magnetic equator at large distances from the planet.
The lack of constraint in the optima for each R eq range in Figure 8 suggests that the values in Table 2 should be used cautiously. The reason for such a shallow minimum in f RMS is found in the large spread of the frequencies in each radial bin. The standard deviations vary between 3 and 9 mHz, between 10 and 21 mHz, and between 10 and 24 mHz, for the first, second, and third harmonics, respectively. This analysis also combined the harmonics observed from all MLT sectors in the dayside magnetosphere, which is likely to have added to the uncertainty in the predicted power law parameters; the field morphology near noon is very different to that near the flanks of the magnetosphere, so harmonic ratios would be different for the same L-shell. It is also possible that the distribution of plasma mass density along those field lines might vary with MLT. Another factor that could be causing the variations in harmonic frequencies is in the fluctuations in field morphology and mass loading due to the changes in activity experienced by Mercury's magnetosphere from aphelion to perihelion as discussed in section 4.3. The assumption that the plasma mass density varies monotonically from the field line footprints to the equator could also be incorrect; previous studies at Earth have used nonmonotonic functions to describe the plasma mass density profile, where a common feature is a bulge of plasma density near the equatorial plane (e.g., Denton et al., 2004Denton et al., , 2006Sandhu et al., 2017). A more reliable model for the plasma mass distribution may be possible if there were enough transverse waves to be split into MLT bins and magnetospheric activity bins.

Conclusions
Using the MAG data collected by the MESSENGER mission, 223 events which resemble the classical picture of terrestrial FLRs are found and analyzed. These waves are used to provide an estimate for the plasma mass density, given an assumption of the plasma mass density profile along the field lines.
The Earth-like FLRs are compared to a further 343 transverse resonances which do not conform to the typical depiction of a terrestrial FLR. Overall, the density predictions obtained from both populations of waves were found to be very similar, where plasma mass density was highest at around 500 amu/cm 3 near the planetary surface, and lowest for the highest L-shells with values around 1 amu/cm 3 .
A power law fit to the equatorial plasma mass density with radial distance reveals that it varies with approximately R −7.5 . If an R −4 variation in number density is assumed, then a further variation of average ion mass with ∼ R −3.5 would be required, which could be realized if the plasma closest to Mercury's surface consists mostly of sodium ions, and plasma farther from the planetary surface is made up of mostly protons.
More work is required to provide a reliable plasma mass density model for the Hermean magnetosphere. If the field lines in the magnetosphere "ring" at their natural frequencies as it has been suggested by Glassmeier et al. (2004), then it may be possible to use these small amplitude oscillations and their harmonics to improve on the predictions made here. Further work using FIPS could also provide more measurements to help elucidate the plasma mass density profile along the field lines. The upcoming BepiColombo mission