Characterising mesoscale magnetopause surface waves within magnetosphere–ionosphere–ground coupling

Disturbances to the magnetopause location driven by upstream pressure variations or flow shear instabilities may be described as surface waves, which act as localised sources of field-aligned currents coupling the magnetosphere to the ionosphere. However, their impacts on the ionosphere and ground across representative ranges of wave and system properties are poorly understood. We, therefore, develop a simple numerical model for dispersionless mesoscale magnetopause surface waves within the coupled magnetosphere–ionosphere–ground system to gain insight into how their amplitudes and spatial scales throughout the system might vary with conditions. In general, the impacts of finite wave packets can be decomposed into periodic fluctuations (with matching wavelength to that directly above in the magnetosphere) along with slowly-varying trends that result from finite wave effects. Finite wave packets act in the far-field like a string of alternating field-aligned currents well described both in the ionosphere and on the ground as a two-dimensional current dipole. In the ionosphere, near-field periodic fluctuations exponentially decay over the reduced wavelength latitudinally away from the projected magnetopause boundary layer flux tubes, which may limit how well they can be resolved by radar. The relationship between the magnetic field above and below the ionosphere becomes more complicated for surface waves than infinite plane Alfvén waves due to the additional spatial structure, which introduces interference across the spectrum of wavenumbers present. This modifies how the ionosphere screens, rotates, and spatially smears magnetic field perturbations across all three components in different ways. For mesoscale wavelengths this importantly results in latitudinal scales of amplitude and polarisation variation smaller than typical ground magnetometer spacings, motivating the need for denser networks. A range of effective skin depths in the ground are applicable to surface waves, meaning ground induction can vary between a near-perfect insulator to a good conductor, affecting both observable ground magnetic fields and resulting geoelectric fields. The predicted peak amplitudes of surface waves' impacts suggest they may act as significant sources of ionospheric/thermospheric Joule heating and geoelectric fields in the ground, thereby contributing to space weather impacts. These are, however, highly localised latitudinally when considering typical mesoscale waves. Our results provide key insight into interpreting ground-based observations, of particular timeliness with the rollout of new digital ionospheric radars and the SMILE mission's planned conjugate ground–space campaigns.

1 Introduction

Field-aligned currents (FACs) are the main physical mechanism mediating magnetosphere–ionosphere (MI) coupling (e.g. Ganushkina et al., 2018), dynamically connecting/transferring energy to and from these geospace regions. Magnetospheric ultra-low frequency (ULF, ∼0.1–100 mHz) waves contribute to variability in this coupling and are often described using magnetohydrodynamics (MHD; Jacobs et al., 1964; Wright et al., 2024). In basic theory the only MHD body wave (free to propagate through the bulk plasma) to exhibit FACs is the shear Alfvén mode. However, discontinuities separating different plasmas allow surface waves (Archer et al., 2024b, and references therein), collective evanescent fast magnetosonic modes that support FACs within the plasmas' interface (Plaschke and Glassmeier, 2011). Magnetopause surface modes at the solar–terrestrial interaction's interface impose a wave pattern of plasma displacements on the boundary. When the wavelength and amplitude of these displacements is larger the boundary layer’s thickness, a single coupled surface mode across the entire boundary layer exists rather than distinct modes at the different interfaces (Kivelson and Pu, 1984). This surface wave compresses/rarefies the magnetic field on either side of the equatorial magnetopause in opposite senses and by different amounts, locally enhancing/diminishing the Chapman-Ferraro current. These perturbation currents within the closed field lines of the boundary layer are diverted, preserving current continuity, from largely perpendicular to magnetospheric field lines to along them, becoming FACs which then close in the ionosphere as longitudinally periodic currents. An illustration of this is shown in Fig. 1. Therefore, magnetopause surface waves directly influence the ionosphere and ground (Kivelson and Southwood, 1988; Archer et al., 2023), unlike other compressive waves which do so secondarily via coupling to Alfvén waves due to non-uniformities in the bulk plasma/magnetic field (Southwood and Kivelson, 1990, 1991).

Figure 1 Illustration of the current patterns caused by a magnetopause surface wave.

Magnetopause surface waves play a global role in filtering, accumulating, and guiding turbulent solar wind driving (Kivelson and Chen, 1995; Archer et al., 2024b). They form on the boundary maintaining total pressure balance, excited either by upstream pressure variations (originating from the solar wind, foreshock, or magnetosheath; e.g. Sibeck et al., 1989; Shue et al., 2009) or grow from seed fluctuations due to velocity shear instabilities (Kelvin–Helmholtz instability, KHI; e.g. Fairfield et al., 2000). The surface wave dispersion relation dictates how the boundary responds to different driving regimes. While the general dispersion relation can only be solved numerically, if incompressibility is assumed it has analytic solution

$$\begin{array}{rl}{\displaystyle }\mathit{\omega }=& {\displaystyle \frac{{\mathit{k}}_t\cdot ({\mathit{\rho }}_{\mathrm{0},\text{msh}}{\mathit{u}}_{\mathrm{0},\text{msh}}+{\mathit{\rho }}_{\mathrm{0},\text{msp}}{\mathit{u}}_{\mathrm{0},\text{msp}})}{{\mathit{\rho }}_{\mathrm{0},\text{msh}}+{\mathit{\rho }}_{\mathrm{0},\text{msp}}}}\\ {\displaystyle }& {\displaystyle }\pm \sqrt{\begin{array}{l}{\displaystyle \frac{{\mathit{\rho }}_{\mathrm{0},\text{msh}}({\mathit{k}}_t\cdot {\mathit{v}}_{\mathrm{A},\text{msh}}{)}^{\mathrm{2}}+{\mathit{\rho }}_{\mathrm{0},\text{msp}}({\mathit{k}}_t\cdot {\mathit{v}}_{\mathrm{A},\text{msp}}{)}^{\mathrm{2}}}{{\mathit{\rho }}_{\mathrm{0},\text{msh}}+{\mathit{\rho }}_{\mathrm{0},\text{msp}}}}\\ -{\displaystyle \frac{{\mathit{\rho }}_{\mathrm{0},\text{msh}}{\mathit{\rho }}_{\mathrm{0},\text{msp}}}{({\mathit{\rho }}_{\mathrm{0},\text{msh}}+{\mathit{\rho }}_{\mathrm{0},\text{msp}}{)}^{\mathrm{2}}}}[{\mathit{k}}_t\cdot ({\mathit{u}}_{\mathrm{0},\text{msh}}-{\mathit{u}}_{\mathrm{0},\text{msp}}){]}^{\mathrm{2}}\end{array}}\end{array}$$

for wave angular frequency ω, transverse wave vector k_t, velocities u, mass densities ρ, and magnetic fields B, where subscript 0's indicate equilibria, “msp” is the magnetosphere, and “msh” the magnetosheath (Chandrasekhar, 1961; Pu and Kivelson, 1983). Given typical conditions at the magnetopause, a simpler approximate dispersion relation

$$\begin{array}{}\text{(1)}& {\displaystyle }\mathit{\omega }& {\displaystyle }\simeq {\mathit{k}}_t\cdot {\mathit{u}}_{\mathrm{0},\text{msh}}\pm {\mathit{\omega }}_{\mathrm{0}},\left\{\begin{array}{l}{u}_{\mathrm{0},\text{msh}}\gg u_{\mathrm{0},\text{msp}}\\ {\mathit{\rho }}_{\mathrm{0},\text{msh}}\gg {\mathit{\rho }}_{\mathrm{0},\text{msp}}\\ B_{\mathrm{0},\text{msh}}\ll B_{\mathrm{0},\text{msp}}\end{array}\right.\text{(2)}& {\displaystyle }{\mathit{\omega }}_{\mathrm{0}}& {\displaystyle }\equiv \sqrt{{\displaystyle \frac{({\mathit{k}}_t\cdot {\mathit{B}}_{\mathrm{0},\text{msp}}{)}^{\mathrm{2}}+({\mathit{k}}_t\cdot {\mathit{B}}_{\mathrm{0},\text{msh}}{)}^{\mathrm{2}}}{{\mathit{\mu }}_{\mathrm{0}}({\mathit{\rho }}_{\mathrm{0},\text{msp}}+{\mathit{\rho }}_{\mathrm{0},\text{msh}})}}}\simeq {\displaystyle \frac{|{\mathit{k}}_t\cdot {\mathit{B}}_{\mathrm{0},\text{msp}}|}{\sqrt{{\mathit{\mu }}_{\mathrm{0}}{\mathit{\rho }}_{\mathrm{0},\text{msh}}}}}\end{array}$$

holds when $({\mathit{k}}_t\cdot {\mathit{u}}_{\mathrm{0},\text{msh}}{)}^{\mathrm{2}}{\mathit{\rho }}_{\mathrm{0},\text{msp}}/{\mathit{\rho }}_{\mathrm{0},\text{msh}}\ll {\mathit{\omega }}_{\mathrm{0}}^{\mathrm{2}}$ (Archer et al., 2024b). Equation (1) consists of a natural surface wave frequency for no flow shear due to standing field-aligned structure, ω₀ (Chen and Hasegawa, 1974), along with an advective Doppler shift by the background magnetosheath flow. Along most of the magnetopause, frequencies are dictated by the magnetosheath velocity ($u_{\mathrm{0},\text{msh}}\gg {\mathit{\omega }}_{\mathrm{0}}/k_{\parallel }$) and perpendicular wavelengths ($k_{⟂}\gg k_{\parallel }$), with negligible natural frequency (${\mathit{\omega }}_{\mathrm{0}}\ll |{\mathit{k}}_t\cdot {\mathit{u}}_{\mathrm{0},\text{msh}}|$) (Miura and Pritchett, 1982; Kozyreva et al., 2019). Typical equatorial wavelengths ∼2–14 R_E (azimuthal wavenumbers m∼4–30 for exp (imϕ) dependence) and periods ∼1–7 min are observed (Lin et al., 2014). However, around the subsolar point the magnetosheath is slow ($u_{\mathrm{0},\text{msh}}\ll {\mathit{\omega }}_{\mathrm{0}}/k_{\parallel }$), hence the natural frequency dominates (${\mathit{\omega }}_{\mathrm{0}}\gg |{\mathit{k}}_t\cdot {\mathit{u}}_{\mathrm{0},\text{msh}}|$) leading to a resonant magnetopause surface eigenmode (MSE) response to (particularly impulsive) external driving (Archer et al., 2019). This constitutes the lowest-frequency eigenmode of the dayside magnetosphere, typically less than 2 mHz (Plaschke et al., 2009a; Archer and Plaschke, 2015; Hartinger et al., 2015b). While MSE from solar wind forcing are global in scale (Hartinger et al., 2015b; Archer et al., 2021), bow shock/magnetosheath transients of several R_E in size (m∼20–60) also cause MSE (Archer et al., 2019). In the absence of continuous driving or instabilities, surface modes are strongly damped, likely resulting in only a few cycles at a time (Chen and Hasegawa, 1974; Kozyreva et al., 2019).

Magnetopause surface waves penetrate the magnetosphere over evanescent e-folding scale $k_t^{-\mathrm{1}}$, the reduced/angular wavelength. Thus global and mesoscale waves significantly affect geospace and subsequently couple to other modes (Horvath and Lovell, 2021; Archer et al., 2021, 2022). This leads to mass, momentum, and energy transport into geospace, contributing to the viscous-like solar–terrestrial interaction in absence of global magnetic reconnection (Axford, 1964). Understanding magnetopause surface waves' energy pathways, particularly with regards to the ionosphere and ground, is therefore important, especially within the context of space weather. It is being appreciated more that ULF waves significantly contribute to ionospheric energy input (e.g. Hartinger et al., 2015a). This results in ionospheric/thermospheric dissipation via Joule heating (e.g. Shi et al., 2025a), and induced geoelectric fields (e.g. Hartinger et al., 2020; Shi et al., 2022, 2025b) that drive currents in power grids (e.g. Belakhovsky et al., 2019; Heyns et al., 2021).

Theoretically ULF waves' effects in the ionosphere and on the ground have focused on Alfvén waves. Seminal analytic results on the ionosphere's effects of rotating ground magnetic field fluctuations by 90° (Hughes, 1974; Hughes and Southwood, 1974) and screening small-scale (compared to the ionospheric altitude) waves (Hughes and Southwood, 1976) all rely on the assumption of plane Alfvén waves vertically incident on a homogeneous ionosphere. Given limitations with analytic theory, later works have simulated Alfvén wave transmission using the linearised wave equations in more representative setups, e.g. incorporating dipole magnetic geometries (Lysak, 2004), improved boundary conditions (Lysak and Song, 2006), asymmetric ionospheric conductances (Lysak et al., 2020), or including multiple altitude-varying ion and neutral species (Sydorenko and Rankin, 2012). Theoretical models of the MI-response to upstream pressure variations have not directly invoked surface modes on the magnetopause, instead considering fast–Alfvén coupling Earthward of the boundary (Sibeck, 1990; Southwood and Kivelson, 1990, 1991). The direct effects of magnetopause surface waves have only recently been appraised through high-resolution general-purpose global MHD simulations (Archer et al., 2023). This pertained to impulsively-excited large-scale azimuthally-stationary MSE across the dayside, which seeded tailward propagating surface waves in the flanks (Archer et al., 2021). Surface wave FACs at the MI-interface peaked at the equatorward edge of the magnetopause boundary layer (MPBL), as depicted in Fig. 1, rather than the open-closed boundary (OCB). Earthward of the MPBL, weaker FACs decaying with distance occurred via fast–Alfvén coupling which exhibited slow poleward phase motion due to damping. The result in the ionosphere was convection vortices circulating FAC structures. The ground magnetic field response was predominantly, but not entirely, controlled by Hall currents. This is because Fukushima's (1976) theorem of the perfect cancellation of field-aligned and Pedersen current effects below the ionosphere only holds for semi-infinite FACs along straight vertical field lines (for either a planar or spherical ground geometry) and uniform ionospheric Pedersen conductance. The ground magnetic field was rotated from the magnetosphere by close to 90° on average, but with significant angular spread interpreted as due to violation of Fukushima's (1976) theorem. While these basic features are expected for any process generating FACS, e.g. Alfvénic field line resonances (Greenwald and Walker, 1980), and have long been associated with magnetospheric dynamics in observations (e.g. Friis-Christensen et al., 1988; Glassmeier and Heppner, 1992; Glassmeier, 1992; Bristow et al., 1995), it is the FACs' structure and the scales of variations which might enable discrimination between different dynamical phenomena. The effect of induction in the ground and thus whether surface waves specifically lead to significant geoelectric fields remains an open question (Archer et al., 2023).

Evidence for magnetopause surface modes from ground-based instruments alone was inconclusive for a long time (Pilipenko et al., 2017, 2018; Kozyreva et al., 2019), largely due to a lack of robust theoretical expectations and the scarcity of ground–space conjunctions during clear surface wave events (Archer et al., 2019). The Archer et al. (2023) simulation results demonstrated qualitative agreement with previous candidate ground-based events (Archer et al., 2019; Kozyreva et al., 2019; He et al., 2020; Hwang et al., 2022). Key aspects could not be tested with the data presented though, motivating both reanalysis of these events and dedicated future observational studies. While global MHD models self-consistently reproduce magnetopause surface waves within a representative environment, quantitative predictions from any single run strictly apply to the specific driver, set of upstream conditions, and geospace plasma configuration used, caveated by the numerics' effects on results (e.g. Zhang et al., 2019; Archer et al., 2024b). More widely-applicable comparisons between theory and observations require a comprehensive characterisation of expected surface wave responses in the ionosphere and on the ground across broad representative ranges of both wave and system conditions. However, many of these are not easily tuneable in such simulations. Additionally, the resolution required to resolve the wave physics throughout geospace make these simulations computationally expensive, hence impractical for parameter searches, rendering results more generally qualitative. Dedicated linear ULF wave simulations with dipole or more generalised geometries also exist which are more configurable (Degeling et al., 2010; Wright and Elsden, 2020; Lysak et al., 2020). Magnetopause surface waves have yet to be simulated in a realistic manner within these codes (Elsden et al., 2025), since such waves serve as the imposed outer boundary condition requiring a priori knowledge of the solution. Given the basic physical insight analytic MHD applied to highly simplified environments have provided (e.g. Southwood, 1974; Chen and Hasegawa, 1974; Hughes and Southwood, 1974), we opt to develop a simple local model to advance understanding of how magnetopause surface waves' effects within magnetosphere–ionosphere–ground (MIG) coupling might scale spatially and in amplitude across key wave and system parameters. We introduce the model and its numerical implementation in Sect. 2, derive some analytic relations in Sect. 3, and then compare these with numerical results in Sect. 4.

2 Method

2.1 Model

Using simple “toy” models initially in physics can provide physical insight and intuitive understanding to a problem, as well as enabling computational efficiency, serving as foundational tools for framing hypotheses, testing limiting cases, and bridging the gap to complex simulations. While not realistic by design, their complexity can be iterated through a step-by-step reduction of assumptions once the physical effects that are important are better understood from the simple model results. Thus “toy” models serve as an essential first step in addressing a new science question.

There are few geometries in which linear MHD surface waves can be derived analytically (Archer et al., 2024b); even those with orthogonal curvilinear coordinates (such as a dipole or cylinder) can only be solved when invoking several assumptions and still yield complicated expressions (Nenovski et al., 2007; Leonovich and Kozlov, 2019). Therefore, we opt to use one of the simplest geometries, the Cartesian box model that has long been leveraged for initial physical insights into ULF waves (e.g. Southwood, 1974) and is often still applied to surface wave problems (Plaschke and Glassmeier, 2011; Archer et al., 2023; Elsden et al., 2025). Employing the electric field and current paradigm, common of the ionosphere (Vasyliunas, 2001; Laundal et al., 2015), we solve a boundary value problem using current continuity through different conducting layers to determine the spatial structure of solutions through the coupled MIG-system given imposed magnetopause currents due to a mesoscale surface wave. In order to make our results more general and independent of any specific background, we only incorporate current systems due to the surface wave in this model, hence all responses presented are solely due to waves. The classical Region 0/1/2 FAC systems and their ionospheric counterparts serve as a time-constant (though spatially varying) background to the linear equations solved, which do not affect the wave response and can easily be separated out. Our setup is illustrated in Fig. 2a, along with the various physical quantities the model outputs applied to a single example (panels b–r).

Figure 2 (a) Diagram of the surface wave magnetosphere–ionosphere coupling model showing the current systems for a finite magnetopause surface wave packet: (field-aligned) magnetopause currents, ionospheric Pedersen and Hall currents, and internal currents via the complex image method. Example model outputs are shown in the low-altitude magnetosphere (b–f), ionosphere (g–l), and on the ground (m–r). Shown are the ionospheric potential (g), electric field (h), drift velocity (i), and Pedersen (j), Hall (k), and total (l) currents; along with the horizontal magnetic field (f, q) and contributions to this from magnetopause (b, m), Pedersen (c, n), Hall (d, o), and internal (e, p) currents; and finally the horizontal geoelectric field (r). Only the imaginary parts are depicted.

At equilibrium the model consists of two uniform half-spaces, the magnetosheath and magnetosphere, separated by the magnetopause around x=0. In the magnetosphere (x>0) the geomagnetic field $B_{\mathrm{0},\text{msp}}\widehat{\mathit{z}}$ is uniform with vertical field lines of length z₀, appropriate for high-latitudes. The full field line is included for self-consistency within the model, but little affects results. This geometry is translationally invariant. Thin-shell ionospheres in the northern and southern hemispheres are located at $z=-h$ and $z=-h-z_{\mathrm{0}}$, respectively. The northern hemisphere ground corresponds to z=0, where $\widehat{\mathit{x}}$ points equatorward, $\widehat{\mathit{y}}$ is westward, and $\widehat{\mathit{z}}$ is down (such coordinates are common for ULF waves, Southwood, 1974; geophysical induction, Thomson and Weaver, 1975; and ground magnetometers, Landal and Gjerloev, 2014). Therefore, in the ionosphere x<0 corresponds to open field lines in the polar cap. This is appropriate since there is little evidence that the magnetosphere is ever entirely closed under even prolonged northward interplanetary magnetic field, with in practice some regions of open flux still remaining (Raeder, 1999; Vennerstrom et al., 2005; Bhattarai and Lopez, 2013). Since scales of uniform box models cannot be globally representative and our focus is on ground-observable effects due to principally ionospheric current systems, it is at the MI-interface where the box model must be most representative of physical scales and conditions. Thus horizontal scales up to 1250 km from the origin are considered, corresponding to the distance to the horizon from ionospheric altitude. This necessarily makes our model local to the field lines which map to the mesoscale magnetopause surface wave, hence we make no comment on the global state of the magnetosphere.

Plaschke and Glassmeier (2011) derived the electrical currents for an infinite wave train of incompressible MHD magnetopause surface waves within a box model by finding solutions on closed field lines and requiring continuity of pressure and normal velocity across the (assumed infinitesimally-thin) boundary. The surface waves in this model exhibit standing structure along closed field lines between conjugate points in the ionosphere, established via ionospheric wave reflection (Kivelson and Southwood, 1988). Their expressions for the amplitudes of the surface currents must also hold for the total current within a finite-width boundary via the integral form of Ampère's law. Archer et al. (2023) demonstrated surface wave currents peak on the closed field lines at the inner/equatorward edge of the magnetopause boundary layer (MPBL), yielding longitudinally periodic FACs at the MI-interface (see Fig. 1). The equatorial magnetopause current sheet is typically ∼500 km thick (Berchem and Russell, 1982), which projected along Tsyganenko (1995) model field lines corresponds to ionospheric width d∼40 km (8 % the magnetospheric value). This is much smaller than projected surface mode wavelengths ∼250–3000 km. Therefore, for simplicity, we can represent magnetopause currents as having infinitesimal-thickness at x=0, strictly corresponding to the closed field lines at the inner/equatorward edge of the projected MPBL flux tubes. This approach does not affect ionospheric currents outside the MPBL due to Gauss's law.

$$\begin{array}{}\text{(3)}& \begin{array}{rl}{\displaystyle }\left(\begin{array}{c}{J}_{\text{mp},y}\\ J_{\text{mp},z}\end{array}\right)& {\displaystyle }={\displaystyle \frac{{\mathit{\xi }}_{\mathrm{0}}{B}_{\mathrm{0},\text{msp}}}{{\mathit{\mu }}_{\mathrm{0}}}}{\displaystyle \frac{k_z}{\sqrt{k_y^{\mathrm{2}}+k_z^{\mathrm{2}}}}}\left(\begin{array}{c}-\mathrm{i}{k}_z\sin \left(k_z\right[z+h\left]\right)\\ k_y\cos \left(k_z\right[z+h\left]\right)\end{array}\right)\\ {\displaystyle }& {\displaystyle }\times \exp \left(\mathrm{i}\right[k_{y}y-\mathit{\omega }t\left]\right)\mathit{\delta }\left(x\right),z\in [-z_{\mathrm{0}}-h,-h]\\ {\displaystyle }& {\displaystyle }\simeq J_{\mathrm{0}}\text{sgn}\left(k_y\right)\left(\begin{array}{c}-\mathrm{i}\frac{k_z}{k_y}\sin \left(k_z\right[z+h\left]\right)\\ \cos \left(k_z\right[z+h\left]\right)\end{array}\right)\\ {\displaystyle }& {\displaystyle }\times \exp \left(\mathrm{i}\right[k_{y}y-\mathit{\omega }t\left]\right)\mathit{\delta }\left(x\right),\left|k_y\right|\gg k_z\end{array}\end{array}$$

Here ξ₀ is the amplitude of boundary displacements in the equatorial magnetosphere and we have neglected the influence of the magnetosheath magnetic field on the currents since $B_{\mathrm{0},\text{msh}}\ll B_{\mathrm{0},\text{msp}}$ (Plaschke and Glassmeier, 2011). The resulting current amplitude into/out of the ionosphere $J_{\mathrm{0}}\equiv {\mathit{\xi }}_{\mathrm{0}}{B}_{\mathrm{0},\text{msp}}{k}_z/{\mathit{\mu }}_{\mathrm{0}}$ is surprisingly independent of azimuthal wavelength for all realistic ionospheric scales ($\left|k_y\right|\gg k_z$).

Since surface modes are highly-damped, we model finite trains of constant-amplitude purely-sinusoidal waves. This is achieved by multiplying the (complex phasor) currents in Eq. (3) by a rectangular function of total width Nλ centred on y=0, where N is the number of cycles, giving zero net FAC at the MI-interface. This simple rectangular modulation is unlike the typical picture of a wave packet, which usually has a smooth envelope. The finite wave packet constitutes a localised mesoscale source of FACs, as depicted in Fig. 2a. The rectangular function introduces a sinc spectrum of wavenumbers about k_y, which in general might lead to dispersion. However, the approximate surface wave dispersion relation (Eq. 1) has two dispersionless limits. The first corresponds to propagating waves where the frequency is controlled by the magnetosheath flow $\mathit{\omega }\simeq {\mathit{k}}_t\cdot {\mathit{u}}_{\mathrm{0},\text{msh}}$, typically valid in the flanks (ω₀ negligible; Chandrasekhar, 1961; Miura and Pritchett, 1982; Kozyreva et al., 2019), with phase and group velocities $\mathit{\omega }/k_y{\widehat{\mathit{k}}}_y$. The second are azimuthally stationary waves at fixed frequency ω≃ω₀, typical of the subsolar magnetopause (u_0,msh negligible; Chen and Hasegawa, 1974; Archer et al., 2021, 2023), which has zero group velocity. Stationary waves are formed from a superposition of waves with opposite azimuthal wavenumbers, equivalent to taking the imaginary part spatially while keeping the complex time variation. This makes our model static, since magnetopause currents form time-invariant structures that either propagate downtail or are stationary. Our focus is determining the spatial structure of the waves transmitted through the system. Note, the rectangular function introduces discontinuities at the wave packet edges in FACs' real parts, hence the imaginary parts (displayed throughout Fig. 2) are more representative of potential measurements.

Electrostatic MI-coupling is applied to the northern hemisphere ionosphere (e.g. Wolf, 1975), solving the ionospheric (isp) potential ψ through current continuity ($\mathrm{\nabla }\cdot \mathit{J}=\mathrm{0}$)

$$\begin{array}{}\text{(4)}& {\displaystyle }& {\displaystyle }\begin{array}{rl}{\displaystyle }{J}_{\text{mp},z}(y,z=-h)& {\displaystyle }={\mathrm{\nabla }}_{⟂}\cdot {\mathit{J}}_{\text{isp}}\\ {\displaystyle }& {\displaystyle }={\mathrm{\nabla }}_{⟂}\cdot (\mathbf{\Sigma }\cdot [-{\mathrm{\nabla }}_{⟂}\mathit{\psi }]{)}_{⟂}\\ {\displaystyle }& {\displaystyle }=-{\mathrm{\nabla }}_{⟂}\cdot \left[\left(\begin{array}{cc}{\mathrm{\Sigma }}_{\mathrm{P}}& -{\mathrm{\Sigma }}_{\mathrm{H}}\\ {\mathrm{\Sigma }}_{\mathrm{H}}& {\mathrm{\Sigma }}_{\mathrm{P}}\end{array}\right)\cdot {\mathrm{\nabla }}_{⟂}\mathit{\psi }\right]\end{array}\text{(5)}& {\displaystyle }& {\displaystyle }=-{\mathrm{\Sigma }}_{\mathrm{P}}\left({\displaystyle \frac{{\partial }^{\mathrm{2}}}{\partial x^{\mathrm{2}}}}+{\displaystyle \frac{{\partial }^{\mathrm{2}}}{\partial y^{\mathrm{2}}}}\right)\mathit{\psi },\end{array}$$

where Σ denotes the height-integrated conductivity tensor (or conductance) consisting of Pedersen (P) and Hall (H) conductances, both assumed uniform. Current continuity is appropriate in the absence of ionospheric induction when ${\mathit{\mu }}_{\mathrm{0}}{\mathrm{\Sigma }}_{\mathrm{P}}\mathit{\omega }(\mathrm{1}+{\mathrm{\Sigma }}_{\mathrm{P}}^{\mathrm{2}}/{\mathrm{\Sigma }}_{\mathrm{H}}^{\mathrm{2}})/k_{⟂}(\mathrm{1}+\coth k_{⟂}h)\ll \mathrm{1}$ (Yoshikawa, 2002; Lotko, 2004; Lysak and Song, 2006), found to be ${\mathrm{10}}^{-\mathrm{7}}\text{–}{\mathrm{10}}^{-\mathrm{5}}$ across the ranges explored. Given the ionospheric potential (example shown in Fig. 2g), the ionospheric electric field E, plasma drift velocity v, and currents are obtained

$$\begin{array}{}\text{(6)}& {\displaystyle }\mathit{E}& {\displaystyle }=-\mathrm{\nabla }\mathit{\psi }\text{(7)}& {\displaystyle }\mathit{v}& {\displaystyle }={\displaystyle \frac{\mathit{E}\times {\mathit{B}}_{\mathrm{0},\text{msp}}}{B_{\mathrm{0},\text{msp}}^{\mathrm{2}}}}\text{(8)}& {\displaystyle }{\mathit{J}}_{\mathrm{P}}& {\displaystyle }={\mathrm{\Sigma }}_{\mathrm{P}}\mathit{E}\text{(9)}& {\displaystyle }{\mathit{J}}_{\mathrm{H}}& {\displaystyle }=-{\displaystyle \frac{{\mathrm{\Sigma }}_{\mathrm{H}}}{B_{\mathrm{0},\text{msp}}}}\mathit{E}\times {\mathit{B}}_{\mathrm{0},\text{msp}}\end{array}$$

with examples displayed in Fig. 2h–l. By symmetry, signatures in the southern hemisphere should be mirror reflections of those in the northern hemisphere.

Often in modelling work the ground is treated as either a perfect insulator (e.g. Tanaka et al., 2020), where no telluric/internal currents are induced in the ground (particularly in global MHD simulations), or a perfect conductor (e.g. Hughes, 1974; Hughes and Southwood, 1974), where there is no geoelectric field and telluric currents lead to no vertical ground magnetic field while doubling the horizontal components from MI-currents. Clearly neither regime is realistic. The Complex Image Method (CIM) provides a simplification to the general induction problem (Weaver, 1971) for oscillatory MI-currents by placing image currents at a depth $z=h+\mathrm{2}\mathit{\delta }$ for complex skin depth δ (Wait and Spies, 1969; Thomson and Weaver, 1975; Boteler and Pirjola, 1998). While this approximation is technically valid for skin depths much smaller/larger than either the characteristic horizontal scale (Thomson and Weaver, 1975) or distance to the source (Wait and Spies, 1969), Pirjola and Viljanen (1998) report in practice CIM works very well across all typical parameter ranges in geomagnetic induction (much wider than considered here). They also detail how to apply CIM to arbitrary horizontal and vertical current systems. Given constant ionospheric conductance this is equivalent to image currents equal and opposite to the ionospheric Hall currents here, as shown in Fig. 2a. For simplicity we consider uniform ground conductivity σ_g, resulting in skin depth

$$\begin{array}{}\text{(10)}& \mathit{\delta }={\displaystyle \frac{\mathrm{1}}{\sqrt{\mathrm{i}\mathit{\omega }{\mathit{\mu }}_{\mathrm{0}}{\mathit{\sigma }}_{\mathrm{g}}}}}\end{array}$$

Given our model does not rely on all parameters in the surface wave dispersion relation (Eqs. 1–2) we are free to set the frequency independently.

Magnetic fields can be evaluated at any location $\mathit{r}=(x,y,z)$ using the Biot–Savart law, the general solution to Ampère's law when retarded time can be ignored (Griffiths and Heald, 1991),

$$\begin{array}{}\text{(11)}& \begin{array}{rl}{\displaystyle }& {\displaystyle }\mathit{B}\left(\mathit{r}\right)=\\ {\displaystyle }& {\displaystyle }{\displaystyle \frac{{\mathit{\mu }}_{\mathrm{0}}}{\mathrm{4}\mathit{\pi }}}\iiint {\mathrm{d}}^{\mathrm{3}}{\mathit{r}}^{\prime }{\displaystyle \frac{\mathit{J}\left({\mathit{r}}^{\prime }\right)\times (\mathit{r}-{\mathit{r}}^{\prime })}{\left[\right(x-x^{\prime }{)}^{\mathrm{2}}+(y-y^{\prime }{)}^{\mathrm{2}}+(z-z^{\prime }{)}^{\mathrm{2}}{]}^{\frac{\mathrm{3}}{\mathrm{2}}}}}\end{array}\end{array}$$

and are computed for the magnetopause, Pedersen, Hall, and internal current systems separately along with their totals, as shown in Fig. 2 on the ground (panels m–q) and above the ionosphere (panels b–f). Finally, the horizontal geoelectric field E_g is found via the Coulomb gauge magnetic vector potential A

$$\begin{array}{}\text{(12)}& {\displaystyle }& {\displaystyle }{\mathit{E}}_{\mathrm{g}}\left(\mathit{r}\right)=-{\displaystyle \frac{\partial \mathit{A}}{\partial t}}=\mathrm{i}\mathit{\omega }\mathit{A}\left(\mathit{r}\right)\text{(13)}& {\displaystyle }& {\displaystyle }\begin{array}{rl}{\displaystyle }& {\displaystyle }\mathit{A}\left(\mathit{r}\right)=\\ {\displaystyle }& {\displaystyle }{\displaystyle \frac{{\mathit{\mu }}_{\mathrm{0}}}{\mathrm{4}\mathit{\pi }}}\iiint {\mathrm{d}}^{\mathrm{3}}{\mathit{r}}^{\prime }{\displaystyle \frac{\mathit{J}\left({\mathit{r}}^{\prime }\right)}{\sqrt{(x-x^{\prime }{)}^{\mathrm{2}}+(y-y^{\prime }{)}^{\mathrm{2}}+(z-z^{\prime }{)}^{\mathrm{2}}}}}\end{array}\end{array}$$

Pirjola's (1998) procedure in our model is equivalent to applying these only to the ionospheric and image Hall currents. Unlike for the ground magnetic field, contributions from each current system have no physical meaning thus only the total is calculated, as shown in Fig. 2r. Since the skin depth is complex, calculations differ for propagating and stationary waves (the currents' imaginary part must first be taken for stationary waves).

The limitations of this deliberately simple model and how results from it may be modified through more sophisticated modelling in the future is discussed in Sect. 5.2.

2.2 Numerical implementation

Given the large number of parameters that can be varied within this model, we choose to investigate effects of the ionospheric wavelength, Pedersen and Hall conductances, ground conductivity, and the number of wave oscillations on the surface wave response throughout the system. Ionospheric conductances control the amount of current carried through the ionosphere, making them crucial to MIG-coupling. Representative values between 1.5–12 S are used (Ridley et al., 2004). Wavelengths also strongly affect Alfvén wave transmission to the ground (e.g. Hughes and Southwood, 1976). Only horizontal scales up to 1250 km from the origin are physical within our model, corresponding to when the line-of-sight vector between a ground point and an ionospheric current source, as applicable in e.g. the Biot–Savart law, is tangent to the Earth's surface (similar reasoning applied to the curvature of field lines suggest scales comparable to the radius of curvature or L-shell, both of which are significantly larger than this at the high-latitudes applicable to the magnetopause). To mitigate edge effects magnetopause/ionosphere grids extend to ±2500 km. Given this domain size, we only consider ionospheric wavelengths between 200–1600 km, corresponding to m∼8–70 or ∼1–9 R_E at the equatorial magnetopause ($\sim \mathrm{35}\times$ factor from the ionosphere), covering the majority of the observed range. Since finite wave effects are also of interest given damping, we also vary the number of oscillations from one up to the maximum possible within our domain for each wavelength. Finally, a wide range of ground conductivities applicable from city (0.1 mS m⁻¹) to salt water (5 S m⁻¹) environments are considered, given the variety of ground magnetometer locations.

All other parameters are fixed as detailed in Table 1. The field line length is representative of the dayside magnetopause (e.g. Archer et al., 2022). The magnetospheric magnetic field is that of Earth's dipole at ionospheric altitude, with 70° magnetic latitude chosen for reference (e.g. Smith and Sojka, 2019). Given the E–j paradigm, it is crucial currents (the only physical quantity coupling different regions here) at the MI-interface are correct, dictating the choice of field strength. Note, these currents are unchanged when typical equatorial magnetospheric field and scales are instead used with dipole flux tube scaling subsequently applied (Goodman, 1995), since this rescales the current by the magnetic field ratio with J₀ proportional to the local field and independent of k_y. This further justifies the use of the magnetic field at the MI-interface here. We consider a fundamental magnetopause surface wave with 1 R_E equatorial displacement, of similar magnitude to observed waves (e.g. Archer et al., 2019), which does not require rescaling to the ionosphere since its influence on the currents is accounted for by the magnetic field ratio discussed earlier. Changes to any of these fixed parameters simply rescale the wave's current amplitude J₀ proportionally, hence it is easy to adjust results presented for different scenarios. The ionospheric altitude is set at the typically used E region value of 110 km.

Table 1 Fixed parameters of the surface wave model and values used.

The ionospheric potential is solved using the Green's function (G) for the 2D Laplacian L

$$\begin{array}{ll}{\displaystyle }L& {\displaystyle }={\mathrm{\nabla }}_{\mathrm{2}D}^{\mathrm{2}}={\displaystyle \frac{{\partial }^{\mathrm{2}}}{\partial x^{\mathrm{2}}}}+{\displaystyle \frac{{\partial }^{\mathrm{2}}}{\partial y^{\mathrm{2}}}}\\ {\displaystyle }G& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathit{\pi }}}\ln R={\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathit{\pi }}}\ln \left(\sqrt{x^{\mathrm{2}}+y^{\mathrm{2}}}\right),\end{array}$$

where $R=\sqrt{x^{\mathrm{2}}+y^{\mathrm{2}}}$ denotes the perpendicular radius in cylindrical coordinates, such that

$$\begin{array}{ll}{\displaystyle }\mathit{\psi }(x,y)& {\displaystyle }=\iint \mathrm{d}{x}^{\prime }\mathrm{d}{y}^{\prime }\left[-{\displaystyle \frac{J_{mp,z}(x^{\prime },y^{\prime },z=-h)}{{\mathrm{\Sigma }}_{\mathrm{P}}}}\right]\\ {\displaystyle }& {\displaystyle }\times {\displaystyle \frac{\ln \left(\sqrt{(x-x^{\prime }{)}^{\mathrm{2}}+(y-y^{\prime }{)}^{\mathrm{2}}}\right)}{\mathrm{2}\mathit{\pi }}}\\ {\displaystyle }& {\displaystyle }=\int \mathrm{d}{y}^{\prime }\left[-{\displaystyle \frac{J_{mp,z}(y^{\prime },z=-h)}{{\mathrm{\Sigma }}_{\mathrm{P}}}}\right]\\ {\displaystyle }& {\displaystyle }\times {\displaystyle \frac{\ln \left(\sqrt{x^{\mathrm{2}}+(y-y^{\prime }{)}^{\mathrm{2}}}\right)}{\mathrm{2}\mathit{\pi }}}\end{array}$$

Integration is performed numerically using the trapezium rule, since this converges faster than Simpson's rule for periodic functions (Rahman and Schmeisser, 1990; Weideman, 2002). The ionospheric electric field is calculated with a second-order centred difference. Given the discontinuity at x=0, potentials at the midpoints between the discontinuity and adjacent gridpoints are calculated for improved gradients. The Biot–Savart law is also evaluated using the trapezium rule. These calculations are the most sensitive to grid resolution, hence we performed grid convergence tests on the ground magnetic field. We repeated calculations for 200 km wavelength (the shortest considered) in the range $\mathrm{0}\le x\le \mathrm{75}\mathrm{km}$ and fit a power law to the differences from the highest resolution (Δx=0.5 km). The average fractional error in the magnitude varied as $\mathrm{0.0021}(\mathrm{\Delta }x/\mathrm{km}{)}^{\mathrm{1.09}}$ and the angular error as $\mathrm{0.60}\mathit{°}(\mathrm{\Delta }x/\mathrm{km}{)}^{\mathrm{0.99}}$, from which we chose 1 km resolution. While the same resolution was used for the magnetopause currents in y, the large scales along the magnetospheric magnetic field render a uniform grid inefficient. Given the integrand falls off $\sim |z{|}^{-\mathrm{3}}$ at its fastest, we employ a cubically-spaced grid. Similar tests on the number of points along the field N_z showed errors converged quickly, with fractional error in magnitude of $\mathrm{244}{N}_z^{-\mathrm{2.05}}$ and angular error $\mathrm{2120}\mathit{°}{N}_z^{-\mathrm{2.00}}$, leading us to settle on N_z=200.

3 Analytic theory

Before presenting numerical results we first attempt to derive analytic relations to the surface wave effects in the ionosphere and ground, which is performed for two limiting cases. The first is the near-field, which should be well approximated by the simpler case of an infinite wave. This is particularly relevant for long wave trains, e.g. resulting from KHI (e.g. Fairfield et al., 2000) or solar wind periodic density structures (e.g. Viall et al., 2009). The second case is the far-field, pertinent to localised sources of FACs which result from a limited number of cycles, e.g. due to damping (Chen and Hasegawa, 1974; Kozyreva et al., 2019).

3.1 Ionosphere

3.1.1 Near-field approximation

For an infinite wave, Eq. (5) for the ionospheric potential becomes

$$\begin{array}{rl}{\displaystyle }& {\displaystyle }\left({\displaystyle \frac{{\partial }^{\mathrm{2}}}{\partial x^{\mathrm{2}}}}+{\displaystyle \frac{{\partial }^{\mathrm{2}}}{\partial y^{\mathrm{2}}}}\right)\mathit{\psi }(x,y,t)\\ {\displaystyle }& {\displaystyle }=-{\displaystyle \frac{J_{\mathrm{0}}}{{\mathrm{\Sigma }}_{\mathrm{P}}}}\text{sgn}\left(k_y\right)\exp \left(\mathrm{i}\right[k_{y}y-\mathit{\omega }t\left]\right)\mathit{\delta }\left(x\right)\end{array}$$

By separation of variables $\mathit{\psi }(x,y,t)=\mathit{\psi }\left(x\right)\exp \left(\mathrm{i}\right[k_{y}y-\mathit{\omega }t\left]\right)$ we have

$$\left({\displaystyle \frac{{\mathrm{d}}^{\mathrm{2}}}{\mathrm{d}{x}^{\mathrm{2}}}}-k_y^{\mathrm{2}}\right)\mathit{\psi }\left(x\right)=-{\displaystyle \frac{J_{\mathrm{0}}}{{\mathrm{\Sigma }}_{\mathrm{P}}}}\text{sgn}\left(k_y\right)\mathit{\delta }\left(x\right)$$

This is the Green's function problem for the “screened Poisson” or “modified Helmholtz” operator which leads to solutions (Abramowitz and Stegun, 2000; McLean, 2000; Schulz, 2001)

$$\begin{array}{ll}{\displaystyle }\mathit{\psi }\left(x\right)& {\displaystyle }={\displaystyle \frac{J_{\mathrm{0}}}{\mathrm{2}{\mathrm{\Sigma }}_{\mathrm{P}}}}{\displaystyle \frac{\mathrm{1}}{k_y}}\exp (-|k_{y}x\left|\right)\\ {\displaystyle }\mathit{\psi }(x,y,t)& {\displaystyle }={\displaystyle \frac{J_{\mathrm{0}}}{\mathrm{2}{\mathrm{\Sigma }}_{\mathrm{P}}{k}_y}}\exp (-|k_{y}x|+\mathrm{i}[k_{y}y-\mathit{\omega }t\left]\right)\end{array}$$

from which it is straightforward to derive the ionospheric electric field, drift velocity, and Pedersen and Hall currents

$$\begin{array}{ll}{\displaystyle }\mathit{E}=& {\displaystyle }\left(\begin{array}{c}{E}_x\\ E_y\end{array}\right)=-{\displaystyle \frac{J_{\mathrm{0}}}{\mathrm{2}{\mathrm{\Sigma }}_{\mathrm{P}}}}\left(\begin{array}{c}\text{sgn}\left(k_y\right)\\ -\mathrm{i}\end{array}\right)\exp (-|k_{y}x|+\mathrm{i}[k_{y}y-\mathit{\omega }t\left]\right)\\ {\displaystyle }\mathit{v}=& {\displaystyle }\left(\begin{array}{c}{v}_x\\ v_y\end{array}\right)=-{\displaystyle \frac{J_{\mathrm{0}}}{\mathrm{2}{\mathrm{\Sigma }}_{\mathrm{P}}{B}_{\mathrm{0},\text{msp}}}}\left(\begin{array}{c}\mathrm{i}\\ \text{sgn}\left(k_y\right)\end{array}\right)\\ {\displaystyle }& {\displaystyle }\times \exp (-|k_{y}x|+\mathrm{i}[k_{y}y-\mathit{\omega }t\left]\right)\end{array}$$

$$\begin{array}{}\text{(14)}& {\displaystyle }& {\displaystyle }\begin{array}{rl}{\displaystyle }{\mathit{J}}_{\mathrm{P}}=& {\displaystyle }\left(\begin{array}{c}{J}_{\mathrm{P},x}\\ J_{\mathrm{P},y}\end{array}\right)={\displaystyle \frac{J_{\mathrm{0}}}{\mathrm{2}}}\left(\begin{array}{c}\text{sgn}\left(k_{y}x\right)\\ -\mathrm{i}\end{array}\right)\\ {\displaystyle }& {\displaystyle }\times \exp (-|k_{y}x|+\mathrm{i}[k_{y}y-\mathit{\omega }t\left]\right)\end{array}\text{(15)}& {\displaystyle }& {\displaystyle }\begin{array}{rl}{\displaystyle }{\mathit{J}}_{\mathrm{H}}=& {\displaystyle }\left(\begin{array}{c}{J}_{\mathrm{H},x}\\ J_{\mathrm{H},y}\end{array}\right)={\displaystyle \frac{J_{\mathrm{0}}}{\mathrm{2}}}{\displaystyle \frac{{\mathrm{\Sigma }}_{\mathrm{H}}}{{\mathrm{\Sigma }}_{\mathrm{P}}}}\left(\begin{array}{c}\mathrm{i}\\ \text{sgn}\left(k_{y}x\right)\end{array}\right)\\ {\displaystyle }& {\displaystyle }\times \exp (-|k_{y}x|+\mathrm{i}[k_{y}y-\mathit{\omega }t\left]\right)\end{array}\end{array}$$

3.1.2 Far-field approximation

A finite FAC distribution centred on the origin with zero net current into/out of the ionosphere is a localised source of ionospheric dynamics. We can employ a multipole expansion of the electrostatic potential far from the extent of the current sources ($R\gg N\mathit{\lambda }/\mathrm{2}$) to understand the large-scale impact of this localised source

$$\begin{array}{ll}{\displaystyle }\mathit{\psi }(x,y)& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{{\mathrm{\Sigma }}_{\mathrm{P}}}}\mathit{p}\cdot \mathrm{\nabla }G\\ {\displaystyle }& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{{\mathrm{\Sigma }}_{\mathrm{P}}}}\mathit{p}\cdot \mathrm{\nabla }{\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathit{\pi }}}\ln R\\ {\displaystyle }& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathit{\pi }{\mathrm{\Sigma }}_{\mathrm{P}}}}{\displaystyle \frac{\mathit{p}\cdot \mathit{r}}{R^{\mathrm{2}}}},\end{array}$$

where the current dipole moment p (analogous to the electric dipole moment in electrostatics) is

$$\begin{array}{ll}{\displaystyle }\mathit{p}& {\displaystyle }=-\int {\mathit{r}}^{\prime }{J}_z\left({\mathit{r}}^{\prime }\right){\mathrm{d}}^{\mathrm{2}}{\mathit{r}}^{\prime }\\ {\displaystyle }& {\displaystyle }=J_{\mathrm{0}}\text{sgn}\left(k_y\right)\underset{-N\mathit{\pi }/k_y}{\overset{N\mathit{\pi }/k_y}{\int }}{\mathit{y}}^{\prime }\exp \left(\mathrm{i}{k}_y{y}^{\prime }\right)\mathrm{d}{y}^{\prime }\\ {\displaystyle }& {\displaystyle }=J_{\mathrm{0}}\text{sgn}\left(k_y\right){\displaystyle \frac{\mathrm{2}\mathrm{i}N\mathit{\pi }(-\mathrm{1}{)}^{N-\mathrm{1}}}{k_y^{\mathrm{2}}}}\widehat{\mathit{y}}\end{array}$$

The concept of a current dipole is more readily used in biophysics (e.g. Cohen and Hosaka, 1976; Sarvas, 1987) and geophysics (Poikonen et al., 1997; GeoSci.xyz Project, 2015), describing a localised current source within a conducting medium. Indeed, our finite wave packet can be thought of as a string of current dipole sources (see Fig. 2a), hence the scaling with N. This leads to ionospheric potential in polar coordinates

$$\mathit{\psi }={\displaystyle \frac{p_y}{\mathrm{2}\mathit{\pi }{\mathrm{\Sigma }}_{\mathrm{P}}}}{\displaystyle \frac{\sin \mathit{\theta }}{R}}=\mathrm{i}{J}_{\mathrm{0}}\text{sgn}\left(k_y\right){\displaystyle \frac{N(-\mathrm{1}{)}^{N-\mathrm{1}}}{{\mathrm{\Sigma }}_{\mathrm{P}}{k}_y^{\mathrm{2}}}}{\displaystyle \frac{\sin \mathit{\theta }}{R}}$$

from which the ionospheric electric field, drift velocity, and currents are derived (Eqs. 6–9)

$$\begin{array}{ll}{\displaystyle }\mathit{E}& {\displaystyle }={\displaystyle \frac{p_y}{\mathrm{2}\mathit{\pi }{\mathrm{\Sigma }}_{\mathrm{P}}}}\left(-{\displaystyle \frac{\sin \mathit{\theta }}{R^{\mathrm{2}}}}\widehat{\mathit{R}}+{\displaystyle \frac{\cos \mathit{\theta }}{R^{\mathrm{2}}}}\widehat{\mathit{\theta }}\right)\\ {\displaystyle }& {\displaystyle }=\mathrm{i}{J}_{\mathrm{0}}\text{sgn}\left(k_y\right){\displaystyle \frac{N(-\mathrm{1}{)}^{N-\mathrm{1}}}{{\mathrm{\Sigma }}_{\mathrm{P}}{k}_y^{\mathrm{2}}}}{\displaystyle \frac{\mathrm{1}}{R^{\mathrm{2}}}}\left(\begin{array}{c}-\sin \mathrm{2}\mathit{\theta }\\ \cos \mathrm{2}\mathit{\theta }\end{array}\right)\\ {\displaystyle }\mathit{v}& {\displaystyle }={\displaystyle \frac{p_y}{\mathrm{2}\mathit{\pi }{\mathrm{\Sigma }}_{\mathrm{P}}{B}_{\mathrm{0},\text{msp}}}}\left({\displaystyle \frac{\cos \mathit{\theta }}{R^{\mathrm{2}}}}\widehat{\mathit{R}}+{\displaystyle \frac{\sin \mathit{\theta }}{R^{\mathrm{2}}}}\widehat{\mathit{\theta }}\right)\\ {\displaystyle }& {\displaystyle }=\mathrm{i}\text{sgn}\left(k_y\right){\displaystyle \frac{J_{\mathrm{0}}N(-\mathrm{1}{)}^{N-\mathrm{1}}}{{\mathrm{\Sigma }}_{\mathrm{P}}{B}_{\mathrm{0},\text{msp}}{k}_y^{\mathrm{2}}}}{\displaystyle \frac{\mathrm{1}}{R^{\mathrm{2}}}}\left(\begin{array}{c}\cos \mathrm{2}\mathit{\theta }\\ \sin \mathrm{2}\mathit{\theta }\end{array}\right)\\ {\displaystyle }{\mathit{J}}_{\mathrm{P}}& {\displaystyle }={\displaystyle \frac{p_y}{\mathrm{2}\mathit{\pi }}}\left(-{\displaystyle \frac{\sin \mathit{\theta }}{R^{\mathrm{2}}}}\widehat{\mathit{R}}+{\displaystyle \frac{\cos \mathit{\theta }}{R^{\mathrm{2}}}}\widehat{\mathit{\theta }}\right)\\ {\displaystyle }& {\displaystyle }=\mathrm{i}\text{sgn}\left(k_y\right){\displaystyle \frac{J_{\mathrm{0}}N(-\mathrm{1}{)}^{N-\mathrm{1}}}{k_y^{\mathrm{2}}}}{\displaystyle \frac{\mathrm{1}}{R^{\mathrm{2}}}}\left(\begin{array}{c}-\sin \mathrm{2}\mathit{\theta }\\ \cos \mathrm{2}\mathit{\theta }\end{array}\right)\\ {\displaystyle }{\mathit{J}}_{\mathrm{H}}& {\displaystyle }=-{\displaystyle \frac{p_y}{\mathrm{2}\mathit{\pi }}}{\displaystyle \frac{{\mathrm{\Sigma }}_{\mathrm{H}}}{{\mathrm{\Sigma }}_{\mathrm{P}}}}\left({\displaystyle \frac{\cos \mathit{\theta }}{R^{\mathrm{2}}}}\widehat{\mathit{R}}+{\displaystyle \frac{\sin \mathit{\theta }}{R^{\mathrm{2}}}}\widehat{\mathit{\theta }}\right)\\ {\displaystyle }& {\displaystyle }=\mathrm{i}\text{sgn}\left(k_y\right){\displaystyle \frac{J_{\mathrm{0}}N(-\mathrm{1}{)}^N}{k_y^{\mathrm{2}}}}{\displaystyle \frac{{\mathrm{\Sigma }}_{\mathrm{H}}}{{\mathrm{\Sigma }}_{\mathrm{P}}}}{\displaystyle \frac{\mathrm{1}}{R^{\mathrm{2}}}}\left(\begin{array}{c}\cos \mathrm{2}\mathit{\theta }\\ \sin \mathrm{2}\mathit{\theta }\end{array}\right)\end{array}$$

Note these solutions are static in the frame of the wave packet travelling at $\mathit{\omega }/k_y$ for propagating waves, hence are not oscillatory.

3.2 Ground magnetic field

The Biot–Savart law (Eq. 11) applied to ionospheric currents for an infinite wave involves complicated integrals over sources in both x^′ and y^′ (Eqs. 14–15). While it is possible to evaluate the y^′ integral, recovering that the ground magnetic field is periodic in y with wavenumber k_y as expected, this leaves integrals of the form

$$\int \mathrm{d}{x}^{\prime }\exp (-|k_y{x}^{\prime }\left|\right){\displaystyle \frac{K_{\mathit{\alpha }}\left(k_y\sqrt{(x-x^{\prime }{)}^{\mathrm{2}}+h^{\mathrm{2}}}\right)}{\left(\right(x-x^{\prime }{)}^{\mathrm{2}}+h^{\mathrm{2}}{)}^{\frac{\mathit{\beta }}{\mathrm{2}}}}}$$

where K is the modified Bessel functions of the second kind and different orders/exponents (α,β) apply to each component. We are not aware of any closed form for these integrals, hence cannot provide analytic solutions for the ground magnetic field in the near-field limit.

In contrast, the far-field approximation can be derived by considering an ideal current dipole – a pair of opposite semi-infinite FAC sources separated by distance d along the y axis. If we assume Fukushima's (1976) theorem is approximately valid, given that in our model field lines are vertical and ionospheric conductances are uniform, then one needs only consider the contributions from Hall currents. Using Eq. (13), Pirjola and Viljanen (1998) gave solutions to the magnetic vector potential due to Pedersen currents from a single FAC I at the origin

$${\mathit{A}}_{\mathrm{P}}=-{\displaystyle \frac{{\mathit{\mu }}_{\mathrm{0}}I}{\mathrm{4}\mathit{\pi }}}{\displaystyle \frac{R}{\sqrt{R^{\mathrm{2}}+(z+h{)}^{\mathrm{2}}}+z+h}}\widehat{\mathit{R}},z>-h$$

Since in our model ${\mathit{J}}_{\mathrm{H}}={\mathrm{\Sigma }}_{\mathrm{H}}/{\mathrm{\Sigma }}_{\mathrm{P}}{\mathit{J}}_{\mathrm{P}}\times \widehat{\mathit{z}}$, the vector potential due to Hall currents is

$${\mathit{A}}_{\mathrm{H}}={\displaystyle \frac{{\mathit{\mu }}_{\mathrm{0}}I}{\mathrm{4}\mathit{\pi }}}{\displaystyle \frac{{\mathrm{\Sigma }}_{\mathrm{H}}}{{\mathrm{\Sigma }}_{\mathrm{P}}}}{\displaystyle \frac{R}{\sqrt{R^{\mathrm{2}}+(z+h{)}^{\mathrm{2}}}+z+h}}\widehat{\mathit{\theta }},z>-h$$

with ground (z=0) magnetic field

$$\begin{array}{}\text{(16)}& \begin{array}{rl}{\displaystyle }& {\displaystyle }{\mathit{B}}_{\mathrm{H}}=\mathrm{\nabla }\times {\mathit{A}}_{\mathrm{H}}={\displaystyle \frac{{\mathit{\mu }}_{\mathrm{0}}I}{\mathrm{4}\mathit{\pi }}}{\displaystyle \frac{{\mathrm{\Sigma }}_{\mathrm{H}}}{{\mathrm{\Sigma }}_{\mathrm{P}}}}\\ {\displaystyle }& {\displaystyle }\times \left[{\displaystyle \frac{R}{\sqrt{R^{\mathrm{2}}+h^{\mathrm{2}}}\left(\sqrt{R^{\mathrm{2}}+h^{\mathrm{2}}}+h\right)}}\widehat{\mathit{R}}+{\displaystyle \frac{\mathrm{1}}{\sqrt{R^{\mathrm{2}}+h^{\mathrm{2}}}}}\widehat{\mathit{z}}\right]\end{array}\end{array}$$

Considering sources at $y=\pm d/\mathrm{2}$ and series expanding about d=0 yields

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{B}_{\mathrm{H},x}={\displaystyle \frac{{\mathit{\mu }}_{\mathrm{0}}I}{\mathrm{4}\mathit{\pi }}}{\displaystyle \frac{{\mathrm{\Sigma }}_{\mathrm{H}}}{{\mathrm{\Sigma }}_{\mathrm{P}}}}[{\displaystyle \frac{x}{R^{\mathrm{2}}+h^{\mathrm{2}}+h\sqrt{R^{\mathrm{2}}+h^{\mathrm{2}}}}}\\ {\displaystyle }& {\displaystyle }\pm {\displaystyle \frac{d}{\mathrm{2}}}{\displaystyle \frac{xy(\mathrm{2}\sqrt{R^{\mathrm{2}}+h^{\mathrm{2}}}+h)}{(R^{\mathrm{2}}+h^{\mathrm{2}}{)}^{\frac{\mathrm{3}}{\mathrm{2}}}(\sqrt{R^{\mathrm{2}}+h^{\mathrm{2}}}+h{)}^{\mathrm{2}}}}+\mathcal{O}\left(d^{\mathrm{2}}\right)]\\ {\displaystyle }& {\displaystyle }{B}_{\mathrm{H},y}={\displaystyle \frac{{\mathit{\mu }}_{\mathrm{0}}I}{\mathrm{4}\mathit{\pi }}}{\displaystyle \frac{{\mathrm{\Sigma }}_{\mathrm{H}}}{{\mathrm{\Sigma }}_{\mathrm{P}}}}[{\displaystyle \frac{y}{R^{\mathrm{2}}+h^{\mathrm{2}}+h\sqrt{R^{\mathrm{2}}+h^{\mathrm{2}}}}}\\ {\displaystyle }& {\displaystyle }\pm {\displaystyle \frac{d}{\mathrm{2}}}{\displaystyle \frac{y^{\mathrm{2}}(\mathrm{2}\sqrt{R^{\mathrm{2}}+h^{\mathrm{2}}}+h)}{(R^{\mathrm{2}}+h^{\mathrm{2}}{)}^{\frac{\mathrm{3}}{\mathrm{2}}}(\sqrt{R^{\mathrm{2}}+h^{\mathrm{2}}}+h{)}^{\mathrm{2}}}}\\ {\displaystyle }& {\displaystyle }\mp {\displaystyle \frac{d}{\mathrm{2}}}{\displaystyle \frac{\mathrm{1}}{R^{\mathrm{2}}+h^{\mathrm{2}}+h\sqrt{R^{\mathrm{2}}+h^{\mathrm{2}}}}}+\mathcal{O}\left(d^{\mathrm{2}}\right)]\\ {\displaystyle }& {\displaystyle }{B}_{\mathrm{H},z}={\displaystyle \frac{{\mathit{\mu }}_{\mathrm{0}}I}{\mathrm{4}\mathit{\pi }}}{\displaystyle \frac{{\mathrm{\Sigma }}_{\mathrm{H}}}{{\mathrm{\Sigma }}_{\mathrm{P}}}}[{\displaystyle \frac{\mathrm{1}}{\sqrt{R^{\mathrm{2}}+h^{\mathrm{2}}}}}\pm {\displaystyle \frac{d}{\mathrm{2}}}{\displaystyle \frac{y}{(R^{\mathrm{2}}+h^{\mathrm{2}}{)}^{\frac{\mathrm{3}}{\mathrm{2}}}}}+\mathcal{O}(d^{\mathrm{2}}\left)\right]\end{array}$$

Summing over both sources, converting back to cylindrical coordinates, and remembering the current dipole moment p_y=Id here, we then have

$$\begin{array}{}\text{(17)}& {\displaystyle }& {\displaystyle }\begin{array}{rl}{\displaystyle }& {\displaystyle }{B}_{\mathrm{H},\text{hor}}={\displaystyle \frac{{\mathit{\mu }}_{\mathrm{0}}{p}_y}{\mathrm{4}\mathit{\pi }}}{\displaystyle \frac{{\mathrm{\Sigma }}_{\mathrm{H}}}{{\mathrm{\Sigma }}_{\mathrm{P}}}}{\displaystyle \frac{\mathrm{1}}{R^{\mathrm{2}}+h^{\mathrm{2}}+h\sqrt{R^{\mathrm{2}}+h^{\mathrm{2}}}}}\\ {\displaystyle }& {\displaystyle }\times \left\{\sin \mathit{\theta }\left[{\displaystyle \frac{R^{\mathrm{2}}\left(\mathrm{2}\mathit{\{}{R}^{\mathrm{2}}+h^{\mathrm{2}}\mathit{\}}+h\sqrt{R^{\mathrm{2}}+h^{\mathrm{2}}}\right)}{\left(R^{\mathrm{2}}+h^{\mathrm{2}}+h\sqrt{R^{\mathrm{2}}+h^{\mathrm{2}}}\right)(R^{\mathrm{2}}+h^{\mathrm{2}})}}-\mathrm{1}\right]\right.\\ {\displaystyle }& {\displaystyle }\left.\phantom{\left[{\displaystyle \frac{R^{\mathrm{2}}\left(\mathrm{2}\mathit{\{}{R}^{\mathrm{2}}+h^{\mathrm{2}}\mathit{\}}+h\sqrt{R^{\mathrm{2}}+h^{\mathrm{2}}}\right)}{\left(R^{\mathrm{2}}+h^{\mathrm{2}}+h\sqrt{R^{\mathrm{2}}+h^{\mathrm{2}}}\right)(R^{\mathrm{2}}+h^{\mathrm{2}})}}-\mathrm{1}\right]}\times \widehat{\mathit{R}}-\cos \mathit{\theta }\widehat{\mathit{\theta }}\right\}\end{array}\text{(18)}& {\displaystyle }& {\displaystyle }\simeq {\displaystyle \frac{{\mathit{\mu }}_{\mathrm{0}}{p}_y}{\mathrm{4}\mathit{\pi }}}{\displaystyle \frac{{\mathrm{\Sigma }}_{\mathrm{H}}}{{\mathrm{\Sigma }}_{\mathrm{P}}}}{\displaystyle \frac{\mathrm{1}}{R^{\mathrm{2}}}}\mathit{\{}\sin \mathit{\theta }\widehat{\mathit{R}}-\cos \mathit{\theta }\widehat{\mathit{\theta }}\mathit{\}},R\gg h\text{(19)}& {\displaystyle }& {\displaystyle }{B}_{\mathrm{H},z}={\displaystyle \frac{{\mathit{\mu }}_{\mathrm{0}}{p}_y}{\mathrm{4}\mathit{\pi }}}{\displaystyle \frac{{\mathrm{\Sigma }}_{\mathrm{H}}}{{\mathrm{\Sigma }}_{\mathrm{P}}}}{\displaystyle \frac{R\sin \mathit{\theta }}{(R^{\mathrm{2}}+h^{\mathrm{2}}{)}^{\frac{\mathrm{3}}{\mathrm{2}}}}}\text{(20)}& {\displaystyle }& {\displaystyle }\simeq {\displaystyle \frac{{\mathit{\mu }}_{\mathrm{0}}{p}_y}{\mathrm{4}\mathit{\pi }}}{\displaystyle \frac{{\mathrm{\Sigma }}_{\mathrm{H}}}{{\mathrm{\Sigma }}_{\mathrm{P}}}}{\displaystyle \frac{\sin \mathit{\theta }}{R^{\mathrm{2}}}},R\gg h\end{array}$$

The limit R≫h is valid for typical magnetopause wavelengths when considering the far-field.

The Biot–Savart law applied to internal (image) currents takes the same form as for MI-currents, thus it is not possible to solve analytically in the near-field limit. In the far-field limit, disturbances due to a propagating wave packet are not periodic, instead varying slowly compared to ULF timescales leading to negligible telluric currents. On the other hand, stationary waves' disturbances will be periodic and one could use a similar approach to that above to derive approximate expressions. Given the more limited scope and that fields in this limit are likely small, we shall not consider them further.

3.3 Geoelectric field

As with the ground magnetic field, near-field solutions to the geoelectric field have no closed form we are aware of. Far-field solutions are again only appropriate for stationary waves, which for similar reasons we shall not assess here.

3.4 Outlook

It is clear that, even for the simplest model setup under several approximations, analytic solutions of surface wave transmission throughout the MIG-system are highly limited. The near-field for even an infinite wave has no closed form on the ground. While far-field solutions for finite FAC sources can be derived throughout the system, finite wave effects in the near-field (where they are most important) cannot be expressed in terms of elementary functions. This motivates numerical results within our model, guided by analytic theory, to elucidate scaling relations for surface wave transmission across the input parameters of interest.

4 Numerical results

We now present numerical results from our model, comparing these to theory from Sect. 3 where possible.

4.1 Ionosphere

The ionospheric impacts show distinct behaviours in the near- and far-field limits, as seen in Fig. 2g–l (imaginary parts are shown). Near the wave packet centre the potential (panel g) is similarly periodic in y (East–West) to the FACs, peaking at their maxima/minima (white circles), and decaying with distance perpendicular to the MPBL-projection, consistent with Sect. 3.1.1. These result in periodic ionospheric electric fields, velocities, and currents that all exibit maximum magnitudes along the MPBL. Pedersen currents (panel j) emanating from FAC sources close at their adjacent opposite polarity sinks, but are not simply confined to the East–West direction spreading out laterally also. An identical pattern occurs in the electric field (panel h) since the conductance is uniform. Hall currents form vortices around FAC maxima/minima (panel k), which are clockwise for downwards FACs and anticlockwise for upwards. Ionospheric flow vortices (panel i) with identical patterns but the opposite rotation sense are a result of the uniform conductances and magnetic field. All these features are qualitatively consistent with MHD simulations (Archer et al., 2023). Far-field ionospheric signatures are instead strongly influenced by the potentials near the ends of the wave train, which exhibit larger maxima/minima than in the centre. These result in electric fields and Pedersen currents that resemble a dipole field, as predicted in Sect. 3.1.2, and two large-scale vortices in both the flow and Hall currents. While this was not clear in the simulation of Archer et al. (2023), likely due to the larger wavelengths making spherical effects important, flow patterns at mid-latitudes certainly exhibited scales larger than the distance between adjacent opposite-polarity FACs.

Figure 3 Example ionospheric cuts along x=0 (a–f) and y=0 (g–l) of the field-aligned current (a, g), potential (b, h), electric field (c, i), drift velocity (d, j) and Pedersen (e, k) and Hall (f, l) currents. Both imaginary (thick, darker) and real (thin, lighter) parts of the solutions are shown. Background trends are shown as dashed lines where the Spearman rank correlation >0.5.

Figure 3 shows cuts along and perpendicular to the MPBL for the same model run, depicting both real (thinner/lighter) and imaginary (thicker/darker) parts. It is clear the surface wave's sinusoidal FACs at the MI-interface (panel a) result in similarly sinusoidal ionospheric properties within the extent of the wave train (panels b–f). However, fluctuations are not about zero as slight offsets and/or trends occur throughout the surface wave, most clearly seen in the potential (panel b) where the real part has a positive offset and the imaginary part exhibits a slope. Mathematically these arise from the Green's function integration resulting in a modulated complex phasor superposed with the logarithm of a polynomial. We apply a spline-based filtering method detailed in appendix A to extract the offsets/trends, shown by the dashed lines. Periodic fluctuations are then defined as the difference between the two, which we denote with tilde marks. These offsets/trends occur within the surface wave due to its finite extent. For an infinite wave an equally infinite number of FAC oscillations about zero exist to the left and right of any point along the MPBL, which leads to overall cancellation in the Green's function integration. In contrast, for a finite wave imbalances occur due to the applied rectangular function leading to unequal net currents to the left and right of points along the MPBL. While the discontinuity in the FACs' real part from the rectangular function leads to spikes in the real parts of all ionospheric quantities except the potential, only the imaginary parts are treated as physical. Along the MPBL, all signatures slowly reduce in magnitude with distance outside the wave packet. In contrast, perpendicular to the MPBL (panels h–l), ionospheric quantities appear to quickly exponentially decay with distance (as expected from near-field theory) with perturbations becoming negligible well within a single wavelength's distance. This means even relatively large-scale surface waves have highly localised impacts latitudinally. The model suggests ionospheric velocities up to several km s⁻¹ and Joule heating rates up to hundreds of mW m⁻² might result from large-amplitude magnetopause surface wave activity, those these peak values will be spatially concentrated. Ionospheric flows of tens of km s⁻¹ have been reported in global simulations due to large-scale magnetopause motion (Slinker et al., 1999).

Figure 4 Comparison of the ionospheric potential with near-field theory. (a, b) Example comparison across a single run. The fluctuating part of the potential is divided by a complex phasor, with real (a) and imaginary (b) parts shown as a function of x for various y-values (colours). (c–h) Comparison across runs varying model input parameters. The real part of the fluctuating potential along y=0 is shown (c–e) along with a normalised version of this (f–h) according to theory, such that curves should collapse on to one another and the theoretical solution (black dashed lines).

We now quantitatively compare ionospheric surface wave signatures to theory. Figure 4 examines the potential in the near-field limit. For a single set of model inputs, we extract the fluctuations in the potential and divide these by the expected complex phasor variation, with panels (a) and (b) depicting the real and imaginary parts, respectively, of the result as a function of x at different points along the y axis. This shows excellent agreement with the theoretical exponentially-decaying solution, given by the black dashed lines. Note, the e-folding scales for mesoscale waves (∼30–250 km) will not be well-resolved by typical ionospheric grid resolutions in global simulations ($\mathrm{1}\mathit{°}/\mathrm{110}\mathrm{km}$; Rastätter et al., 2014). For the dominant real part, the amplitude at the MPBL only significantly differs from theory when within a quarter-wavelength of the wave packet edge, where the infinite wave assumption breaks down. There is also a slight increase in discrepancies with distance from the MPBL. Close to the MPBL the potential will be dominated by the local FAC oscillations, whereas further away the influence from other adjacent oscillations will proportionally increase leading to differences between finite and infinite waves. Differences in the imaginary part tend to peak at/near the MPBL and grow with distance from the wave packet centre. We now also investigate how the potential changes when varying the model inputs. Panels (c)–(e) show the real part of the fluctuating potential along y=0, whereas panels (f)–(h) normalise this based on theory such that curves should collapse on to one another and the theoretical solution (black dashed). First we vary the number oscillations (panels c and f) from one up to the maximum possible within our model domain for this wavelength (later denoted N_max), keeping all other inputs fixed. The number of oscillations does not significantly affect the amplitude or shape of the potential, with differences between all solutions and the theoretical curve only apparent when N is two or less. Increasing the wavelength (panels d and g) results in the fluctuating potential having increased amplitude and larger lateral extent, though these remain very localised with Half Width at Half Maximum (HWHM) values ranging between only ∼20–175 km for the wavelength range considered. The normalisation, which rescales both distances and amplitudes by k_y, leads to excellent agreement between the different runs and the theoretical solution, with only slight differences as k_yx becomes large. Finally, the Pedersen conductance (panels e and h) affects only the amplitude of the fluctuating potential with normalisation giving agreement between runs to within machine accuracy, as the constant conductance is simply a factor. Overall, the ionospheric signatures of a finite magnetopause surface wave are well approximated by the infinite wave solutions for small $k_{y}R/N$.

Figure 5 Comparison of the ionospheric potential with near-field theory in a similar format to Fig. 4. (a, b) Example comparison across a single run. The real (a) and imaginary (b) parts of the potential are divided by sin θ and shown as a function of perpendicular radial distance for various azimuths. (c–h) Comparison across runs varying model input parameters. The imaginary part of the potential along x=0 is shown (c–e) along with its theory-normalisation version (f–h). Theoretical solutions (blacked dashed lines) are also shown. Vertical dotted lines indicate the radial extent of the current distribution.

For the far-field, we initially compare the ionospheric potential to theory for a single set of model inputs. Figure 5a–b show the real and imaginary parts, respectively, of the potential divided by sin θ as a function of perpendicular radial distance for a range of polar angles. While the red curves are along the periodic surface wave FACs, resulting in undulating potentials within the finite wave packet's extent (denoted by the vertical dotted line), the yellow curves are perpendicular to the MPBL hence do not exhibit this feature. For distances greater than the wave packet size, applicable to the far-field, the potential's real part tends towards zero and the imaginary part asymptotes to the theoretical $\mathrm{1}/R$ dependence given by the black dashed line. The numerical solutions appear to converge faster as the polar angle approaches 90°, i.e. along the MPBL. As before, we then vary model inputs individually, comparing the potential's imaginary part along the MPBL. The sign of the far-field potential (panel c) changes with N, with odd N yielding positive values and even N giving negative ones. The normalisation used here rescales distances according to N, such that the wave train sizes become identical, given by the vertical dotted line in panels (f)–(h). This results in the different model solutions tending towards one another and far-field theory outside of the wave packet. It is clear that, similar to the near-field case, increasing the wavelength (panel d) leads to both larger scales and amplitudes. Accounting for these in the normalisation (panel g) leads to excellent agreement between the different potentials, which again asymptote towards theory. Finally, the Pedersen conductance (panels e and h) again serves as a simple dividing factor. We thus conclude that far-field theory agrees well with our numerical solutions for large $k_{y}R/N$.

4.2 Ground

Model outputs for the ground, as shown in Fig. 2m–r, are the contributions to the ground magnetic field from magnetospheric, ionospheric, and internal currents, as well as the geoelectric field.

4.2.1 Magnetic field (MI-currents)

We first consider the ground magnetic field only due to magnetospheric and ionospheric currents. The horizontal field patterns associated with the magnetopause and Pedersen currents appear similar but opposite to one another, meaning the ground field is dominated by Hall currents. Such cancellation is only perfect when field lines are vertical, ionospheric conductances are uniform, and FACs are of semi-infinite extent (Fukushima, 1976), whereas only the first two criteria strictly hold here. As with the ionosphere there are distinct near- and far-field patterns. The far-field, which resembles two vortices for the magnetopause and Pedersen current systems and a dipole for the Hall currents, dominate the ground magnetic field compared to near-field perturbations, which are difficult to discern in the streamlines but cause clear modulation of the horizontal field strength. This is likely the ionosphere screening the smaller scales (Hughes and Southwood, 1976). Further evidence of this can be seen comparing the magnetic field computed at the same distance above the ionosphere (panels b–f), which exhibit overall larger field strengths principally due to magnetopause currents (ionospheric currents yield equal but opposite contributions to those on the ground). Finally, ground field orientation appears close to a systematic 90° rotation from that above the ionosphere due to magnetospheric currents. While expected for an infinite plane Alfvén wave when the background field is vertical and conductances are uniform (Hughes, 1974; Hughes and Southwood, 1974), Archer et al. (2023) showed this does not exactly hold for a surface wave. Indeed, closer inspection of the rotation angle yields a bimodal distribution with distinct closely located peaks at 85° and 95° and overall standard deviation 14° for this example. At distinct wave phases along the MPBL-projection the sense of the rotation even becomes reversed, i.e. −90°, returning to normal within half a wavelength perpendicular to the MPBL. This further demonstrates the rule-of-thumb becomes altered for even simple changes to wave structure, such as with a finite surface wave.

Figure 6 Example ground cuts along y=0 (a–e) and x=0 (h–l) in a similar format to Fig. 3. Contributions to the ground magnetic field from magnetopause (a, h), Pedersen (b, i), Hall (c, j) and internal (d–e, k–l) currents, and the geoelectric field (f–g, m–n) are shown. Where required these are shown separately for propagating and stationary surface waves.

Figure 6 shows cuts along (panels a–c) and perpendicular to (h–j) the projected MPBL. Similar to with the ionosphere, the ground magnetic field shows sinusoidal behaviour along y about background trends/offsets (estimated by the dashed lines). These trends/offsets appear more significant than in the ionosphere, e.g. the fluctuation amplitude only constitutes half the total field perturbation near the MPBL. Perpendicular to the MPBL the magnetic field is smoothly varying, unlike the sharp exponential decay in the ionosphere, and the characteristic decay scale appears larger too. All these effects are likely a result of Biot–Savart spatial integration smearing wave effects out further and/or ionospheric screening of smaller-scale features. However, despite the smearing/screening, perturbation amplitudes are still significant at several hundreds of nT, in agreement with observed candidate events (e.g. Kozyreva et al., 2019). The figure also demonstrates that the fields due to magnetopause (a and h) and Pedersen (b and i) currents cancel to within 1 nT accuracy. Fukushima's (1976) theorem is thus an excellent approximation in our model. While not strictly valid, since the surface wave's magnetopause currents are not semi-infinite FACs, the scale of variation along the field in the magnetosphere are much longer than those across the field (k_y≫k_z), likely making it a reasonable approximation. In contrast, Archer et al. (2023) found greater discrepancies in a box model when using several R_E perpendicular wavelengths, representative of magnetospheric scales. Given, however, that the ground field rotation by the ionosphere still significantly differed from 90° means this cannot result from the violation of Fukushima's (1976) theorem, but may instead be due to the non-planar wave structure (Hughes, 1974; Hughes and Southwood, 1974). In global MHD simulations of large-scale magnetopause surface waves, while the ground magnetic field is principally due to Hall currents, only a complex superposition of all current systems fully describes the total variance (Archer et al., 2023).

Figure 7 Comparison of the ground magnetic field with far-field theory in a similar format to Fig. 5. (a–f) Example comparison across a single run in cylindrical coordinates. Shown as a function of perpendicular radial distance are the real (a, c, e) and imaginary (b, d, f) parts of the radial (a, b), azimuthal (c, d), and vertical (e, f) components divided by trigonometric functions based on theory for various azimuths (colours). (g, h) Comparison across runs varying model input parameters. The imaginary parts of the field are shown along y=0. Inputs varied (indicated by colour) are the number of oscillations (g, k), wavelength (h, l), and Pedersen (i, m) and Hall conductances (j, n). Theoretical solutions are shown (dashed lines) in their full form (grey) and the R≫h asymptotic limit (black).

We compare the numerical results with far-field theory (Sect. 3.2) in Fig. 7. Firstly we compare results in cylindrical coordinates from a single model run across a range of polar angles. Plotted as a function of perpendicular radius are the three components of the ground magnetic field due to Hall currents, with panels (a), (c), (e) depicting the real part and (b), (d), (f) the imaginary. For distances greater than the wave's extent, the field tends towards zero in the real and a $\mathrm{1}/R^{\mathrm{2}}$ dependence in the imaginary. This is in good agreement with theory for R≫h (black dashed), necessarily the case here as there is little difference when this assumption is relaxed (grey dashed). As was found in the ionosphere, the ground magnetic field converges on the theoretical curves most quickly along the y axis. We then study how model inputs affect the far-field ground magnetic perturbations, again varying each input individually across panels (h)–(j) and normalising results based on theory in panels (k)–(n). Throughout the imaginary part of the y component along the x axis is shown. Like with the ionospheric potential, the direction of the magnetic field alternates between odd and even N (panel g). We again normalise the distance by N, such that wave packets' extents are identical and depicted by the vertical dotted line in panel (k)–(n). Theory dictates that, unlike with the ionosphere, this should leave a residual $\mathrm{1}/N$ dependence on the ground magnetic field. Panel (k) accounts for this and shows the curves for all runs do tend towards theory (black dashed). Ground magnetic perturbations increase in amplitude with wavelength (panel h) and are well-accounted for by normalising distances (panel l), again giving good agreement with theory. Finally reducing the Pedersen (panel i) and/or increasing the Hall (panel j) conductances lead to increased ground fields, with these being simply reciprocal and multiplicative factors (panels m–n) as expected. Thus our numerical results are shown to be in excellent agreement with far-field theory.

Figure 8 Characterising the external ground magnetic field for pseudoinfinite waves. (a) Amplitude at x=0 of the x (blue crosses), y (red circles), and z (yellow squares) components as a function of wavenumber. (b) Half width at half maximum along the x axis for the three components as a function of wavelength. In both panels solid lines indicate fits to the data with shaded areas showing their 95 % prediction intervals. Bottom panels show the data along the x axis for all wavenumbers/wavelengths and normalised by both amplitude and half width at half maximum for the x (c), y (d), and z (e) components as coloured dots. For the x component the analytic envelope from a Hilbert transform is also shown (grey). Illustrative bell-shaped functional fits to the data are depicted by dashed lines.

We now move to the near-field, where no analytic form could be determined, modelling pseudoinfinite waves by setting the number of cycles to the maximum possible in the domain N_max for each run. Throughout Fig. 8 the effect of the wavelength on the (Hall) ground magnetic field fluctuations is investigated. Panel (a) shows through markers the amplitudes of all three components at the origin multiplied by the conductance ratio. For the x component, which follows an odd function passing through zero at x=0 (Fig. 6j), we use the analytic signal envelope via a Hilbert transform in x to extract a suitable amplitude (see Fig. 8c). Across all wavelengths considered, amplitudes are lower than those in the magnetosphere (given by Plaschke and Glassmeier, 2011, and depicted as the horizontal dashed line), even when considering typical conductances (Ridley et al., 2004). All three components appear to follow an exponential decay law with k_y, akin to the relation of Hughes and Southwood (1976) $\exp (-kh){\mathrm{\Sigma }}_{\mathrm{H}}/{\mathrm{\Sigma }}_{\mathrm{P}}$ for infinite plane Alfvén waves where the e-folding scale is the ionospheric altitude. We, therefore, fit a decaying exponential to the data for each component using nonlinear least squares with robust iterative bisquare weighting (Huber, 1981). To account for systematic errors in the numerical modelling, a constant term is included in the fit. The results are shown as the solid lines along with their 95 % prediction intervals (shaded areas) and fit coefficients are listed in the figure. These reveal ground field amplitudes of up to a few thousand nT may be possible in the limit of large-scale waves. Note, while k_y=0 results in no FACs hence no ionospheric or ground signatures, the FAC amplitude quickly becomes within 5 % of J₀ for k_y≈3k_z, only slightly altering the upper limit of ground amplitudes. Signals are strongest in the vertical component for all wavelengths (∼20–40 % larger than ${\stackrel{\mathrm{̃}}{B}}_{\mathrm{H},y}$). This component is generally associated with FACs if ground inductive effects can be neglected – Eq. (16) demonstrates the magnetic field due to Hall currents resulting from an isolated FAC are significantly stronger in the vertical, exhibiting a prominant peak directly below the FAC. Therefore, given a surface wave constitutes a string of alternating FACs, it is not surprising that the vertical field component is strongest. Interestingly, the e-folding scale with k_y is larger than the ionospheric altitude and varies for each component, ranging between ∼140–180 km. The Hughes and Southwood (1976) relation was derived requiring electric and magnetic fields in the atmosphere, ionosphere, and magnetosphere follow the same plane wave ansatz, clearly not the case for a surface wave. In theory one can decompose any problem into its consituent wavenumbers and still apply their result. Based on this, a larger decay scale would be expected, given that surface waves necessarily exhibit a spectrum of total wavenumbers $k=\sqrt{k_x^{\mathrm{2}}+k_y^{\mathrm{2}}}$. However, our numerical work highlights that the interference pattern between the different wavenumbers, even for this relatively simple example, can yield unexpected results. Heuristically, the magnitude of these decay scales can be rationalised through the Fourier transform of the ionospheric solutions

$${\mathcal{F}}_x\left[{\mathrm{e}}^{-\left|k_{y}x\right|}\right]\left(k_x\right)=\sqrt{{\displaystyle \frac{\mathrm{2}}{\mathit{\pi }}}}{\displaystyle \frac{\left|k_y\right|}{k_y^{\mathrm{2}}+k_x^{\mathrm{2}}}}$$

Calculating an average k_x for this distribution as the Fourier-weighted root-mean-square gives $\sqrt{\mathit{\pi }/\mathrm{2}}{k}_y$ corresponding to ∼176 km total e-folding scale from the magnitude of the full wave vector, broadly consistent with the numerical results. The surface wave has three distinct directions (along the magnetospheric field, and both along and perpendicular to the MPBL), whereas an infinite plane wave has only two (along either the field and the perpendicular wave vector component) and a single FAC only has one (the field direction). This means different behaviours might be expected for a surface wave in all three components. Indeed, this has already been seen in the ionosphere with sign changes across the MPBL in only one component of the two-dimensional currents. Therefore, it is perhaps not surprising that the three ground field components behave differently. An unexpected consequence though is the polarisation of waves as observed on the ground must become altered by the ionosphere from the circular polarisation in the magnetosphere (Plaschke and Glassmeier, 2011). Ground observed ULF waves are often assumed to provide remote sensing of the same waves originating in the magnetosphere, with the only ionospheric considerations being a 90° rotation and screening of smaller-scale waves (e.g. Ozeke et al., 2009). We have shown the introduction of even simple structure beyond a plane wave can dramatically change the one-to-one association from above to below the ionosphere.

For the same runs we also determine the Half Width at Half Maximum (HWHM) along the x-direction of the ground field perturbations, with results shown in Fig. 8b (for the x component this is again applied to the analytic signal envelope). While in both the magnetosphere and ionosphere a linear relation between HWHM and wavelength exists, as depicted by the dotted line, we find that not only are the scales larger on the ground but they appear to follow power laws. Given the form of the latitudinal Biot–Savart integral, this scaling is likely due to the Bessel function multiplication theorem (Abramowitz and Stegun, 2000)

$$k^{-\mathit{\alpha }}{K}_{\mathit{\alpha }}\left(kx\right)=\sum _{n=\mathrm{0}}^{\mathrm{\infty }}{\displaystyle \frac{(-\mathrm{1}{)}^{\mathit{\alpha }}}{\mathit{\alpha }\mathrm{!}}}{\left({\displaystyle \frac{(k^{\mathrm{2}}-\mathrm{1})x}{\mathrm{2}}}\right)}^{\mathit{\alpha }}{K}_{\mathit{\alpha }+n}\left(x\right),$$

a form of power law scale invariance preserving the fundamental spectral behaviour. Such behaviour is perhaps not surprising given dispersionless waves are scale invariant and power laws often result in self-similar natural models. We again apply fits, revealing the power law exponents are all less than unity and significantly differ between the three components. The latitudinal HWHM, the largest of the three components, while still greater than those in the magnetosphere and ionosphere range only ∼100–200 km even for reasonably large-scale waves up to ∼1000 km wavelengths (equivalent to m∼14 at 70° magnetic latitude). The vertical component has the smallest scales of ∼70–150 km. The different exponents for the two horizontal components means they approach one another for large wavelengths. Different latitudinal scales of amplitude variation across components of the ground magnetic field were also found in a global MHD simulation of large-scale surface waves (Archer et al., 2023), again with the vertical being smallest. Note, these resulting latitudinal scales comparable to the ionospheric altitude are not in violation of Hughes and Southwood (1976), since the longitudinal wavelengths considered are greater than h and only the total wavenumber $k=\sqrt{k_x^{\mathrm{2}}+k_y^{\mathrm{2}}}$ applies in ionospheric screening. Our results are a consequence of the interference between the different attenuated wavenumbers in x. Interestingly, these latitudinal HWHM are smaller than the typical ∼200–300 km closest separations of ground magnetometer chains (Archer et al., 2024a) and typical global MHD ground-based output grids (Brenner et al., 2025), meaning such networks/simulations might not well resolve mesoscale surface waves. The different power laws for each component result in further polarisation changes on the ground compared to space.

For all the runs depicted in Fig. 8 we normalise the magnetic field profiles according to their amplitude and HWHM, shown as the dots in panels (c)–(e), which collapses them onto one another. More suprisingly, this normalisation renders all three components following the same profile (though in the case of the x component this applies to the analytic signal envelope).

Figure 9 Characterising the external ground magnetic field horizontal polarisation for pseudoinfinite waves. (a–c) The three polarisation Stokes parameters normalised by the local intensity plotted along y=0 for different wavenumbers (colours). Inset diagrams indicate corresponding polarisations for each sign. (d) Ellipticity of the wave polarisation along y=0. (e) Half width at half maximum of the ellipticity as a function of wavelength. The solid line shows a power law fit to the data with the shaded area depicting the 95 % prediction interval.

We quantify the, thus far inferred, variation in the horizontal ground polarisation through the four Stokes (1852) parameters

$$\begin{array}{rl}{\displaystyle }{S}_{\mathrm{0}}=& {\displaystyle }|{\stackrel{\mathrm{̃}}{B}}_{\mathrm{H},x}{|}^{\mathrm{2}}+|{\stackrel{\mathrm{̃}}{B}}_{\mathrm{H},y}{|}^{\mathrm{2}}\\ {\displaystyle }{S}_{\mathrm{1}}=& {\displaystyle }|{\stackrel{\mathrm{̃}}{B}}_{\mathrm{H},x}{|}^{\mathrm{2}}-|{\stackrel{\mathrm{̃}}{B}}_{\mathrm{H},y}{|}^{\mathrm{2}}\\ {\displaystyle }{S}_{\mathrm{2}}=& {\displaystyle }\mathrm{2}\text{Re}\left({\stackrel{\mathrm{̃}}{B}}_{\mathrm{H},x}{\stackrel{\mathrm{̃}}{B}}_{\mathrm{H},y}^{*}\right)\\ {\displaystyle }{S}_{\mathrm{3}}=& {\displaystyle }-\mathrm{2}\text{Im}\left({\stackrel{\mathrm{̃}}{B}}_{\mathrm{H},x}{\stackrel{\mathrm{̃}}{B}}_{\mathrm{H},y}^{*}\right)\end{array}$$

S₀ is the intensity, S₁ corresponds to linearly polarized waves aligned with either coordinate axis, S₂ relates to linearly polarized waves oriented at 45° to the coordinate axes, and finally S₃ is associated with circularly polarised waves (Collett, 2005), as depicted by the red arrows in Fig. 9a–c. The degree of polarisation, $\sqrt{S_{\mathrm{1}}^{\mathrm{2}}+S_{\mathrm{2}}^{\mathrm{2}}+S_{\mathrm{3}}^{\mathrm{2}}}/S_{\mathrm{0}}$, was reassuringly unity to machine precision. Figure 9a–c plots the Stokes parameters normalised by intensity along y=0 for a range of wavenumbers, where we omit values where S₀<1 nT² due to potential numerical errors. Panel (c) shows S₂ is negligible, indeed calculating the polarisation ellipse's orientation angle $\mathrm{1}/\mathrm{2}\arctan (S_{\mathrm{2}}/S_{\mathrm{1}})$ yields zero everywhere, i.e. it is aligned with the x axis. $S_{\mathrm{1}}/S_{\mathrm{0}}$ has a minimum of −1 at x=0, demonstrating that near the MPBL the surface wave ground signatures are east–west linearly polarised, and slowly tends towards zero away from the MPBL with only slight wavelength variations. While $S_{\mathrm{3}}/S_{\mathrm{0}}$ is zero at the MPBL, it grows in magnitude towards ±1 with distance, again with minor differences across wavelengths. This indicates a reversal in polarisation occurs across the MPBL, as expected for a surface wave. We quantify this reversal through the ellipticity, the ratio of the minor to major axes (or equivalently the tangent of the ellipticity angle) of the polarisation ellipse calculated as $\mathrm{2}{S}_{\mathrm{3}}\sqrt{S_{\mathrm{1}}^{\mathrm{2}}+S_{\mathrm{2}}^{\mathrm{2}}}/(S_{\mathrm{1}}^{\mathrm{2}}+S_{\mathrm{2}}^{\mathrm{2}}-S_{\mathrm{3}}^{\mathrm{2}})$, shown in panel (d). This takes values from −1 (right-hand circular polarisation), through zero (linear), up to +1 left-hand circular), where handedness is with respect to the vertical, i.e. looking down on the ground. While this indeed shows a reversal in handedness across the MPBL and linear polarisation at x=0, this transition is gradual unlike in the magnetosphere where the infinitesimal current sheet leads to the discontinuous change indicated by the grey dashed line (Plaschke and Glassmeier, 2011). We quantify the scale over which this transition occurs, which appears to have a small dependence on wavelength, using the HWHM corresponding here to when the ellipticity crosses ±0.5 as indicated by the dotted horizontal lines. The inset panel (e) shows the HWHM to be ∼60–70 km, much smaller than typical ground magnetometer spacing and even comparable to those of localised magnetotelluric survey arrays (Archer et al., 2024a). As with the amplitudes, the ellipticity's HWHM follows a power law, though the exponent is rather small explaining why the dependence with wavelength is weak. The exponent is comparable to the difference in exponents of the two horizontal components' amplitude HWHM. It is worth noting that while the ellipticity has only a weak wavelength dependence, the scale over which the ellipticity changes is comparable to that of the amplitude for smaller-scale waves meaning the polarisation on the ground never appreciably reaches that in the magnetosphere. Therefore, as previously inferred, the ionosphere significantly modifies the polarisation of surface waves on the ground from that in the magnetosphere. This was qualitatively mentioned in the global MHD simulation of Archer et al. (2023), though the magnetometer output-resolution was insufficient to resolve the scale over which transitions in polarisation occurred across the ground. This result is not expected for infinite plane waves (Hughes and Southwood, 1976), hence is another example of how breaking this assumption leads to more complicated wave behaviour on the ground.

We tested all these results for pseudoinfinite waves against finite wave packets also (e.g. Fig. 6). After the fluctuations were extracted by the spline-based filtering, we found the ground fields agreed very well with the pseudoinfinite wave results, particularly when close to the centre of the wave packet. This is likely a result of the fluctuating part of the ionospheric currents for finite waves also being well approximated by the infinite wave case, as was shown earlier. Consequently, we have completed a comprehensive characterisation of magnetopause surface waves' ground magnetic field in our model when considering only MI-currents.

4.2.2 Magnetic field (telluric currents)

We use CIM to calculate the contribution to the ground magnetic field from induced telluric currents, which only depends on the skin depth (Eq. 10). While skin depth varies with two model inputs, surface wave frequency and ground conductivity, we only need consider ground-induced effects over a representative range of skin depths. Figure 10 indicates how skin depth varies with both conductivity and frequency, where typical ranges for different ground types and wave phenomena are also given. This shows that magnetopause surface eigenmodes' (MSE) low frequencies (Archer and Plaschke, 2015) result in skin depths much larger than the ionospheric altitude, making telluric current contributions likely small for most ground types as previously commented by Archer et al. (2023). In contrast, telluric currents may be more important for Kelvin–Helmholtz waves, due to higher frequencies, as skin depths become comparable to h. Finally, the well-known phenomenon that the presence of salt water greatly affects ground magnetic field observations is apparent in this figure given skin depths of only a few kilometers across all frequencies.

Figure 10 Variation of skin depth magnitude with ground conductivity and surface wave frequency depicted as both colours and logarithmically-spaced contours. Ranges of typical conductivities and frequencies for different ground types and magnetopause surface wave phenomena, respectively, are also indicated.

CIM assumes oscillatory currents, which introduces differences between propagating and stationary surface wave packets. With stationary waves all perturbations oscillate in time. However, only the periodic fluctuations are oscillatory for propagating waves, since background trends lead to much slower temporal variations as the wave packet propagates with speed $\mathit{\omega }/k_y$. Our implementation of CIM does not take this into account, simply integrating over all currents within the model. However, given our model is linear with dispersionless waves, we can spatially filter following CIM to correct the internal ground magnetic fields for propagating waves. The results for a finite wave packet are shown in panels (d)–(e) and (k)–(l) of Fig. 6. Amplitudes of fluctuations are identical between propagating and stationary waves, with a 90° phase difference present. However, since the trends also contribute for stationary waves they have significantly larger internal field contributions, which to our knowledge has not been suggested before. While these finite wave effects appear much larger than the fluctuations in the near-field, we leave modelling them to future work focusing only on pseudoinfinite propagating waves.

Figure 11 Characterising the induced ground magnetic field for pseudoinfinite propagating waves as a function of skin depth for different wavelengths (colours). A similar format of pairs of rows to that in Fig. 5 is used. The top rows (a–c, g–i, m–o) show the model output, whereas bottom rows (d–f, j–l, p–r) normalise results based on the external ground magnetic field such that curves collapse on to one another. Empirical fits (black) and their 95 % prediction intervals are also shown for these normalised results. Columns depict results for each magnetic field component separately. (a–f) Amplitudes at x=0. (g–l) Half width at half maximum along the x axis. (m–r) Phase difference from the external ground field at x=0.

Figure 11a–c show internal magnetic field amplitudes of the three components on the ground as a function of skin depth for five different wavelengths. As before, the x component corresponds to an analytic envelope. All three components appear to follow exponential decays with different decay constants for each wavelength. However, as the field becomes small results level off, suggesting numerical errors begin to dominate. We wish to normalise these curves on to one another in order to characterise the response. An obvious choice of normalisation for the distance is the wavenumber, as used previously. For the amplitude though, we recognise in the limit of infinite ground conductivity (zero skin depth) all components of the internal field must be equal in magnitude to that from Hall currents. This is because in a perfect conductor the induced currents act to cancel out the vertical magnetic field at the surface, in turn doubling the horizontal field. Therefore, we divide the internal field amplitudes by those due to Hall currents, with the results shown in panels (d)–(f). This succeeds at normalising the numerical results onto one another and we fit an exponential decay law to all normalised data. While the amplitudes and offsets from the three fits are consistent with unity and zero, respectively, we again find decay coefficients differs for each component.

We next investigate the HWHM along x for each component. Given numerical errors for large skin depths, we only consider when the ratio of internal to Hall amplitudes was greater than $\mathrm{5}\times {\mathrm{10}}^{-\mathrm{3}}$ (cf. Fig. 11d–f). The results in Fig. 11g–i show in the limit of zero skin depth the HWHM plateau to a minimum value (different for each wavelength), while as the skin depth increases so too do the HWHM. This makes sense as in CIM image Hall currents get placed further away from the ground as the skin depth increases, thereby widening the scale of their ground signatures. The minimum HWHM should, for similar reasons as the amplitude, be equal to that due to Hall currents, therefore, in panels (j)–(l) we plot the ratio of the two scales minus one. Given the HWHM for the Hall currents followed a power law, we show a log–log plot of this normalised HWHM against skin depth. This again causes the different wavelength curves to collapse onto one another and fitting reveals the power law exponent for each component is consistent with unity, i.e. linear relationships.

Finally, we show the phase difference between the internal field and that from Hall currents in Fig. 11m–o, again only when the ratio of the internal to Hall amplitudes was sufficiently large. For small skin depths phase differences are close to zero for the horizontal components and −π for the vertical, as expected in the perfectly conducting limit. As skin depth initially increases, so too do the phase differences, though at different rates for the various wavelengths with this being greater for larger wavelengths. Normalising the skin depth using the wavenumber and unwrapping the phase, we find in panels (p)–(r) that the phase differences collapse onto one another and appear to follow linear relations.

Similarly to the ground field due to MI-currents, and likely for the same reason, the fluctuating part of the internal field from finite wave packets (Fig. 6) agree well with the infinite wave case. Thus concludes our characterisation of the ground magnetic field contributions due to telluric currents.

4.2.3 Geoelectric field

CIM is again used to calculate the geoelectric field, with example cuts shown in Fig. 6f–g and m–n for both propagating and stationary wave cases. The geoelectric field is generally greatest in the latitudinal x-direction peaking at the MPBL-projection, with the longitudinal component reversing sign across the MPBL and being significantly weaker in amplitude. Like with the ground magnetic field, the geoelectric field has the same amplitudes of periodic fluctuations between both propagating and stationary cases with a 90° phase shift between them (note y=0 is an E_g,y node for stationary waves). Furthermore, stationary waves have significantly larger overall amplitudes due to finite wave effects' contributions, i.e. the trends away from periodic behaviour. Here we focus only on the periodic geoelectric field fluctutations for pseudoinfinite propagating waves.

The geoelectric field depends on both wave frequency and ground conductivity independently, i.e. the dependence cannot be reduced to a single joint parameter. However, like with the induced magnetic field, the vector potential from which the geoelectric field is derived is a function of only skin depth (Eq. 13), meaning for fixed skin depth the quantity ${\mathit{E}}_{\mathrm{g}}/\mathit{\omega }$ is invariant (Eq. 12), simplifying the analysis.

Figure 12 Characterising the geoelectric field for pseudoinfinite propagating waves in the limit of zero ground conductivity. Format is similar to Fig. 8. Amplitudes of the x (a) and y (b) electric field components for different surface wave frequencies are depicted by colours, with the frequency-normalised amplitudes shown below (c). Asymptotic limits of the half width at half maximum for infinite ground conductivity is also shown in panel (d) as the lighter colours.

We first consider the zero ground conductivity (infinite skin depth) limit, which will result in the largest geoelectric fields. Figure 12a–b shows how amplitudes of the two geoelectric field components at the MPBL vary with wavenumber for a range of representative frequencies. For the y component an analytic signal envelope has been taken. The amplitude increases with frequency in this limit, demonstrating significant geoelectric fields of up to several mV m⁻¹ may be possible, particularly for large-scale waves and higher frequencies. Both components show approximately exponential decay with wavenumber, which is parameterised in panel (c) where the wave frequency and ionospheric conductance ratio dependences have been removed. The e-folding scale is different for the two components and is significantly greater than both the ionospheric altitude and the e-folding scales of the ground magnetic field. The HWHM again appears to follow a power law with wavelength, with the scales of the geoelectric field again being greater than that in the magnetosphere/ionosphere as well as that of the ground magnetic field, being in the range ∼80–300 km. This highlights that the local $\partial {\mathit{B}}_{\mathrm{g}}/\partial t$ is insufficient in prescribing the geoelectric field for mesoscale phenomena. The longitudinal geoelectric field, E_g,y, while weaker in amplitude is the more widespread geoelectric field component. In panels (e)–(f) all profiles are normalised by their respective amplitudes and HWHM. Curiously, the geoelectric field follows the same normalised profile as the ground magnetic field. This is not obvious mathematically as while the geoelectric field is derived from the same current distribution, the weighting by distance of these currents is different.

Figure 13 Characterising the geoelectric field for pseudoinfinite propagating waves as a function of skin depth for different wavelengths (colours) in a similar format to Fig. 11. Subscripts of zero and infinity refer to the ground conductivity limits. Illustrative fits (black) use a sigmoid-shaped function.

Having characterised the geoelectric field for the limiting case of zero ground conductivity, and knowing that infinite ground conductivity will yield no geoelectric field, we wish to understand the transitional behaviour between these limits. As intimated earlier, we do this by normalising the geoelectric field by the frequency and conductance ratio such that the result is only dependent on skin depth. We first investigate amplitudes, which are shown in Fig. 13a–b for the two components as a function of skin depth. This shows for small skin depths the geoelectric field is small, but as skin depth increases it quickly rises in magnitude. Once the skin depth becomes of order of the ionospheric altitude the geoelectric field plateaus towards the wavelength-dependent zero ground conductivity limits. We find in panels (c)–(d) that normalising the geoelectric field to its aysymptotic limit and the skin depth by the wavenumber causes all the curves to collapse onto a single sigmoid-shaped function.

Next we investigate the HWHM perpendicular to the MPBL. Panels (e)–(f) show this transitions between a slightly smaller scale for zero skin depth (infinite ground conductivity) to a larger one for infinite skin depth (zero ground conductivity), both of which vary with wavelength. While we have already parameterised the latter, the former cannot be ascertained directly as zero skin depth corresponds to no geoelectric field. Therefore, we numerically determine the aymptotic HWHM by assuming a geometric series, found to hold well for the smallest three skin depths shown. The results are displayed in Fig. 12b as the lighter colours. While the HWHM for infinite ground conductivity also follow power laws in wavelength, these have different coefficients and exponents to the zero ground conductivity case. In Fig. 13g–h we subtract this HWHM offset and divide by the difference in HWHM between the ground conductivity limits to yield normalised HWHM between zero and one. When also normalising the skin depth by the wavenumber the curves collapse onto another sigmoid-shaped form.

Finally, we compare the phase differences between the geoelectric field for finite skin depth and that in the limit of zero ground conductivity (infinite skin depth), as shown in Fig. 13i–j. While as expected these all tend towards zero for large skin depth, there is also a wavelength-independent limit for small skin depths of $-\mathit{\pi }/\mathrm{4}$. However, how quickly the phase converges on these limits is wavelength-dependent. Normalising the skin depth through the wavenumber again causes the data to collapse onto a single curve for each component.

Again, finite waves' fluctuations (Fig. 6) show good agreement with the infinite wave case. We, therefore, have fully parameterised the geoelectric field due to magnetopause surface waves.

5 Discussion

5.1 Summary

We have performed a comprehensive characterisation of expected effects of mesoscale magnetopause surface waves (ionospheric wavelengths ranging 200–1600 km) within the ionosphere and on the ground using a highly simplified local numerical model, detailing how the response's amplitudes and spatial scales vary with key wave and system properties. We have focused on the near- and far-field limits of finite constant-amplitude sinusoidal wave packets of limited extent, motivated by the existence of both quasi-continuous (e.g. Fairfield et al., 2000; Viall et al., 2009) and impulsive (e.g. Shue et al., 2009; Kozyreva et al., 2019) surface wave drivers. These have resulted in numerous novel insights into the problem of surface wave coupling across the MIG system. We summarise the major findings of this study:

• 1.

  Finite surface wave packets' ionospheric and ground signatures can be described as the sum of periodic fluctuations and slowly-varying background trends.

  • a.

    Within the wave packet, fluctuations are modulated sinusoids along the wave vector direction with the same wavelength as the original surface wave, which are not present far outside the wave packet longitudinally.

  • b.

    Trends are a result of the finite nature of the wave packet, as imbalances in source terms on either side along the wave vector direction result in imperfect cancellation.

  • c.

    These trends for finite wave packets introduce differences in effects for propagating and stationary surface waves. With propagating waves, trends cause slower-than-ULF time variations meaning only the fluctuations contribute to the geoelectric field and ground-induced contribution to the surface magnetic field. In contrast, for stationary waves the trends also oscillate at the surface wave frequency, leading to larger overall amplitudes in both quantities.

• 2.

  Close to the wave packet (the near-field) ionospheric fluctuations agree well with theory for an infinite surface wave when $k_{y}R/N$ is small, showing exponentially decaying signatures with latitudinal distance from the projected Magnetopause Boundary Layer (MPBL) flux tubes. Outside the wave packet where $k_{y}R/N$ is large (the far-field), surface waves acts as a localised string of field-aligned current sources causing large-scale ionospheric electrodynamics. Their influence on the ionosphere and ground magnetic field can be described as a two-dimensional current dipole, where the dipole moment is larger for greater surface wave amplitudes, wavelengths, and number of oscillations.

• 3.

  The relationship between the magnetic field above and below the ionosphere becomes more complicated for surface waves than infinite plane Alfvén waves (Hughes, 1974; Hughes and Southwood, 1974, 1976) due to interference patterns between the spectrum of wavenumbers introduced by the additional spatial wave structure.

  • a.

    The ionosphere similarly screens small-scale surface waves from the ground, however the e-folding scale (which varies with each component of the field) of this effect with azimuthal wavenumber is larger than the ionospheric altitude, likely due to a spectrum of total wavenumbers present.

  • b.

    While the rotation of the ground magnetic field is on average close to (but not exactly) 90°, significant spread in rotation angles occurs spatially. Instances of the opposite rotation sense also occur periodically along the MPBL-projection.

  • c.

    The ionosphere significantly modifies the polarisation of the surface wave, only gradually transitioning between handedness across the MPBL-projection, unlike the discontinuous change in the magnetosphere. For smaller wavelengths this means ground-observed polarisations may never match that in space.

• 4.

  Surface wave effects are highly localised latitudinally. While in the ionosphere their latitudinal scales vary linearly with wavelength as in the magnetosphere, the Half Width at Half Maxima (HWHM) for the ground fields appear to instead follow power laws with exponents that differ by component.

  • a.

    Ionospheric HWHM range ∼20–175 km meaning even relatively large-scale surface waves should have very concentrated impacts latitudinally.

  • b.

    Ground magnetic field HWHM, though larger than in the ionosphere and magnetosphere, are comparable to or smaller than typical spacings of ground magnetometer chains for even large wavelength surface waves at ∼100–200 km.

  • c.

    Ground magnetic polarisation changes occur over distances comparable to the spacings of even localised magnetotelluric survey arrays at ∼60–70 km.

  • d.

    Geoelectric field scales are larger than those of the ground magnetic field at ∼80–300 km, highlighting that local $\partial {\mathit{B}}_{\mathrm{g}}/\partial t$ is insufficient in prescribing the geoelectric field for mesoscale phenomena.

• 5.

  Inductive effects in the ground and their contribution to the ground magnetic and geoelectric fields depend on the magnitude of the wavenumber multiplied by the skin depth, i.e. $\left|k_y\mathit{\delta }\right|$. Across the ranges of wavenumbers and skin depths applicable to magnetopause surface waves, the ground can act either like:

  • a.

    A perfect insulator ($\left|k_y\mathit{\delta }\right|\gg \mathrm{1}$), where the vertical component of the ground magnetic field acts a strong detector of surface wave's FACs and geoelectric fields are maximised.

  • b.

    A perfect conductor ($\left|k_y\mathit{\delta }\right|\ll \mathrm{1}$), where telluric currents cancel the vertical component of the magnetic field, double the horizontal ground field components, and no geoelectric field is induced.

  • c.

    A finite conductor, where properties smoothly transition between the above limiting cases with normalised skin depth $\left|k_y\mathit{\delta }\right|$.

• 6.

  The predicted amplitudes of magnetopause surface wave effects suggest they may act as a significant, but highly latitudinally localised, source of Joule heating in the ionosphere (up to hundreds of mW m⁻²) and the geoelectric fields (up to several mV m⁻¹) that drive geomagnetically-induced currents, particularly for large-scale waves, thereby contributing to space weather impacts.

5.2 Limitations

These results have been determined from a deliberately simple local numerical model, enabling such comprehensive characterisation of how surface wave effects may vary with key wave and system parameters. However, as we outlined from the outset, the simplicity of this model is designed to enable initial physical insight and intuitive understanding to the coupling of surface waves across the MIG-system. It is, therefore, important to consider the validity of the model assumptions along with outlining further refinements that could test incorporating more realistic effects in the future to iterate these results.

We have used a box model magnetosphere with vertical field lines. While we justified this does not affect current amplitudes at the MI-interface, meaning that ionospheric current systems are reliable, the assumption may affect the ground fields that integrate the magnetospheric currents. Our results for the ground fields pertain solely to ionospheric Hall currents and should only be treated as such. This is because Fukushima's (1976) theorem was shown to be a good approximation due to the model setup. In reality though this approximation is only considered reasonable for high-latitudes (Untiedt and Baumjohann, 1993; Laundal et al., 2015) and, in general, one would also expect a contribution due to the non-cancellation of magnetospheric and Pedersen currents, with several factors to consider. Firstly, the field-aligned phase variation in our model is linear with distance. However, inhomogeneity in magnetic field and plasma condition along field lines may modify parallel wave structure. Archer and Plaschke (2015) demonstrated surface wave phase speeds are significantly larger close to Earth, primarily due to the dipole field, such that very slow phase variation are expected along the field. Even in our simple model the ground fields were insensitive to field-aligned phase variations, given $k_{⟂}\gg k_{\parallel }$ and currents closest to Earth dominate the Biot–Savart law (only 5 % of the ground field contribution arises from magnetopause currents above ∼0.7–4 R_E, dependent on wavelength). Thus one might similarly expect little sensitivity to slow field-aligned phase variations, though this could be checked. Secondly, current density variations along field lines are also expected due to the changing area of flux tubes along the field, not accounted for in a uniform field. Given the product of the current density and flux tube area, jdA, remains invariant along flux tubes though, this aspect of the Biot–Savart integrand is little affected to a first approximation. However, the magnetic geometry does alter the integrand via the position vectors. Fukushima (1976) outlines a correction term for realistic geometries using an additional current system along these magnetic field lines and cancelling out the vertical ones. Examining this correction for surface waves and its dependence on latitude and wavelength should be investigated.

In our model we only consider wavelengths up to of order the distance to the horizon from ionospheric altitude, for which the Cartesian model used is appropriate. For larger wavelengths a spherical model would likely be required, which again would introduce latitudinal-dependence to results.

We justified using magnetopause currents with infinitesimal width in terms of the ionospheric currents outside the thin boundary layer projection. A simple modification to account for finite thickness is to linearly interpolate ionospheric currents in the range $-d/\mathrm{2}<x<d/\mathrm{2}$, corresponding to a uniform magnetopause current layer of typical thickness d∼40 km at the ionosphere. Given the highly localised ionospheric response, maxima in the ionosphere become significantly capped for small-scale waves (by ∼50 % for 200 km wavelength), whereas large-scale waves are little affected. In terms of the ground field fluctuations, we find accounting for the boundary thickness overall rescales the entire perturbation strengths by a factor of order $\exp (-d/\mathit{\lambda })$, i.e. less than a 20 % reduction for small-scale waves, and does not seemingly affect the HWHM or latitudinal-form. Background trends also remain unaffected. While in reality the current layer is likely not uniform, instead exhibiting a peak near the inner/equatorward edge (Archer et al., 2023), a uniform current layer provides the maximum possible difference in results from our assumed infinitesimal layer when considering boundary layer structure.

The magnetopause surface waves currents were derived from Ideal MHD (Plaschke and Glassmeier, 2011). Therefore, we include no parallel electric fields between the magnetosphere and ionosphere. The result was FAC amplitudes independent of wavelength. Johnson et al. (2021), in contrast, employed a linear Knight (1973) relation in their MI-coupling, finding a nonlinear Kelvin–Helmholtz vortex size which maximises FACs. Following the framework of Yoshikawa et al. (2011), it can be shown in our model for dusk-sector linear waves with small current amplitudes (compared to the background Region-1 currents) the FAC amplitude is slightly modified becoming $J_{\mathrm{0}}(\mathrm{1}+k_y^{\mathrm{2}}{\mathrm{\Sigma }}_{\mathrm{P}}/\mathit{\kappa })$, where κ is the Knight parameter. The resistive-coupling auroral scale length $\sqrt{{\mathrm{\Sigma }}_{\mathrm{P}}/\mathit{\kappa }}$ is typically ∼50 km (Lyons, 1980), implying auroral acceleration should not significantly modify linear surface wave results for the wavelengths presented. If the resistive-coupling length scale becomes comparable to or larger than the wavelength, so long as the linear regime still holds our results are only slightly modified given the predicted current amplitude modulation is weaker than the exponential and power law dependencies on wavenumber already found. However, if waves become nonlinear and/or develop into rolled-up vortices, then the parallel electric field likely would become important to the wave coupling, requiring a different approach.

While a uniform magnetic field and plasmas were employed in each half-space, perpendicular non-uniformities in the magnetic field or plasma necessarily couple compressional waves to Alfvén modes throughout the magnetosphere (Radoski, 1971). Archer et al. (2023) showed in high-resolution global MHD simulations this leads to additional FACs on closed field lines throughout the magnetosphere that contribute to the ionospheric and ground responses, which our model does not capture. Nonetheless, the simulation showed FAC amplitudes exhibited distinct peaks within the projection of the magnetopause boundary layer, with those from the non-resonant coupling being around an order of magnitude weaker. Therefore, one might expect the overall response to be dominated by the magnetopause currents that our model does capture. Future work could investigate introducing these additional currents and their effects.

Only uniform ionospheric conductances were used. It is known that solar illumination and auroral precipitating particles contribute to the spatial distribution of conductances (Ridley et al., 2004), which will affect current patterns. Sunlight effects result in small variations over scales comparable to the planet's surface, i.e. much larger than a surface mode wavelength, hence will likely have little effect on predictions. The auroral oval typically has a latitudinal width of 10–17° (Walach et al., 2017), i.e. of similar size to or greater than the model domain and thus larger than the scales of surface waves considered. Proportional changes in conductance are also typically smaller on the dayside than the nightside (Ridley et al., 2004). This suggests mesoscale dayside currents due to surface waves might be little affected by the additional auroral conductance. Laundal et al. (2015) also argue that at the dayside polar cap boundaries ionospheric conductance gradients are approximately perpendicular to equipotential contours, meaning only Hall currents should contribute to the ground magnetic field as per Fukushima (1976). Further work could modify the model to include conductance gradients and investigate their effects on surface waves over different scales.

We have also only considered uniform ground conductivities, which results in a simple relationship for the skin depth in terms of both wave frequency and conductivity. However, CIM can also be used for vertically-varying ground conductivities, i.e. a multi-layered Earth (Thomson and Weaver, 1975). In this case, the effective skin depth used in CIM for placing image currents becomes altered from the uniform case with a more complicated frequency dependence. Since we showed our results for surface waves are a function only of the skin depth, they can readily be applied to these cases also. However, where three-dimensional conductivities are thought to be important (e.g. Bedrosian and Love, 2015) a full inductive treatment would likely be required (Weaver, 1971).

Finally, our model contains no self-consistent dissipation, nor additional non-Ideal physics such as ionospheric viscosity (Yerg, 1952), turbulence (Guio and Pécseli, 2021), or kinetic plasma effects (e.g. Hasegawa, 1976; Lee et al., 1994). This means amplitudes presented are likely overestimates. In contrast, global MHD simulations usually suffer from strong numerical diffusion artificially damping waves and smearing their effects over larger scales (Claudepierre et al., 2009; Hartinger et al., 2015b), leading to underestimates. Given the large amplitudes and short scales predicted from this simple model, it is clear future simulations using numerical schemes and grid geometries that limit diffusion at the magnetopause/cusps (e.g. Sorathia et al., 2017; Zhang et al., 2019) run at even higher spatial resolution throughout the magnetosphere, ionosphere and on the ground are required to better test results derived from such a simple model against a much more representative environment.

5.3 Implications

The results presented in this work provide insight into the ionospheric and ground-based impacts magnetopause surface waves may have, particularly how their amplitudes and spatial scales may vary. Archer et al. (2024a) recently laid out future observational directions for unveiling magnetopause dynamics and their geospace impacts, with the wealth of current ground-based instrumentation being highlighted in particular. Correctly interpreting ground-based measurements as being due to magnetopause surface waves though requires comprehensive theoretical underpinning. Reanalysis of previously examined events in light of this work is thus needed, as well as dedicated future campaigns taking advantage of recent and upcoming capabilities.

Doppler radars, such as the SuperDARN network (e.g. Nishitani et al., 2019), provide measurements of the line-of-sight component to the ionospheric velocity. Historically ULF wave studies have been limited to measurements from a single radar (e.g. Fenrich et al., 1995) due to ∼1–2 min full scan time. However, SuperDARN's range gates have typical spatial resolution ∼45 km (Chisham et al., 2007), similar to surface waves' latitudinal scales for smaller wavelengths (∼200–400 km at the MI-interface, ∼1–2 R_E at the equatorial magnetopause, m∼30–70). It is thus questionable whether SuperDARN might resolve these features. The discretionarily-available “myopic” mode can provide ∼15 km resolution but with smaller field-of-view (Lester et al., 2004) and may be required for dedicated surface wave studies. On the other hand, availability of only the line-of-sight velocity component is limiting since the look direction might not always be best aligned with wave phenomena or provide sufficient discrimination against different wave modes. Ongoing upgrades to the SuperDARN network to digital systems are now enabling full scans over only 3.5 s, which removes the restriction of only single radar (McWilliams et al., 2023). While typical spatial resolution is unchanged, the higher ∼15 km resolution is also still available. Therefore, it may soon be possible, thanks to overlapping fields-of-view from different digital radar and applying high-order local/regional fitting (Fenrich et al., 2019), to reconstruct the ionospheric convection vortices and current dipole patterns of magnetopause surface waves, though this would likely need to be as part of a dedicated campaign. The ability, or not, of radar to resolve surface waves' ionospheric effects has implications not only on remote sensing magnetospheric processes, but on overall estimates of Joule heating in the ionosphere/thermosphere and hence global space weather. Hurd and Larsen (2016) previously estimated small-scale fluctuating ionospheric electric fields, poorly resolved by radar measurements, might contribute more than twice the contributions to Joule heating of the large-scale steady fields. Our model suggests magnetopause surface waves could be a significant source of highly latitudinally-localised oscillatory Joule heating. This strengthens the need for high spatio-temporal resolution ionospheric coverage to better understand how mesoscale phenomena contribute to global MI-coupling and space weather (Heelis and Maute, 2020; Pakhotin and Mann, 2024).

Our ground magnetic field results suggest mesoscale structures, like magnetopause surface waves, can have significant amplitudes but may be poorly resolved by the typical spacing of ground magnetometer networks. This may be why conclusive ground-based evidence for surface modes has proven challenging (Pilipenko et al., 2017, 2018; Kozyreva et al., 2019). These results come from a simple homogeneous model, without the complication of spatially-varying ground conductivities that are often invoked when closely-spaced (<100 km) magnetometers measurements differ (e.g. Juusola et al., 2020; Shi et al., 2025b). Our work, therefore, motivates the need for denser networks of ground measurements in order to identify the polarisation and wave power variations predicted. Another possibility might be placing magnetotelluric survey networks, smaller dense arrays taking simultaneous geoelectric and geomagnetic field measurements (Hartinger et al., 2020; Shi et al., 2022), close to the expected position of the OCB as part of dedicated campaigns. Our work has also highlighted the value of studying the vertical component of the ground magnetic field. Firstly, this component acts as a clear signifier of field-aligned currents, when ground inductive effects can be neglected. However, given for surface waves the expected horizontal and vertical field components are of comparable magnitude, this suggests that the ratio of the vertical to horizontal components' amplitudes could be used as a measure of the relative contribution of ground induction in ground-based magnetic field observations.

The ESA/CAS Solar Wind Magnetosphere Ionosphere Link Explorer (SMILE) mission aims to image the dayside magnetopause in soft X-rays, intending to uncover the fundamental modes of the dayside solar wind–magnetosphere interaction through near-continual remote sensing (Wang et al., 2025). While signal-to-noise is expected to be difficult under typical driving, it is hoped large-scale magnetopause motion such as directly-driven surface waves might be resolvable (Connor et al., 2025), though this may require advanced processing techniques (Archer et al., 2024a). As part of the SMILE mission, coordinated efforts to complement space-based measurements with ground-based instruments are being undertaken (Carter et al., 2024), including software tools for planning dedicated ground-based campaigns (Walach et al., 2024). These might address open science questions such as the scale sizes and structuring of ULF oscillations in the ionospheric cusp region and how these map out and relate to wave phenomena in space, such as magnetopause surface waves. Near-continual remote-sensing of magnetopause surface waves in soft X-rays with simultaneous conjugate measurements of the ionosphere and ground would likely provide the best test of these theoretical results, eliminating the need for rare in situ conjunction events.

Finally, the framework and scope of the numerical model presented here might be extended. Plasmapause surface wave (Chen and Hasegawa, 1974; He et al., 2020) impacts require only scaling some of the fixed model parameters, therefore, many of the results already presented might be readily applied. Adjusting the form of the currents in the magnetosphere to other mesoscale, time-varying phenomena, e.g. field line resonances (Southwood, 1974; Wright and Elsden, 2020) or flux transfer events (Lee and Fu, 1985), would enable them to be studied in a comparative way to better understand the different sources of ionospheric and ground perturbations in general, as well as their relative importance within the context of space weather.

Appendix A: Spline-based filtering

Here we describe a spline-based filtering method for splitting data, a(y), from our model into periodic fluctations, $\stackrel{\mathrm{̃}}{a}\left(y\right)$, plus a background trend, 〈a(y)〉. It was inspired by Empirical Mode Decomposition (EMD; Huang et al., 1971), which algorithmically decomposes non-stationary, non-linear data into well-behaved oscillatory components plus a non-oscillatory residual/trend. Our filtering method, depicted in Fig. A1, is performed along the y-direction only, where the periodicity of fluctuations is known a priori to be λ (a marked difference from EMD where such information is not known). For a given starting position, points of the same phase, i.e. separated by λ, are Akima (1970) spline interpolated across (red line in Fig. A1a). The same is performed for points in antiphase, i.e. distance $\mathit{\lambda }/\mathrm{2}$ from the originals (blue). The midpoint of this pairs of interpolants are an estimate of the background trend (green). For increased robustness, this procedure is repeated across different starting points (we use 16 equally-spaced points between $-\mathit{\lambda }/\mathrm{4}\le y<\mathit{\lambda }/\mathrm{4}$), with the background trend defined to be the average over all resulting midpoints (Fig. A1b). It follows that the fluctuations are then the difference between the original data and the trend (Fig. A1c).

Figure A1 Example demonstrating the spline-based filtering used.