An Empirical Model for Mercury’s Field‐Aligned Currents Derived From MESSENGER Magnetometer Data

Mercury is the only planetary magnetosphere which does not possess a significantly conducting ionosphere, yet magnetic field observations by the Mariner 10 and MErcury Surface, Space ENvironment, GEochemistry, and Ranging (MESSENGER) missions revealed a stable presence of large‐scale field‐aligned currents (FACs). Several empirical magnetic field models have been developed to describe Mercury’s magnetospheric currents and their associated magnetic fields, but none attempted to include an explicit description of large‐scale FACs. Here, we describe a dynamic FAC magnetic field model for Mercury based on the analytical solution to conical currents, which are shielded by a dynamic magnetopause boundary. The model contains free parameters setting the FACs' latitudinal extent and amplitude as functions of magnetospheric activity. The parameters are fit by minimizing the root‐mean‐square (RMS) differences between the model and MESSENGER Magnetometer data. During magnetically quiet conditions, the modeled FACs have an intensity of ∼10 nA/m2 and extend in latitude from ∼83°N to ∼40°N. They intensify to ∼20 nA/m2 and expand equatorward to ∼28°N during the most active times. The inclusion of the FAC model reduces low‐altitude RMS residuals by 8.1% when compared to prior models. The model effectively captures the azimuthal component of the magnetic field present in the MESSENGER low‐altitude data, but largely misses the radial and co‐latitudinal components, indicating other large‐scale physics remain missing in the mathematical descriptions of Mercury's magnetosphere.


Introduction
The Mariner 10 spacecraft performed three flybys of Mercury during the years 1974 and 1975 and discovered a relatively weak dipolar planetary magnetic field and concomitant magnetosphere (Ness et al., 1974(Ness et al., , 1975)).The observations acquired by this spacecraft demonstrated that Mercury's magnetosphere is qualitatively similar to that of other magnetized planets in that it possesses a magnetopause boundary and a magnetotail, where the northern and southern lobes are separated by a cross-tail current sheet.The magnetic field observed during the Mercury flybys further revealed signatures of field-aligned currents (FACs), which were interpreted to correlate with substorm like particle injections (Slavin et al., 1997).FAC systems play a fundamental role in mediating stresses within planetary magnetospheres (Cowley, 2000).Of all the bodies in our solar system known to have magnetospheres, Mercury is unique in that it does not possess a significantly conducting ionosphere.Given that this region provides for closure of FACs at Earth, the presence of large-scale FACs in absence of an ionosphere was somewhat unexpected (Glassmeier, 1997(Glassmeier, , 2000)).The morphology and dynamical evolution of magnetospheric current systems, including FACs, can be quantified and analyzed using empirical magnetic field models (Tsyganenko, 2013).These models enable the extraction of system-level information from data by fitting analytical representations of planetary magnetic fields to collections of magnetometer observations.A spacecraft's magnetometer inherently observes the total magnetic field, B tot = B int + B ext , which is composed of contributions from internal, B int , and external sources, B ext .B int and B ext are qualitatively different from one another in that the former is created in the dynamo region in the planet's interior and from crustal magnetization (Hood et al., 2018;Johnson et al., 2015Johnson et al., , 2018) ) whereas the latter is generated by electric currents flowing in space.B int changes over relatively long timescales (years to decades) and co-rotates with the planet.Meanwhile, ion kinetic processes drive variations in B ext on the order of 1 Hz at Mercury (Boardsen et al., 2012(Boardsen et al., , 2015) ) and global convection timescales, supported by substorm observations, are on the order of 2 min (Imber & Slavin, 2017;Siscoe et al., 1975;Slavin et al., 2010;Sun et al., 2015).Further, B ext displays 88-day periodic signatures (Wardinski et al., 2019) driven by variations in the solar wind dynamic pressure caused by Mercury's eccentric orbit (Baker et al., 2013;Dewey et al., 2015).The geometry and amplitude of the magnetospheric current systems generating B ext are thereby largely dominated by the solar wind speed and flow direction.The internal field can be described by a gradient of a scalar potential function, B int = ∇ψ int , which is, however, not possible for the external field because the magnetic field in the region of space traversed by the spacecraft is not curl free due to the presence of magnetospheric current systems.For this reason, the descriptions of the internal and external fields must be separated in order to understand both the dynamo mechanisms in the planet's interior and the processes that dominate in the magnetosphere.Several attempts were made to characterize B int based on the Mariner 10 flybys, but due to the lack of spatial coverage, the analysis is not unique and resulting models were inconsistent with one another (Connerney & Ness, 1988;Genova et al., 2021).
A more detailed understanding of Mercury's magnetosphere would await the arrival of the MESSENGER spacecraft, which entered orbit around the planet on 18 March 2011.MESSENGER orbited Mercury in a highly eccentric (periapsis of 200-600 km and an apoapsis of 10,000-15,000 km) and high inclination (82°) orbit until running out of propellant and impacting the surface on 30 April 2015, after completing over 4,000 orbits.Among others, the MESSENGER spacecraft was equipped with a tri-axial fluxgate magnetometer (Anderson et al., 2007).Even prior to Mercury orbit insertion, the Magnetometer on MESSENGER was revealing new insights into the structure of Mercury's magnetic field.Fitting a paraboloidal empirical magnetic field model to a combined magnetometer dataset consisting of Mariner 10 flyby data and data from the first two MESSENGER flybys, which occurred in 2008, Alexeev et al. (2010) determined that B int may be well described as an offset magnetic dipole.This study determined Mercury's dipole to have a magnetic moment of 196 nT R M

3
, where R M is the mean radius of Mercury, 2,440 km, with an offset 405 km northward from the planet's geometric center along the rotation axis.It also found the dipole axis was shifted ∼4°from the rotation axis.The orbital data would later confirm and refine this finding.An examination of MESSENGER crossings of the magnetic equator constrained the offset value to 479 km (0.196 R M ) and indicated that the tilt angle between the magnetic dipole axis and the spin axis is < 0.8° ( Anderson et al., 2012).Meanwhile, a parabolic model derived from orbital data confirmed an axially aligned offset dipole with a magnetic moment of 190 nT R M 3 (Johnson et al., 2012).More recently, attempts to characterize B int include localized models of the northern hemisphere (Oliveira et al., 2015;Thébault et al., 2018) and global spherical harmonic models (Wardinski et al., 2019(Wardinski et al., , 2021)).These studies found markedly lower dipole offsets ranging from 0.106 R M (Oliveira et al., 2015) to 0.155 R M (Wardinski et al., 2021).The upcoming BepiColombo mission will enable a more precise determination of the internal field (Toepfer et al., 2022).
Because the geometry of the magnetosphere is defined with respect to the magnetic dipole center rather than the planet's geometrical center, it is convenient to adopt a coordinate system with its origin at the dipole center.A widely used coordinate system in magnetospheric studies is the Mercury Solar Magnetospheric (MSM) coordinate system, in which the +x axis points from the dipole center to the center of the sun and its +z axis aligned along the spin/dipole axis, and the +y axis completes the right-handed system.Another commonly used coordinate system is the Mercury Solar Orbital (MSO) coordinate system, in which the axes are parallel to the MSM system and centered at the planet's geometrical center.Due to the nearly zero obliquity of Mercury's rotation axis, the +z axis of the MSO and MSM coordinate systems also coincides with the normal of the planet's orbital plane.
The MESSENGER Magnetometer observations also confirmed the presence of FACs at Mercury (Anderson et al., 2014) and revealed several new details about their large-scale distribution.By examining segments of the MESSENGER orbit when the spacecraft traverse the northern pole at low-altitudes (altitude <1,000 km), it was found that the FACs show a stable presence at Mercury and are not solely transient events associated with individual substorms as had been suggested earlier (Glassmeier, 1996).Indeed, the FAC signatures were quite steady during the ∼10 min traversals of the low-altitude northern hemisphere by the MESSENGER spacecraft compared with the typical ∼3 min duration of a substorm at Mercury (Imber & Slavin, 2017).The magnetic perturbations associated with the FACs do not rotate with the planet, are consistently oriented sunward at latitudes northward of 60°, and are qualitatively similar to Region-1 FACs at Earth (Anderson et al., 2018).However, no evidence of Region-2 type FACs was observed at Mercury.In absence of a conducting ionosphere, these authors concluded that the FACs must close within the planet, adopting a closure scenario proposed by Janhunen and Kallio (2004).In this scenario, Mercury's FACs penetrate vertically through the planet's resistive mantle before closing laterally along a higher conducting layer near the core-mantle boundary.The geometry of current closure was simulated using magnetohydrodynamic (MHD) (Dong et al., 2019) and hybrid MHD-kinetic simulations (Exner et al., 2020), although the latter study demonstrated that a dense enough sodium exosphere allows for a partial closure of the FACs through Pederson and Hall currents.Further, although the FAC signatures were steady on ∼10 min timescales, their intensity changed notably from orbit to orbit, that is, on an 8 12 hr or faster cadence.
In order to understand the statistical FAC patterns, Anderson et al. (2014) spatially averaged the MESSENGER magnetometer vectors on polar grids and binned them by the magnetic disturbance level (Anderson et al., 2013), which allowed the computation of the spatial distribution of the currents using Ampère's law.The analysis confirmed a Region-1-type FAC distribution with current densities of 50-200 nA/m 2 corresponding to magnetic field signatures with magnitudes of 10-100 nT, with lower current intensities occurring during less active times and the larger values during times of higher activity (Anderson et al., 2018).
The availability of an increasing amount of magnetic field data over the course of the MESSENGER mission enabled the construction of more complex empirical magnetic field models.These models are useful tools for describing the global magnetic configuration of planetary magnetospheres (Tsyganenko, 2013).The models are analytical representations of planetary magnetic fields, which are fit to collections of magnetometer observations.A common methodology for constructing these models is a modular approach.In this approach, the external magnetic field, B ext , is decomposed into a superposition of fields generated by individual sources, This introduces a set of linear amplitude coefficients, a i , which can then be fit to the available magnetometer data.The individual source fields, B i (Stern, 1994), are chosen to be of the form, B i (r;β j ), whereby the parameters, β j , define physical properties of the underlying current systems, such as the location or the thickness of a current system.These parameters are either set based on a priori knowledge or are fit to magnetometer data.Korth et al. (2015) developed a modular, empirical magnetic field model for Mercury using MESSENGER Magnetometer observations acquired while orbiting the planet for 20 months.To remove the internal contributions in the observed magnetic field vectors, the model adopted an axially aligned offset dipole to represent B int , with a magnetic moment of 190 nT R M 3 (Johnson et al., 2012) and an offset value of 0.196 R M (Anderson et al., 2012).
This model sought to analytically describe the magnetic field generated by two magnetospheric current systems, that is, the cross-tail current, B t , and the Chapman-Ferraro currents flowing on the magnetopause, B cf .B t was defined using two separate modules.Closer to the planet (r = 1.5-5RM ), a disk-like current system, B d , was used to model the cross-tail current (Tsyganenko & Peredo, 1994) while further away in the magnetotail (r > 5R M ), a solution resembling a Harris current sheet, B s , was utilized (Harris, 1962).The Chapman-Ferraro system was obtained by independently shielding these modules, along with the offset, dipolar planetary field, B int , on a static magnetopause boundary defined in Winslow et al. (2013).The total Chapman-Ferraro shielding field is thus s , and the total modeled magnetic field is given by , where t 1 and t 2 are the linear coefficients for the intensity of the disk and Harris tail current systems, respectively.
A useful property of these modular models is that they can readily account for magnetospheric dynamics by considering the linear amplitude coefficients and nonlinear parameters as functions of time, a i = a i (t), β j = β j (t).For models of the Earth's magnetosphere, a i (t) and β j (t) are parameterized using metrics known to correlate with magnetospheric activity.These include the solar wind velocity, dynamic pressure, and the interplanetary magnetic field (IMF) measured near-continuously by upstream solar wind monitors, as well as geomagnetic indices computed from data sampled by networks of ground-based magnetometers.Such approaches have allowed for detailed global empirical reconstructions of geomagnetic storms (e.g., Sitnov et al., 2008;Stephens et al., 2016) and more recently magnetospheric substorms (Stephens et al., 2019;Tsyganenko et al., 2021).Modeling the dynamics in Mercury's magnetosphere is more challenging because MESSENGER could not concurrently sample the solar wind and the magnetosphere.Further, the generation of geomagnetic indices at Earth with high temporal resolution is possible due to the widespread distribution of groundbased magnetometer stations.In contrast, during the MESSENGER era, this spacecraft was the sole source of magnetic field information at Mercury, so that the magnetospheric disturbance level could only be inferred from the MESSENGER data.In order to achieve sufficient spatial sampling in the magnetosphere to derive a measure for the disturbance level, one must use sufficiently large time windows.Following these considerations, Anderson et al. (2013) derived the Mercury disturbance index (DI) from the MESSENGER Magnetometer data on an orbit-by-orbit basis.In brief, the DI was constructed as follows.The standard deviations of B ext were computed for three different frequency bands for an entire orbital pass through the magnetosphere.These standard deviations were then summed, normalized to account for the different orbital configurations and Heliocentric distances, and then placed on a scale ranging from 0 to 100, with DI = 0 corresponding to the least active orbits and DI = 100 being the most active.The resulting index allowed Korth et al. (2017) to generalize their earlier magnetospheric magnetic field model (Korth et al., 2015) by parameterizing the linear amplitude coefficients as function of DI: t 1 = t 1 (DI) and t 2 = t 2 (DI).Further, the magnetopause boundary was also defined as a function with dependence on DI.In order to determine the functional form of this dependence, the MESSENGER Magnetometer data was binned into five DI quintiles.For each quintile, t 1 , t 2 , and the magnetopause sub-solar distance, R ss , were fit independently.The resulting fit showed a linear dependence of t 1 and t 2 on DI, see Figure 1 from Korth et al. (2017), implying a likewise linear increase of the tail current with rising magnetic activity.Additionally, the factor f = R ss r H 1/ 3 linearly decreased with increasing DI, where r H is the average Heliocentric distance of the MESSENGER spacecraft within the DI quintile.The factor r H 1/ 3 is included to account for the theoretical dependence of the magnetopause location on solar wind dynamic pressure.This finding shows that not only does Mercury's magnetopause expand and contract in response to changes in the solar wind dynamic pressure (Winslow et al., 2013;Zhong et al., 2015), which varies statistically with r H , but also the magnetosphere is more compressed during active periods, presumably because this activity is initiated by enhanced dayside reconnection (Slavin et al., 2014(Slavin et al., , 2019;;Slavin & Holzer, 1979).The final model, termed KT17, was then constructed by linearly fitting the values of t 1 , t 2 , and f across the DI quintiles.Since DI is a function of time, a time-dependent model of the external field, B ext (r,t), is obtained, which can be combined with the offset dipole model (Anderson et al., 2012;Johnson et al., 2012) to yield a dynamical model of the total magnetospheric magnetic field B tot (r,t).
Here, we seek to expand upon the modular KT17 empirical magnetic field model (Korth et al., 2017) by adding a module for Mercury's field-aligned currents.In Section 2, we describe the mathematical framework used to construct the module, adopting the approaches developed for Earth-based empirical magnetic field models (Tsyganenko, 1991(Tsyganenko, , 2002)).The fitting of the module's nonlinear parameters and amplitude coefficients is also detailed in Section 2. The resulting fit FAC model and its comparison to the magnetometer data is presented in Section 3. Finally, the results are discussed and summarized in Section 4. The source code of the FAC module has been implemented in the FORTRAN programming language and is included in the Supporting Information S1.

Field-Aligned Current Model
In this section, we detail the mathematical structure of the FAC model.First, we describe the conical current system that is the backbone of the FAC model (Section 2.1) followed by its shielding field (Section 2.2).In Section 2.3, we detail how dynamical effects, such as the response of the system to changes in the solar wind ram pressure and the level of magnetic activity, are accounted for following the approach taken in the construction of the KT17 model.The last section (2.4) then briefly describes how the FAC model is fit using the MESSENGER Magnetometer data.

Field-Aligned Current Structure
The backbone of the FAC model is a set of radially directed "conical" currents (Tsyganenko, 1991) J (N)  cone = J r r, where the superscript (N) indicates that this is the current system for the northern hemisphere.The geometry of these currents, shown in Figure 1, resembles a cone with its apex emanating from the origin (location of the magnetic dipole center) with the z axis serving as the axis of symmetry.The width of the cone is defined by the polar angle θ = θ 0 .We assume J r to have the form of a separable function J r = R(r) P(ϕ) F(θ).The conical currents are restricted to flow within a layer centered about θ 0 and of angular thickness 2Δθ, where Δθ is the halfthickness.Thus, F(θ) is non-zero when (θ 0 Δθ) < θ < (θ 0 + Δθ), and F(θ) = 0 applies outside this layer.Imposing current continuity (∇ • J = 0) requires R(r) = r 2 .The functional dependence of the azimuthal angle is given by a Fourier expansion: where ϕ is the magnetic local time.Considering the typical dawn-to-dusk asymmetry of the FAC system (Cowley, 2000), only the sine terms are used.Here, utilizing only the first Fourier term yields: (1) Consistent with Equation 1, the magnetic field, B(r, θ, ϕ), can be found by identifying a vector potential that satisfies μ 0 J = ∇ × ∇ × A and thus B = ∇ × A which also ensures a divergenceless magnetic field.The vector potential solution found by Tsyganenko (1991) and adopted here takes the form: A (N) cone = T(θ) sin ϕr.This allows F(θ) to be represented as: Note that finding an analytical expression for T(θ) that satisfies the constraint that F(θ) is non-zero within the current layer θ 0 ± Δθ and zero outside of it is not straightforward, and we refer to the prior work for details.The expression for T(θ) was found to be The result is that the magnetic field due to the conical current system takes the form: Equation 4 only yields a solution for the northern hemispheric FACs, but the southern FACs are readily modeled by mirroring the solution in the x-y plane (Tsyganenko, 2002).In principle, the northern and southern hemispheres could be treated independently, but MESSENGER's strongly elliptical orbit results in a dearth of lowaltitude data in the southern hemisphere so that the FACs in the southern hemisphere cannot be independently determined and are thus assumed to be a mirror of the northern FACs.An x-y plane reflected field is readily modeled by: The total conical current field is then the sum of the north and south hemispheric contributions: cone , and the spatial structure of the FACs is determined by two free nonlinear parameters, θ 0 and Δθ.We note, the conical current system described here assumes the FACs flow radially and are thus not aligned to the orientation of the magnetic field.This assumption is likely only valid at low-altitudes and becomes increasingly inaccurate at higher altitudes.Earth based models of the FACs account for this by bending the conical currents so they approximately flow along dipolar field lines (Tsyganenko, 2002).However, as detailed below, the limitations of the MESSENGER orbit result in substantial challenges for determining the FACs' 3D structure and their closure paths from the MESSENGER Magnetometer data.

Shielding the Conical Current System
As with many empirical magnetic field models (e.g., Tsyganenko, 2013), the KT17 model adopted a closed magnetosphere.In a closed magnetosphere, magnetic field lines cannot cross the magnetopause boundary.This characteristic is imposed by enforcing that the normal component of the total magnetic field on the magnetopause surface, S, be zero to within numerical limitations: B tot • n| S = 0.In order to close the magnetic field lines of the conical current system, a shielding field, B (sh)  cone , is introduced.Following the typical procedure for shielding fields, the domain of validity for the model is limited to just within the magnetopause boundary, where magnetopause currents are zero, thus allowing B (sh)  cone to be irrotational and written as the gradient of a magnetic scalar potential: B (sh)  cone = ∇U (sh) cone .Given the divergenceless nature of B, U (sh) cone can be represented by one of the many general solutions to Laplace's equation.In Cartesian coordinates, a solution to U (sh)  FAC is an expansion of "Box" harmonics: where a ik are linear coefficients, p i and p k are nonlinear parameters, and the coordinates are in units of R M (e.g., Tsyganenko, 2013).The parameters p i and p k were assigned the static values p i = 0.2i, p k = 0.2k.N 2 determines the number of expansions, here set to 25 (N = 5).
The values for a ik are found by randomly sampling locations on S and minimizing the root-mean-square of the normal component of the combined field, on that surface, where r i is the position of the ith sample of S, and n i is normal to S at that position.We randomly selected I = 25,000 samples on the modeled magnetopause out to a distance of r = 12 R M , taken to be the size of Mercury's magnetosphere, and utilized the singular-value-decomposition pseudo-inversion method (Press et al., 1992) to solve for a ik .Thus, the domain over which the model is valid thereby extends to x MSM = 12 R M down the tail.Note, the non-linear parameters p i and p k could have been included as free parameters when minimizing Equation 7, however, this would make the solution of a ik non-unique.The prescribed values are comparable to the exact solution of a boundary-value problem where the length of the system is 15 R M .
The magnetopause model adopted here is identical to that used in the KT17 model.Based on the terrestrial magnetopause model of Shue et al. (1997), the functional form of the magnetopause boundary is given by: The axis of symmetry for the model magnetopause is the planet-Sun line, here defined as originating at the center of Mercury's dipole, which is displaced 0.196R M northward of the planet's geometrical center, and the center of the Sun after accounting for the aberration of the solar wind (Korth & Stephens, 2022).The distance between the dipole center and the magnetopause surface, r, is parameterized by the angle between the planet-Sun line and the magnetopause surface, ε. α is the flaring parameter and R SS is the average subsolar standoff distance.An analysis of Mercury magnetopause crossings found α = 0.5 and an average value of R SS = 1.42RM (Korth et al., 2015;Winslow et al., 2013).The resulting FAC model for the average sized magnetosphere is where the coefficient, A FAC , determines the magnitude of the current in the conical current system.The resulting FAC model thus includes two free nonlinear parameters (θ 0 and Δθ), an amplitude coefficient (A FAC ), as well as 25 coefficients (a ik ) for its shielding field.Also note, the naught subscripts used in the left-hand-side of Equation 9signify that the FAC model is for the static model and does not yet account for dynamical effects that may expand and contract the system as discussed below.

Dynamical Effects
Mercury's magnetosphere is embedded within a highly dynamic solar wind and the IMF it carries.Variations in the solar wind and IMF drive large-scale changes in planetary magnetospheres.For example, at the magnetopause boundary, to first order the magnetic pressure of the planetary field is in balance with the ram pressure, P ram , in the magnetosheath.This forces Mercury's magnetosphere to compress as P ram increases and, inversely, to expand as P ram decreases (Winslow et al., 2013).The solar wind also initiates a reconnection-driven Dungey cycle (Dungey, 1961) within Mercury's magnetosphere (Slavin et al., 2021).The Dungey cycle begins when the IMF reconnects with Mercury's magnetic field along the dayside magnetopause.The flow of the solar wind then drags the newly connected field lines, eroding the dayside magnetosphere (Heyner et al., 2016;Slavin et al., 2014;Slavin & Holzer, 1979) while transporting plasma, magnetic flux, and energy into the magnetotail (Sun et al., 2015).Magnetic flux accumulates in the lobes, squeezing the central plasma sheet while intensifying the cross-tail current, before initiating tail reconnection at x MSM ≈ 3R M (Poh et al., 2017), generating a Mercury substorm (Imber & Slavin, 2017;Siscoe et al., 1975).Tail reconnection triggers fast sunward flows (Dewey et al., 2020), which dipolarizes the tail magnetic field and returns magnetic flux back to the dayside magnetopause, completing the Dungey cycle.During times of extremely high values of P ram along with an intense and southward IMF, Mercury's dayside magnetosphere can disappear entirely (Slavin et al., 2019).
The solar wind sweeps past Mercury's magnetosphere in less than a minute, assuming a typical velocity of 400 km/s.Meanwhile, the Dungey cycle and substorm timescales are on the order of 2 min (Imber & Slavin, 2017;Siscoe et al., 1975;Slavin et al., 2010Slavin et al., , 2021;;Sun et al., 2015).As noted previously, these short timescales and the lack of concurrent upstream observations of the solar wind, makes it unfeasible to empirically model specific events.However, their aggregate statistical effects can be incorporated (Korth et al., 2017).Here, we briefly summarize the methodological approach for how these dynamical effects were built into the KT17 model, and, following the same approach, how they are accounted for in this FAC model.
The expansion and contraction of Mercury's magnetopause can be modeled by treating the free parameters of Equation 8, R SS and α, as functions of P ram .Empirical studies found a clear dependence of R SS on P ram but little to no dependence on the flaring parameter α except for times with extreme driving (Winslow et al., 2013).Thus, α can be treated as a constant in time resulting in a modeled magnetopause that expands and contracts in a "selfsimilar" fashion, that is, its size changes without varying its shape.Because concurrent observations of the solar wind were not available while MESSENGER was in the magnetosphere, P ram is unknown.However, equating the ram pressure with the magnetic pressure of a pure dipole provides a theoretical dependence of R SS ∝ P ram and R SS is thus f = R SS r H 1 3 .Generally, f is not a constant in time because other factors impact R SS .For example, P ram is highly dynamic due to the natural variability of the solar wind which modulates the location of pressure balance of the magnetopause.Additionally, dayside reconnection decreases R SS (Slavin et al., 2014), and currents induced in the planet's core may generate magnetic fields that strengthen the planetary field and increase R ss (Johnson et al., 2016).The aggregate of these effects, averaged over the several hours it takes the MESSENGER spacecraft to traverse through the magnetosphere, are presumably captured by the aforementioned magnetic disturbance index, DI(t) (Anderson et al., 2013).
To determine the dependence of f on DI, Korth et al. (2017) identified the inbound and outbound magnetopause crossings for each orbit.R SS was then evaluated for each crossing using Equation 8.The magnetopause crossings were then binned into five DI quintiles, and the average values of R ss and r H were computed across each quintile.The result indicated a linear decrease in f as a function of DI.These values were fit linearly to obtain a continuous dependence f = (2.06870.0028 • DI).To summarize, the modeled magnetopause surface used in the KT17 model, which is also used here, is given by Equation 8 with a constant flaring parameter α = 0.5 and a dynamic subsolar standoff distance R SS = (2.06870.0028 As is common with most empirical magnetic field models (e.g., Tsyganenko, 2013), the KT17 model represents all magnetospheric current systems using the same self-similar rescaling property as the magnetopause.The primary reason for this treatment is to avoid complications with the derivation of their shielding fields.For example, if the magnetopause and tail current expand and contract together, then the tail current shielding field can be developed for single nominal magnetosphere size.This then allows both systems to be expand and contract together by simply rescaling the position vector, κr MSM = r 0 , where the position vector in the physical coordinate system, r MSM , is related to the position vector in the nominal coordinate system, r 0 , via the scalar κ. κ is then the ratio between the nominal standoff distance and the standoff distance in the physical coordinate system κ = R 0 / R SS , where R 0 = 1.42RM (Korth et al., 2015;Winslow et al., 2013).To keep the analytical structure of the FAC model consistent with the KT17 model, the rescaling of the position vector is retained, thus the FAC model in the MSM coordinate system is In addition to the varying magnetopause, the KT17 model was also made dynamic by parametrizing the linear amplitude coefficients that determine the relative strength of the disk current, t 1 , and quasi-Harris tail current modules, t 2 , as a function of DI.This was accomplished by binning the MESSENGER Magnetometer vector data into five DI quintiles, with 0 ≤ DI < 20 being the least active and 80 ≤ DI ≤ 97 being the most active quintiles.Note, data from the most active 3% of orbits (DI > 97) were excluded from the fit.These time periods correspond to extreme solar wind driving during which some assumptions, for example, a time-constant flaring parameter α, are no longer valid.The coefficients t 1 and t 2 were independently fit to the data contained within each quintile.
The results indicated a linear increase in t 1 and t 2 as a function of DI.Linear fits were performed for these t 1 and t 2 values across the DI quintiles to obtain their continuous dependence.Detailed in the next section, we apply a similar methodology for determining the dynamics of the above-described FAC model.The primary difference is that, in the construction of the KT17 model, each module's linear coefficients were made dynamic, while here we treat both the FAC amplitude coefficient, A FAC , along with the nonlinear parameters, θ 0 and Δθ, as functions of DI.As shown below, this allows the FAC system to expand to lower latitudes with increasing magnetic activity in contrast to the contracting magnetopause.Thus, the FAC conical current system and its shielding field are not self-similar with one another.We mitigate this issue by also treating the FAC shielding coefficients, a ik , as linear functions of DI.

Fitting the FAC Model to MESSENGER Magnetometer Data
The MESSENGER spacecraft orbited Mercury from March 2011 until April 2015.During this time, the MESSENGER Magnetometer (Anderson et al., 2007) collected three-axis measurements of the magnetic field nearly continuously.To create a model of Mercury's magnetospheric magnetic field, these data were averaged in 1-min intervals and filtered to remove observations when the spacecraft was outside of the magnetopause.The earlier static model (Korth et al., 2015) was fit to orbital data from the period 24 March 2011 through 28 November 2012.The dataset was expanded through the end of the mission for the development of the KT17 model.The fit magnetometer vectors were represented in the MSM coordinate system and corrected for solar wind aberration.The aberration correction accounts for lagging motion of the magnetosphere during the planet's orbit around the Sun.Thus, in aberrated MSM coordinates, the +x axis does not point toward the Sun but instead antiparallel to the apparent solar wind flow.Specifically, the aberration correction is implemented by rotating the MSM coordinate system by the aberration angle, χ = tan 1 (v M /v SW ) where v M is the tangential component of Mercury's velocity as it orbits the Sun (47 km/s average) and v SW = 400 km/s, about the MSM +z axis.The smaller impact of Mercury's non-zero radial velocity is ignored.The data was then binned into five approximately equally sized DI quintiles.For example, data from the least active 20% of orbits (0 ≤ DI < 20) were placed into the first quintile, the next 20% of orbits in the second quintile, and so on.This binned dataset is again used here to fit the FAC model.The MESSENGER Magnetometer dataset along with further details on the coordinate systems and the aberration correction procedure are archived at the Planetary Data System (Korth, 2022;Korth & Stephens, 2022).
The FAC model described above, B FAC (r MSM ,t), contains two nonlinear parameters (θ 0 , Δθ), that determine its spatial location and size, an amplitude coefficient, A FAC , that defines the intensity of the system, and a set of 25 coefficients, a ik , for the shielding field U (sh) FAC , all of which have to be determined.In order to isolate the impact of the FACs, the KT17 model, including the internal field, are subtracted from each observation contained within the MESSENGER Magnetometer dataset δB n = B mag,n -B KT17,n (r n ).Here, B mag,n is the nth observed magnetic field vector within a DI quintile; B KT17,n (r n ) is the modeled KT17 field evaluated at the location of the MESSENGER spacecraft, r n , during that nth observation; and δB n is the resulting KT17 residual for that nth observation.For each DI quintile, θ 0 , Δθ, A FAC , and a ik are independently fit by minimizing the B RMS residuals: Here, δB i,n are the component-wise KT17 model residuals for the nth observation in that activity quintile; B FAC (r n ) is the modeled FAC field evaluated at the location of the MESSENGER spacecraft, r n , during that nth observation; and N then represents the total number of MESSENGER observations in the respective DI quintile.
A FAC is determined by minimizing Equation 11 using linear least squares and the shielding coefficients, a ik , are fit using the procedure described in Section 2.2.The nonlinear parameters, θ 0 and Δθ, are identified by minimizing Equation 11 using the Nelder-Mead downhill simplex method (Nelder & Mead, 1965).In order to solve the nonlinear parameters and linear coefficients concurrently, the linear solvers are nested within the Nelder-Mead algorithm.That is, when the Nelder-Mead solver adjusts the value of a nonlinear parameter, a new set of corresponding shielding coefficients are found.Then, given that set of θ 0 , Δθ, and a ik , A FAC is then resolved.The resulting value for Equation 11 is communicated back to the Nelder-Mead algorithm, which continues to adjust the nonlinear parameters until the solution converges or the maximum number of iterations is reached.During the course of this investigation, a variety of initial values for θ 0 and Δθ were tested, including those approximating the distributions of Anderson et al. (2014Anderson et al. ( , 2018)).Regardless of the set of initial values, the fit values consistently converged to those shown in the next section.The final set of values was found by restarting the simplex close to the minimum location (Press et al., 1992), and allowing the fit to reach 250 iterations by setting a small convergence tolerance of 1.0 × 10 8 .

Results
The fits for θ 0 , Δθ, and A FAC for each of the five magnetic activity bins are displayed in Figure 2. Except for the highest activity bin for A FAC , they show a linear relationship on DI indicating that the FACs expand to lower latitudes and increase in intensity with increasing magnetospheric magnetic activity.As discussed above, allowing the shape of the FAC systems to change in time violates the self-similarity assumption complicating the shielding field as there are separate sets of shielding coefficients for each DI bin.To rectify this, the shielding coefficients, a ik , were also fit independently for each activity bin.Fortunately, these fits also exhibit a clear linear dependence on DI (similar to those shown in Figure 2), the values of which can be found in the Supporting Information S1.This allows θ 0 , Δθ, A FAC , and a ik to be made functions of DI(t) by simply linearly fitting the values in each activity bin as is shown in Figure 2, thus giving the complete model B FAC (r MSM ,t).
In order to demonstrate the effectiveness of the FAC model in capturing information contained within the magnetometer dataset, the RMS residuals, evaluated using Equation 11, for each of the five DI activity bins are listed in Table 1.Here and for the remainder of this section, the computation of residuals is limited to data contained within the northern hemisphere (λ ≥ 0°N) and low-altitudes (≤1,000 km) because this is the region in the magnetosphere where the FAC densities in MESSENGER's orbit are largest and thus affect the model most significantly (λ is the body-fixed latitude and altitude is defined as the height above the planet surface assuming a constant radius, R M = 2,440 km).Note however, this low-altitude filter was not applied when fitting the FAC model, that is, the results shown in Figure 2.
For the shielded dipole model, the residuals increase approximately linearly from 26.76 to 35.96 nT for low (0 20) to moderately high (60 80) activity, and to 51.69 nT for the highest activity bin.The KT17 disk and current sheet magnetic field modules along with their shielding fields reduce the residuals by 11.7% across all bins, but their impact is relatively larger for the lower activity quintiles: a 12.6% reduction for the low (0 20) and a 7.6% reduction for the highest activity (80 100) quintiles.The addition of the FAC model further reduces the residuals by another 8.1% over the KT17 model alone.Now, the relative reduction in residuals increases with increasing activity: a 5.0% and 9.5% reduction in the low and highest activity quintiles respectively.Overall, the reduction in residuals afforded by the FAC model is only modest but is of the same order as that of the equatorial magnetospheric currents, and its inclusion becomes more significant at higher activity levels.
To understand the spatial distribution of the magnetic residuals and associated FACs, the spatially binned KT17 residuals, δB, are compared to the output of the model in Figure 3 for the second highest magnetic activity quintile (DI = 60-80).Each spatial bin in the first column is 3°in latitude by 10°in longitude.The resulting set of records  Note.The first five rows of the table correspond to the residuals computed over each of the five DI quintiles individually, while the last row indicates the residuals evaluated across all five bins.Meanwhile, the columns, from left-to-right, demonstrate the impact additional modules have in reducing the residual values.The first column indicates the DI range over which the residuals were computed.The second column shows the RMS residuals computed using only the shielded dipole model as defined in the KT17 model.The third column displays residuals using the full KT17 model, equivalent to the δB terms in Equation 11.The last column displays the residuals using the KT17 model in conjunction with newly constructed FAC model.

Journal of Geophysical Research: Space Physics
10.1029/2023JA031702 STEPHENS AND KORTH within each spatial bin is individually sorted by their x, y, and z components, allowing the median δB x , δB y , and δB z values to be computed.These median values are converted to spherical coordinates (δB r , δB θ , and δB ϕ ), and their spatial distributions are displayed.The radial current density, J r , displayed in the bottom left panel of Figure 3, is computed using the differential form of Ampère's Law: The differential terms, ∂B ϕ ∂θ and ∂B θ ∂ϕ , are estimated using the finite difference method: For example, to evaluate J r for a target spatial bin, the values of δB ϕ for the adjacent northern and southern bins are differenced to yield ΔB ϕ , and the result is then divided by the angular separation between these two bins, which for the bin definition used here is Δθ = 6°.Similarly, ΔB θ values are computed as the difference between δB θ in the adjacent bins to the West and East using a longitudinal separation of Δϕ = 20°.θ and ϕ in Equation 12 are then the colatitude and longitude of the center of the target bin, and r is the planet's surface which equals 1 R M .To further visualize the pattern of residuals, lines have been overplotted on the current density panel indicating the direction and magnitude of the residual magnetic field vectors.
The residual distributions in Figure 3 reveal several distinct patterns.It is evident that the δB signatures best articulated by the FAC model are in the longitudinal component.The longitudinal residuals, δB ϕ , display a dominant dawn-dusk asymmetry and at mid to high latitudes are directed primarily eastward at dawn and westward at dusk.At lower latitudes (λ ≲ 30°N), the azimuthal direction reverses, indicating a curl in the magnetic field, which is synonymous for the presence of a radial current.Northward of 65°N, the residual vectors are predominantly sunward and their magnitude is large, particularly on the dayside.Southward of this latitude on the dayside, the magnitude of the magnetic residuals is smaller and their direction is antisunward.On the nightside, particularly within 2 hours of local midnight, the residual vectors remain sunward oriented at all latitudes but gradually diminish in magnitude from the pole toward the equator.A different pattern is evident along the dawndusk meridian.At higher latitudes, the residual vectors possess large magnitudes and a sunward orientation.The residual magnitudes gradually decrease from the pole to the equator and reverse their orientation to antisunward equatorward of 30°N, again indicating a curl in the residual magnetic field.
At high latitudes, the radial residuals, δB r , point predominantly away from the planet.They change sign and point toward the planet at mid-latitudes ≲60°N on the dayside and ≲70°N on the nightside.On the nightside, they are again positive at low latitudes between 10°N ≲ λ ≲ 30°N.The colatitudinal residuals, δB θ , are mostly oriented poleward except for at dayside latitudes northward of 65°N, where they point equatorward, a signature that was previously attributed to the cusp (Korth et al., 2015).
The respective model representations for B r , B θ , B ϕ , and J r are shown in the center column of Figure 3.The FAC model is evaluated on a sphere with a radius of 1 R M from the planet center using the center value of the DI quintile ranging from 60 to 80%, that is, DI = 70 and the average distance from Mercury to the Sun, r H = 0.39 AU.The bottom panel, J r , confirms that the Mercury FAC system has the same polarity as the terrestrial Region-1 FACs but exhibits a somewhat different morphology.Instead of being restricted to a ∼10°wide latitude band, the modeled FACs are spread over a large range of latitudes, with the equatorward boundary being at ≈30°N and extending nearly to the pole at 83°N.This equatorward extend is beyond the average open-closed field line boundary (Lindsay et al., 2016).The Mercury FAC distribution resembles a half-moon structure instead of the concentric ring structure as seen at Earth.This half-moon morphology of the FACs results in a distinctive magnetic field pattern as shown in B r , B θ , and B ϕ panels.The best qualitative agreement between the FAC model and the residual data is present in the B ϕ component of the magnetic field.It shows a clear dawn-dusk asymmetry, with the magnetic field being westward oriented at dusk and eastward at dawn, until it flips orientation equatorward of 35°N . This matches the pattern seen in the KT17 residual data, except that the flip in B ϕ occurs at a somewhat lower latitude (λ ≈ 30°).The conical currents have no radial component in the coordinate system of their definition (cf.Equation 4), thus the small non-zero values of B r must be imparted from the conversion from MSM to MSO coordinates and/or originate from the shielding field.On the dayside, B θ is directed entirely equatorward and nearly entirely poleward on the nightside.This corresponds to predominantly sunward orientation near the noon midnight meridian.This agrees qualitatively with the residual distribution on the nightside but deviates from the distribution on the dayside equatorward of 65°N, where the residual vectors are antisunward.The differences between the binned KT17 residual data and the FAC model are shown in the right column of Figure 3.These panels demonstrate that the FAC model is largely incapable of describing the residual B r and B θ components.However, the model matches the bulk B ϕ structure contained within the residual data, with the magnitude of the differences being less than 10 nT across the majority of the northern hemisphere.
To confirm that these patterns exist across all magnetic activity levels, the B ϕ panel (third row of Figure 3) is shown for each of the five DI quintiles as is displayed in Figure 4.It shows two trends that correlate with the activity level.First, as magnetic activity increases, the magnitude of the magnetic field residuals increase (left column), as would be expected from the numerical values in Table 1.Second, along the dawn-dusk line, the change in orientation of the residual magnetic field vectors, sunward oriented at higher latitudes to anti-sunward at lower latitudes, occurs at lower latitudes with increasing activity.In the least magnetically active bin (DI = 00-20), this occurs at 35°N, but by the highest activity bin (DI = 80-100) the change in orientation is now seen at 20°N . Thus, the FACs become more intense and extend to lower latitudes with increasing activity.This is also apparent in the FAC model output shown in the center column in Figure 4, which visualizes the linear trends previously described in Figure 2.

Discussion and Conclusions
We have developed an empirical magnetic field model for Mercury's large-scale FACs using an analytical formulation of "conical" currents (Tsyganenko, 1991(Tsyganenko, , 2002)).The model's analytical framework follows the modular approach that has been widely used to describe electric currents in the Earth's magnetosphere and their corresponding magnetic fields (Tsyganenko, 2013 and references therein).Specifically, the FAC model is a direct extension to the modular KT17 magnetic field model developed to describe the magnetic field associated with Mercury's cross-tail and magnetopause current systems (Korth et al., 2015(Korth et al., , 2017)).
Similar to the KT17 model, the FAC module is designed to be time variable in three steps.First, the magnetopause boundary used to shield the conical currents, varies in size as a function of Mercury's heliocentric distance, which reflect systematic variations in solar wind ram pressure, and as a function of the Mercury magnetic disturbance index (DI).Second, to maintain consistency with the KT17 modeling framework, the FAC model utilizes the same spatial rescaling as the magnetopause.However, unlike the KT17 cross-tail current system, the FAC distribution does not vary self-similarly with the magnetopause because these currents expand to lower latitudes with increasing magnetospheric activity.To mitigate this inconsistency, the set of FAC shielding coefficients is made linearly dependent on DI.Third, the FAC model's linear amplitude coefficient, which determines the intensity of the currents, and non-linear parameters, which prescribe the spatial extent of the system, are parameterized as function of DI.The latter dependence is created by binning the MESSENGER Magnetometer data into five approximately equally sized DI quintiles, fitting the values independently for each quintile, and then linearly interpolating across the quintiles to obtain a continuous dependence on DI.
This choice of the shielding field could result in incomplete shielding of the FAC system on the magnetopause surface.To test this, the RMS value, that is, the value for σ in Equation 7, of the normal component of the shielded FAC magnetic field is evaluated by sampling 25,000 random locations on this surface with distances of up to r = 12 R M using the center DI value for each DI quintile.The values of σ ranged from 0.03 nT for DI = 10 and up to 0.07 nT for DI = 90, corresponding to a ∼95% reduction in the unshielded values: σ = 0.65 nT for DI = 10 and σ = 1.55 nT for DI = 90.For comparison, this is much less than the observed values (>10 nT) of the magnetopause-normal component of the field near reconnection sites (Slavin et al., 2009).Figures 5a and 5b display magnetic field lines (black lines) traced from the planet's surface using the superposition of magnetic fields from KT17 and the FAC module.The efficacy of the shielding fields is demonstrated by the field lines populating the volume of the magnetosphere but not crossing the magnetopause boundary (green surface).Note, the magnetic field lines in Figure 5b do not cross the magnetopause boundary; they appear to do so because their  projection in the y-z plane is plotted, while the magnetopause boundary and current density show their y-z planar slice.
Fitting the FAC model to the MESSENGER Magnetometer dataset results in a Region-1 sense FAC system which intensifies and expands to lower latitudes with increasing magnetospheric activity.The magnitude of the modeled FAC densities reaches values of 22 nA/m 2 at the planet's northern surface during the most active times.This is a factor of ∼100 smaller than the current densities observed in Earth's FACs during storm-time substorms (Lukianova, 2020).Even during quiet conditions, the magnitude of the FACs is on the order of 10 nA/m 2 .When DI = 0, the modeled FACs extend to λ = 39.7°N and reach as low as 27.6°N at the highest activity levels (DI = 100).Across all activity ranges, the northern extent of the modeled FACs is ~83°N, which corresponds to the maximum latitudinal extent of MESSENGER's orbit and thus may be a consequence of the spatial sampling of the dataset.Qualitatively, both the magnetic field vectors and the associated current densities resemble the Region-1 FAC at Earth.However, the magnetic field residual vectors in Figure 3 show no indication of a Region-2 system, consistent with earlier findings (Anderson et al., 2014(Anderson et al., , 2018)).The center of Earth's FAC ovals are shifted several degrees toward the nightside (Carter et al., 2016;Iijima & Potemra, 1976).The B ϕ residual panels in Figure 4 indicate no such nightward shift in the FAC pattern and as such was not included in the model structure.However, the modeled distribution of FACs indicated by the radial current densities, J r , in the bottommiddle panel of Figure 3, differs considerably compared to earlier studies.For instance, also using MESSENGER Magnetometer data, Anderson et al. (2018) found much larger current densities, j r ≈ 120 nA/m 2 and j r > 300 nA/ m 2 during quiet and the most active times respectively, concentrated in a narrower range of latitudes from ∼70 to 80°N.Meanwhile, a multifluid MHD simulation of Mercury's magnetosphere during southward IMF conditions revealed an FAC distribution of J r > 100 nA/m 2 that extended equatorward to ∼35°N (Dong et al., 2019; see their Figure 2e).A limitation of the FAC model presented here is the overly rigid latitudinal dependence, F(θ), of the current density defined by Equations 1-3.That, along with effect of statistical averaging over many events, may act to smear the FACs over too large of a range of latitudes.Future models can mitigate this issue by introducing a more flexible form of F(θ), however, finding alternative solutions to Equation 3 is non-trivial (Tsyganenko, 1991).Another approach developed for Earth based models, is to use several different conical current sheets, each with their own amplitude coefficients, located at adjacent or overlapping latitudes (Sitnov et al., 2017).
The addition of the FAC module to the KT17 model only modestly reduces the RMS magnetic field residuals and this reduction is primarily in the B ϕ component.Over all activity levels, the low-altitude residuals decreased from 32.26 to 29.64 nT, an 8.1% reduction.The relative reduction in residuals is greater with increasing magnetic activity.Figures 3 and 4 indicate that this reduction in residuals is principally due to the FAC model capturing the B ϕ component of the observed magnetic field as the MESSENGER spacecraft traversed across Mercury's northern hemisphere at low-altitudes.The distributions of B r and B θ were largely not reproduced, indicating that other large-scale current systems many also be missing from the model.One possibility for the origin of these magnetic fields is their generation by induced currents.Mercury's dayside magnetosphere is highly dynamic, and sudden changes of the magnetospheric magnetic field may induce currents within the planet's interior (Heyner et al., 2021 and references therein).Other potential missing current systems include the closure of these FACs along a spherical cap near the core-mantle boundary (Anderson et al., 2018) and currents associated with a dayside-boundary layer (Anderson et al., 2011;Liljeblad et al., 2017).The band of B r residuals near the nightside equator as well as the high-latitude dayside B r and B θ residuals have been attributed to diamagnetic depressions within the plasma sheet and the magnetospheric cusp respectively (Korth et al., 2015).Another source of the residuals are higher frequency variations in the magnetospheric current systems.The dynamical effects described in Section 2.3 are parameterized by DI, which varies on timescales of 8-12 hr as prescribed by MESSENGER's orbital period.This is much longer than timescales associated with variations in the solar wind and with the Mercury substorm (Imber & Slavin, 2017) and Dungey cycles (Slavin et al., 2021).Thus, the model cannot represent the instantaneous magnetic field configuration during individual events and instead captures their aggregate effect, which limits its capacity to reduce residuals.
The FAC model developed here is quite simplistic, as it only possesses a single amplitude coefficient and two free non-linear parameters.In contrast, similar FAC models for Earth utilize as many as 16 amplitude coefficients (Sitnov et al., 2017).These additional coefficients arise from introducing additional Fourier terms in the azimuthal portion of Equation 1 and by using different conical current systems at different latitudes.Also, the currents of the FAC model flow radially, which only approximates the magnetic field orientation at low-altitudes.Earth based models mitigate this issue by bending the conical current surface to match a realistic field line geometry (Tsyganenko, 2002).The complexity of the derived Mercury FAC model is constrained by limitations of MESSENGER's orbit, which only traverses the planet at low-altitude in the northern polar region, meaning it does not adequately sample low to mid-latitudes near the planet and never samples the southern hemisphere at lowaltitudes.The constraints result in substantial challenges for determining the FACs' 3D structure and their closure paths from the MESSENGER Magnetometer data.Furthermore, the FAC model is an extension to the KT17 model which itself is a time-varying update to an earlier static model (Korth et al., 2015).In the intervening years, much work has been done to further characterize Mercury's internal field (Oliveira et al., 2015;Thébault et al., 2018;Wardinski et al., 2019Wardinski et al., , 2021) ) and magnetopause boundary (Philpott et al., 2020;Zhong et al., 2015Zhong et al., , 2020)).In light of these more recent studies, the requisite models for the internal field and magnetopause used in the KT17 and this FAC model should be revisited in the future.The two-spacecraft BepiColombo mission (Benkhoff et al., 2010(Benkhoff et al., , 2021) ) will begin orbiting Mercury in 2025.The low-altitude low-eccentricity orbit of the Mercury Planetary Orbiter along with independent measurements by its companion Mercury Magnetospheric Orbiter will augment the MESSENGER dataset and lead to a more accurate characterization of Mercury's internal and external magnetic fields.Furthermore, BepiColombo will perform simultaneous observations of the magnetosphere and the solar wind, thereby illuminating more details on how magnetospheric current systems respond to their solar wind drivers (Milillo et al., 2020).

Data Availability Statement
The processed magnetometer data, including the KT17 residuals that were used to fit the FAC model, are archived on the NASA Planetary Data System at https://doi.org/10.17189/1522653(Korth, 2022).The file defining the Mercury disturbance index as a function of MESSENGER orbit number (Anderson et al., 2013), a file of Mercury's heliocentric distance as a function of time, and the digital data used for all figures have been uploaded to a Zenodo dataset at https://doi.org/10.5281/zenodo.8383564(Stephens & Korth, 2023).The source code of the FAC model has been implemented in the FORTRAN programming language and is included in the Supporting Information S1.

Figure 1 .
Figure1.Illustration detailing the conical current system used to model Mercury's FACs.The axes correspond to the MSM coordinate system, with +x pointing toward the sun, +z oriented along the planet's rotation axis, and +y points toward dusk.The model is described in spherical coordinates: (r, θ, ϕ) corresponding to the radial distance, polar angle, and azimuthal angle respectively.The parameters that define the FACs' spatial configuration are shown: θ 0 defines the polar angle corresponding to the center of the conical current sheet while ∆θ defines its half thickness.Upward ( j r > 0) and downward ( j r < 0) flowing currents are colored red and blue, respectively, with the shading corresponding to the magnitude of the current density (darker colors have larger |j r |).

Figure 2 .
Figure2.The dependence of the field-aligned conical current sheet parameters, θ 0 (orange circles) and Δθ (green circles), as well as the amplitude coefficient A FAC (pink circles) on the magnetic disturbance index, DI.Linear fits to these points are shown by the lines and their slope-intercept equations are indicated.The fit for A FAC for the 80-100 DI bin was anomalously high and was not included while fitting its linear dependence on DI.

Figure 3 .
Figure 3.Comparison between the spatially binned KT17 residuals, δB (left column) and the FAC model (middle column) sampled at the surface of Mercury's Northern hemisphere for the medium-high activity bin DI = 060-080.The difference between the binned KT17 residuals and the model are shown in the third column.The first, second, and third rows show the spherical components of the magnetic field, B r , B θ , and B ϕ respectively.The modeled open-closed boundary is overplotted in orange.The last row displays the radial current density, j r , with the black lines indicating the direction of the magnetic field vectors.Negative/positive values of B r , B θ , B ϕ , and j r are indicated by the blue/red color bars on the right.Because these are polar plots, the grids shown as well as the magnetic field vectors are now represented in the planetcentered MSO system.

Figure 4 .
Figure 4. Comparison between the B ϕ components of the binned KT17 residuals (left column) and the fit FAC model (middle column) sampled at the surface of Mercury's Northern hemisphere for each of the five magnetic disturbance index quintiles.The panels are the same as the third row of Figure 3.

Figure 5 .
Figure 5. (a) The modeled magnetic field lines (black) and the y-component of the current density, j y , plotted in the x-z plane.The value of j y is indicated by the blue-red color bar, featuring the cross-tail current in red.(b) The modeled magnetic field lines (black) and the radial current density, j r , plotted in the y-z plane showing the FACs flowing into (blue) and out of (red) the planet.The modeled field is the total magnetic field (Dipole + KT17 + FAC) evaluated at a moderate activity level, DI = 50.0,and Mercury's average distance from the Sun, r H = 0.39AU.The average model magnetopause is shown by the green surfaces.
Beard, 1964).Additionally, P ram diminishes by a factor r H 2 (P ram ∝ r H 2 ) with increasing Heliocentric distance, r H . Thus, to first order, R SS can be parameterized by r H : R SS ∝ r H &

Table 1
Modeled RMS Residuals (nT)With and Without the FAC Model for Low-Altitude Data (λ ≥ 0°N and Altitude ≤ 1,000 km)