Like body MHD waves able to propagate through the bulk plasma, surface waves at the interface of two media can be derived from wave solutions in displacement ξ r , t = ξ k , ω e i k ⋅ r − ω t about equilibrium (subscript 0’s) in the linearised Ideal MHD equation ρ 0 ∂ 2 ξ ∂ t 2 = ∇ ξ ⋅ ∇ p 0 + γ p 0 ∇ ⋅ ξ + 1 μ 0 ∇ × B 0 × ∇ × ξ × B 0 + 1 μ 0 ∇ × ∇ × ξ × B 0 × B 0 , which has corresponding density ρ, pressure p, and magnetic field B perturbations δ ρ = − ρ 0 ∇ ⋅ ξ − ξ ⋅ ∇ ρ 0 δ p = − γ p 0 ∇ ⋅ ξ − ξ ⋅ ∇ p 0 δ B = ∇ × ξ × B 0 , with γ being the adiabatic index.

Often surface waves are considered at discontinuities with the most studied being unbounded tangential discontinuities (TDs; Kruskal and Schwartzschild, 1954; Sen, 1963; Southwood, 1968; Goedbloed, 1971; Walker, 1981; Pu and Kivelson, 1983a; Pu and Kivelson, 1983b), pressure balanced surfaces with no threaded mass/magnetic flux—a reasonable approximation to the magnetopause in the absence of reconnection. Surface waves are, however, also supported by the other MHD discontinuities and shocks (Lubchich and Pudovkin, 1999; Lubchich and Despirak, 2005; Ruderman et al., 2018), as well as transition layers (Chen and Hasegawa, 1974; Lee and Roberts, 1986; Uberoi, 1989; De Keyser et al., 1999).

Figure 1B demonstrates the key features of a surface wave on a planar TD with uniform half-spaces. They are collective modes of vortical plasma motions, mathematically constructed from magnetosonic waves on each side independently obeying dispersion relation k n , a 2 = − k t 2 + ω a 4 v A , a 2 ω a 2 + c S , a 2 ω a 2 − k t ⋅ v A , a 2 , (1) where a represents one half-space (magnetosphere/magnetosheath), n and t denote normal and tangential directions, v _A and c _S are the Alfvén and sound speeds, and ω _a is the rest frame angular frequency. If the plasma has velocity u _0,a in the local frame, the Doppler shift gives ω = ω _a + k _t ⋅u _0,a .

Surface waves are necessarily localised to their discontinuity, hence exhibit evanescence normal to the boundary, i.e., R e ( k n , a 2 ) < 0 . Fluctuations’ decay scale depends on plasma conditions, thus can be different on either side. Amplitudes peak at the discontinuity, necessitating a reversal of polarisation across the boundary (Southwood, 1974).

The two wave solutions are tied together through boundary conditions: the tangential wave vector k _t and wave frequency ω are the same on both sides, and the normal displacement ξ _n (or equivalently velocity) and total pressure perturbation δp _tot = δp + B ₀ ⋅ δ B/μ ₀ are continuous. This leads to general surface wave dispersion relation applied to the magnetopause. k n , m s h ρ 0 , m s h ω msh 2 − k t ⋅ v A , m s h 2 = k n , m s p ρ 0 , m s p ω msp 2 − k t ⋅ v A , m s p 2 , (2) where subscripts msh and msp represent the magnetosheath and magnetosphere sides, respectively. Eq. 2 must be solved numerically and yields quasi-fast and quasi-slow modes (Pu and Kivelson, 1983a). It can be simplified assuming incompressibility (c _S → ∞), where k n 2 = − k t 2 on both sides and Eq. 2 becomes. ω = k t ⋅ ρ 0 , m s h u 0 , m s h + ρ 0 , m s p u 0 , m s p ρ 0 , m s h + ρ 0 , m s p ± ρ 0 , m s h k t ⋅ v A , m s h 2 + ρ 0 , m s p k t ⋅ v A , m s p 2 ρ 0 , m s h + ρ 0 , m s p − ρ 0 , m s h ρ 0 , m s p ρ 0 , m s h + ρ 0 , m s p 2 k t ⋅ u 0 , m s h − u 0 , m s p 2 .

This has forward and backward propagating solutions with respect to k _t , though as the flow shear increases one is reversed becoming a “negative energy” wave (e.g., Mann et al., 1999). Increasing the shear further results in exponential growth in time, corresponding to the classical criterion for the Kelvin-Helmholtz Instability (KHI, Chandrasekhar, 1961). While intimately related to surface waves, we shall not discuss KHI further here (see instead review Masson and Nykyri, 2018). At the magnetopause typically u _0,msh ≫ u _0,msp, ρ _0,msh ≫ ρ _0,msp and B _0,msh ≪ B _0,msp, hence the approximate relation ω ≈ k t ⋅ u 0 , m s h ± ω 0 ω 0 = k t ⋅ B 0 , m s p 2 + k t ⋅ B 0 , m s h 2 μ 0 ρ 0 , m s p + ρ 0 , m s h ≈ k t ⋅ B 0 , m s p μ 0 ρ 0 , m s h (3) holds for k t ⋅ u 0 , m s h 2 ρ 0 , m s p / ρ 0 , m s h ≪ ω 0 2 , consisting of a natural frequency ω ₀ for no flow shear (Chen and Hasegawa, 1974; Plaschke and Glassmeier, 2011) along with an advective Doppler shift.

In reality the magnetopause has finite thickness ∼ 100 – 2500 k m (Berchem and Russell, 1982; Paschmann et al., 2005), whereas the above theory treated it as infinitesimally thin. While this limit is valid for wavelengths much larger than the thickness, a finite-width boundary has been considered through either a continuously-varying transition (Chen and Hasegawa, 1974; De Keyser et al., 1999) or uniform layer bounded by two discontinuities (Lee et al., 1981). Finite thickness introduces inner and outer surface modes, and can allow waves to propagate inside the layer due to the presence of turning points. Surface waves may resonantly convert to Alfvén or slow magnetosonic waves if their resonance conditions, ω 2 = k t ⋅ v A 2 or ω 2 = k t ⋅ v A 2 c S 2 / v A 2 + c S 2 , become locally fulfilled within the transition. This irreversible mode conversion leads to damping of the surface mode, even in the absence of any dissipation in the system (Chen and Hasegawa, 1974; Lee and Roberts, 1986; Uberoi, 1989). Damping (or conversely growth) of surface waves results, through Eq. 1, in the evanescent magnetosonic waves exhibiting phase motion along the normal direction (Pu and Kivelson, 1983a; Archer et al., 2021).

The theory presented in Section 2.1 holds for unbounded TDs, whereas magnetospheric field lines are necessarily terminated at their intersection with the ionosphere, as displayed in Figures 2A, B for simple box and cylindrical magnetospheric models (e.g., Southwood, 1974; Kivelson et al., 1984). The ionospheric boundary conditions are highly reflecting for magnetosonic modes, even more than for Alfvén waves (Kivelson and Southwood, 1988), meaning standing surface waves might form between conjugate ionospheres, known as surface eigenmodes (Chen and Hasegawa, 1974; Plaschke and Glassmeier, 2011). These modes have been considered on either the magnetopause or plasmapause, with quantization condition k t ⋅ B 0 , m s p B 0 , m s p = ± j π S , (4) where j ∈ N and S is the field line length. Surface eigenmodes have conceptual similarities to field line resonances, Alfvén waves locally standing on field lines due to ionospheric reflection (e.g., Southwood, 1974). They are, however, rather different from cavity (for a closed magnetosphere with quantized azimuthal wavenumbers; Kivelson et al., 1984; Kivelson and Southwood, 1985) and waveguide (for an open-ended magnetotail with a spectrum of azimuthal wavenumbers; Samson et al., 1992; Wright, 1994) modes. These eigenmodes instead consist of propagating, rather than evanescent, magnetosonic waves that form approximately radially standing structure from reflection by boundaries or turning points. See Archer et al. (2022) for further comparison.

Magnetopause surface eigenmodes (MSEs) are typically considered around the subsolar meridian, where flow shears are low. Its eigenfrequency in an incompressible box model is thus given by Eq. 3, making it the lowest frequency magnetospheric normal mode and highly penetrating into the magnetosphere due to expected low azimuthal wavenumbers (Plaschke et al., 2009b; Plaschke and Glassmeier, 2011; Archer and Plaschke, 2015). The frequency has an approximately linear dependence on solar wind speed—via balance of magnetospheric magnetic and solar wind dynamic pressures, proportionality of magnetosheath and solar wind densities, and “stiffness” of the magnetosphere (Chen and Hasegawa, 1974; Archer et al., 2013a; Archer and Plaschke, 2015; Nenovski, 2021). Compressibility should modify this only slightly (Pu and Kivelson, 1983a).

It had been suggested fast solar wind might inhibit MSE, as meridional magnetosheath flow could reverse one of the counter-propagating surface waves (Plaschke and Glassmeier, 2011). Later this was considered important only for large dipole tilts, based on time-of-flight calculations within empirical models (Archer and Plaschke, 2015). Away from the subsolar meridian it was thought MSE would be advected tailward by the magnetosheath. However, it was shown theoretically (as well as in observations and simulations) that MSE can stand stationary against the flow across a wide local time range on the dayside, trapping wave energy locally (Archer et al., 2021).

In box models MSE currents flow only within the magnetopause, forming field-aligned currents at low altitudes which close via ionospheric Pedersen currents (Plaschke and Glassmeier, 2011). Whether surface modes directly have significant effects on the ionosphere or ground, and where these might map to, has been debated (Kivelson and Southwood, 1988; Southwood and Kivelson, 1990; Southwood and Kivelson, 1991; Kozyreva et al., 2019; Archer et al., 2023a).

Numerous fundamental conceptual challenges concerning magnetopause surface waves remain, even in the linear theory. These are due to standard box/cylindrical models (Figures 2A, B) being oversimplifications of the magnetospheric environment, neglecting aspects of the full physics.

Standard models feature straight equilibrium field lines, however, field line curvature significantly affects MHD waves. To account for this one requires a magnetic (field-aligned) coordinate system with corresponding metric tensor (Stern, 1970; Stern, 1976; D’haeseleer et al., 1991). Note an orthogonal system may not exist, e.g., in the case of background field-aligned currents (Salat and Tataronis, 2000; Rankin et al., 2006). Such methods applied to the magnetospheric Alfvén continuum have shown eigenfrequencies differ from time-of-flight estimates typically by ∼ 20 – 75 % and vary with polarisation (Singer et al., 1981; Rankin et al., 2006; Elsden and Wright, 2020). The problem is more complex for surface waves, since solutions in each half-space, with likely different magnetic coordinate systems, must be tied together. Nonetheless, a simple model to determine surface waves’ sensitivity to field line curvature is the hydromagnetic wedge (Radoski, 1970) depicted in Figure 2C. Here cylindrical coordinates describe axial field lines, each with a constant radius of curvature (their radial coordinate), confined between two angles of azimuth corresponding to ionospheric boundaries.

Field lines are terminated in the ionosphere on both sides of the discontinuity in standard models. While this is valid for the plasmapause, at the magnetopause field lines in the magnetosheath should be open (Kozyreva et al., 2019). A modified box model with open magnetosheath flux is illustrated in Figure 2D. While the quantization condition on closed magnetospheric field lines (Eq. 4) should be unaffected, continuity of normal displacement across the boundary (Walker, 1981; Plaschke and Glassmeier, 2011) requires zero perturbation above/below the intersection of the magnetopause with the ionosphere, i.e., ξ x x = 0 , y , z , t = r e c t π z / S cos k z z exp i k y y − ω t . The Fourier decomposition of this boundary condition introduces additional magnetosheath wavenumbers. These must also follow the magnetosonic dispersion relation (Eq. 1), hence may consist of propagating components in addition to evanescent waves. The overall effect would be a form of diffraction into the magnetosheath.

The magnetic field models presented have still been highly simplified. Close to Earth the field is reasonably approximated by a dipole, but differs substantially from this in the outer magnetosphere. Dipole equilibrium magnetic field models have been used to explore surface modes due to velocity shears (Leonovich and Kozlov, 2019), though by definition these models cannot include background currents or thermal pressure gradients. Placing an image dipole in the solar wind gives an analytic closed magnetosphere with planar infinitely conducting magnetopause and two magnetic null points at the cusps, as shown in Figure 2E (Chapman and Bartels, 1940). More representative closed magnetosphere models, like in Figure 2G, can be constructed by perturbing a dipole. Introducing (just two) spherical harmonic corrections to dipole Euler potentials produces a reasonable magnetosphere (Stern, 1967), though how to describe magnetosheath field lines in this framework is unclear. Alternatively, expressing the field as a scalar potential in parabolic harmonics with contributions from the dipole and magnetopause currents allows one to confine the geomagnetic field within a paraboloidal magnetopause (Stern, 1985), though magnetic coordinates must be determined numerically (Degeling et al., 2010).

In closed magnetosphere models the magnetopause is a TD. While it is often stated the magnetopause may be treated as a rotational discontinuity for an open magnetosphere (e.g., Sonnerup and Ledley, 1974), this neglects any density/pressure gradients present between the two media. Therefore, for an open magnetosphere the magnetopause must consist of both compressional and rotational boundaries (Dorville et al., 2014). The image dipole model (Figure 2E) can be simply extended to represent an open magnetosphere, as shown in Figure 2F (Kan and Akasofu, 1974). Constructing a realistic open magnetosphere model analytically remains an outstanding challenge (Zaharia and Birn, 2005).

These more representative models (Figures 2E–G) may help understanding effects of the polar cusps. Alfvén wave propagation is significantly affected by local variations in magnetic field strength and/or curvature when wavelengths are comparable to inhomogeneity scales (Pilipenko et al., 1999; Pilipenko et al., 2005), leading to reflection up to ∼ 80 – 90 % . If similar holds for surface waves, surface eigenmodes might stand between reflection points in conjugate cusps rather than ionospheres. Another benefit to these models would be in probing non-resonant wave coupling. While Alfvén and magnetosonic modes are independent in uniform media, when inhomogeneities are introduced waves necessarily have mixed properties (Radoski, 1971; Goossens et al., 2019). Understanding the partial, irreversible conversion of surface waves’ compressional energy into Alfvén waves could help determine potential impacts of surface waves on the system and any filtering/processing the magnetosphere imposes upon them (cf. Pilipenko et al., 1999).

Thus far surface waves on a single boundary have been considered. However, it is clear from Figure 1A numerous boundaries exist within the magnetosphere. Given the large-scale nature of surface eigenmodes across the magnetic field, they cannot exist in isolation. This motivates a multi-boundary approach, such as in Figure 2H. Some progress towards this has considered eigenmodes of the outer magnetosphere, modeled as a slab or annular cylinder bounded by magnetopause and plasmapause (Nenovski et al., 2007; Nenovski, 2021). How the surface eigenmodes of individual boundaries are modified by, and couple to, the existence of another boundary has yet to be explored. Furthermore, introduction of the bow shock and equatorial ionospheric boundaries is required for a more complete description.

Resonant absorption for a finite-width boundary has been poorly explored outside of standard models. Where (either inside or outside) the transition layer mode conversion occurs and how this varies with surface mode harmonics is not well understood in more realistic setups. Since mode conversion provides a means for field-aligned current generation from a purely compressional wave (Southwood and Kivelson, 1990; Itonaga et al., 2000; Plaschke and Glassmeier, 2011) it is an important topic relevant to the ionospheric and ground-based impacts of surface modes (Kivelson and Southwood, 1988; Archer et al., 2023a).

Linear surface wave theory typically solves either initial value (where the boundary is perturbed and allowed to evolve) or eigenvalue problems (where normal modes are sought). Under these approaches it is not possible to self consistently consider variable solar wind forcing. However, solar wind and magnetosheath plasmas are highly turbulent, modifying conditions present adjacent to the magnetopause over timescales comparable to (or shorter than) surface wave periodicities. It has been suggested, by analogy with a driven harmonic oscillator with stochastically varying eigenfrequency, magnetosheath turbulence might suppress the surface mode (Pilipenko et al., 2017). The potential damping factor in such a toy model using representative turbulent amplitudes/spectra have not yet been estimated. Expanding this approach to self-consistently treat turbulent driving in addition to its effect on the eigenfrequency is required to more realistically formulate this problem.

It is finally worth noting only linear wave theory has been discussed. Nonlinear effects have been investigated analytically on planar TDs (cf. Figure 1B) for incompressible plasmas (Hollweg, 1987; Alì and Hunter, 2003). While this is of greatest relevance to a KH-unstable boundary, even in the marginally stable and weakly nonlinear regimes surface waves undergo wave steepening, crest/trough sharpening, and non-local self-interactions, leading to their breakdown in finite time. There is certainly more scope to explore this topic.

Many of these current challenges simply cannot be realistically treated by analytic theory. Nonetheless, modifications to simple models will likely provide insight into the physics. However, to more representatively model them in a complex environment like the magnetosphere necessitates numerical simulations.