Techniques to infer mass density structures are well established in terrestrial and solar seismology. Our idea is to use similar techniques to infer the mass distribution in the terrestrial magnetosphere. The basic principle of magnetoseismology is that ρ is related to the frequency and mode structure of standing Alfvén waves. Figure 1 illustrates this relationship using a string model of magnetic field lines. The frequency of vibrations of the string (blue curves) is determined by the tension (restoring force) of the string and the mass (filled red circles) attached to the string. The discrete mass distribution is only for illustrative purposes. In reality, the mass is distributed continuously.

Figure 1A illustrates the fundamental and second harmonic modes for the case of a uniform mass distribution along the field line. The mode structure is a sine function for either harmonic, and the frequencies are related as f _a2 = 2f _a1. The time-of-flight calculation described below gives exact solutions of the mode frequency and structure for all harmonics.

Figure 1B illustrates the case of a nonuniform mass density distribution with a peak at the equator. In this case, the mode structure of the fundamental mode is not a sine function, and the frequency of the fundamental mode (f _b1) is lower than f _a1. However, the equatorial mass density does not affect the mode structure or the frequency of the second harmonic because the string displacement is zero (node) at the equator. As a consequence, the mode of the second harmonic is a sine function, the same as in Figure 1A. The mass density effects on the wave properties occur for higher harmonics also. This means that we can infer the mass density distribution from the frequencies and mode structures determined using spacecraft data.

To advance the concept illustrated in Figure 1 to magnetoseismology of the real magnetosphere, we obtain the relationship between the waves and ρ using the cold plasma MHD wave equation (e.g., Radoski and Carovillano, 1966) ∂ 2 δ E ∂ t 2 = V A × V A × ( ∇ × ∇ × δ E ) , (1) where δE is the wave electric field and V _A is the Alfvén velocity, which depends on the magnetic field B and ρ as V _A = B ( μ 0 ρ ) − 1 ; μ ₀ is the permeability of free space. A justification for using the cold plasma equation is found in the appendix of Singer et al. (1981). When perfect wave reflection is imposed at the ionosphere, the equation for a dipole magnetosphere yields two purely transverse standing Alfvén wave solutions, which are known as the axisymmetric toroidal mode (Radoski and Carovillano, 1966) and the guided poloidal mode (Radoski, 1967). Magnetic field perturbation and fluid motion are in the azimuthal direction for the toroidal mode and in the meridional plane for the poloidal mode. The toroidal (poloidal) mode corresponds to the limit of m = 0 (∣m∣ = ∞), where m is the azimuthal wave number (<0 for westward propagation). The polarization state is relevant to magnetoseismology because the wave frequency depends on it. In the dipole field, the fundamental toroidal frequency (f _T1) is 1.4 times the fundamental poloidal frequency (f _P1) if ρ is constant along the field line (Cummings et al., 1969). This translates to a factor of ∼2 difference in ρ (see Eq. 4). The idealized wave modes explain many observable features of standing Alfvén waves in the magnetosphere, both waves with toroidally (azimuthally) oscillating magnetic perturbations excited by large-scale solar wind disturbances and waves with poloidally (radially) oscillating magnetic perturbations excited by internal instabilities. Spacecraft detect toroidal waves routinely, and magnetoseismic studies usually rely on these waves.

Poloidal waves are detected also, but these waves have not been used much in magnetoseismology. Considering the fact that poloidal waves are not detected by ground magnetometers because of the ionospheric screening of high-m waves (Hughes and Southwood, 1976) and also the fact that the waves exhibit a different local time occurrence distribution than toroidal waves (Arthur and McPherron, 1981), poloidal waves could be a valuable resource in future magnetoseismic studies. Poloidal waves have the advantage of being excited nearly exclusively at the second harmonic (Cummings et al., 1969; Takahashi and McPherron, 1984; Anderson et al., 1990; Liu et al., 2013), reducing the uncertainty in harmonic mode identification. But a disadvantage to using poloidal waves is that the wave frequency depends on the radial pressure gradient, which is not always well determined (Denton, 2003, and references therein).

The magnetospheric magnetic field significantly differs from the dipole field at large distances or during geomagnetically disturbed periods, making it difficult to exactly solve Eq. 1. Fortunately, MHD-scale disturbances quickly settle to toroidal eigenmode oscillations (Allan et al., 1986; Lee and Lysak, 1989) through the FLR process (e.g., Chen and Hasegawa, 1974), and we can treat each field line to be an independent oscillator in the context of magnetoseismology. For example, we can use the time-of-flight approximation for the fundamental frequency f ₁ on a field line (Warner and Orr, 1979; Wild et al., 2005) f 1 = 2 ∫ S N d s V A − 1 , (2) where s is distance along the field line and the integral is taken between the southern (S) and northern (N) ionospheric footpoints. In this approximation, there is no distinction between toroidal and poloidal frequencies, the frequency of the nth harmonic (f _n ) is equal to nf ₁, and the f _n value is higher than that obtained by more accurate methods (Takahashi and McPherron, 1984). More accurate eigenmode solutions are obtained by numerically solving the equation μ 0 ρ ∂ 2 ( ξ i / h i ) ∂ t 2 = 1 h i 2 B ⋅ ∇ { h i 2 [ B ⋅ ∇ ( ξ i / h i ) ] } , (3) where ξ _i is the field line displacement associated with the wave and h _i is the scale factor vector of the background magnetic field B, with the suffix i indicating the direction within the plane perpendicular to B (Singer et al., 1981). To solve this equation for a general magnetic field geometry, one selects two adjacent model magnetic field lines to specify the direction (polarization axis) of magnetic field perturbation. This flexibility is valuable in consideration of theoretical studies (Lee et al., 2000; Wright and Elsden, 2016) indicating that the polarization axis of toroidal waves is not tangential to the magnetic field L shells when the ρ distribution is not axisymmetric. The two field lines are best chosen at the magnetic equator, where the properties have the strongest effect on the wave frequency. A somewhat more self-consistent approach would be to use the equations of Rankin et al. (2006), who solve for the coupled toroidal and poloidal waves.

A practical procedure to estimate the mass density (denoted ρ _est) from the observed wave frequency f _obs is as follows. In the first step, we obtain the reference eigenfrequency f _ref by solving the wave equation for a reference mass density ρ _ref (e.g., 1 amu cm⁻³) at a reference point (e.g., magnetic equator) after choosing models for the magnetic field and mass density variation along the field line. In the second step, we obtain ρ _est using the relationship ρ est = ρ ref ( f ref / f obs ) 2 . (4)

The mass density values at other locations along the field line are obtained using the model field line mass distribution function.

The quality of the models for the magnetic field and field line mass density distribution determines the accuracy of ρ _est. For the magnetic field, several models are available (e.g., Tsyganenko, 1989; Tsyganenko and Sitnov, 2005; Sitnov et al., 2008), and it is even possible to use magnetic fields obtained by global MHD simulation (e.g., Claudepierre et al., 2010). We can choose the best field model by comparing model fields with the field that is observed by the same satellite used for wave observation. This is an advantage of using spacecraft data instead of ground magnetometer data. An even greater advantage relates to determination of the equatorial location of the field line, and hence the L shell and magnetic local time (MLT). That is much less of a problem for spacecraft data, especially for spacecraft with near-equatorial orbits. For field lines mapping from the ground to geostationary orbit (L ∼ 7) or beyond, inaccuracies in mapping can be severe, where L is the magnetic shell parameter.

As for the mass density field line distribution model, we cannot impose many constraining conditions because we have a small number of observable eigenmodes, unlike in terrestrial or solar magnetoseismology. Therefore, we adopt models with a small number of free parameters. The most frequently used mass density model has only two free parameters (ρ _eq and α) ρ = ρ eq ( L R E / r ) α , (5) where ρ _eq is the equatorial mass density, R _E is the Earth’s radius, r is geocentric distance to the field line, and the power law index α specifies how ρ varies along the field line (Radoski and Carovillano, 1966; Cummings et al., 1969). We can add more flexibility to the model density by using polynomial expansion in terms of a parameter related to s (Denton et al., 2001, 2004; Takahashi and Denton, 2007). In the Takahashi and Denton (2007) study, the parameter is defined to be τ = ∫ d s / B , (6) where B is the magnitude of the magnetic field. The mass density model is expressed as log 10 ρ = c 2 τ 2 + c 4 τ 4 + c 6 τ 6 . (7)

Only even terms appear in this equation because of the assumption that the mass density distribution is symmetric about the equator. Although this model has only three free parameters, it is capable of modeling an equatorial enhancement of mass density in a way that the power-law model (Eq. 5) cannot.

Any magnetospheric spacecraft carrying science experiments is a good data source for magnetoseismology. Three types of spacecraft data have been used. They are E and B fields and the plasma bulk velocity (V). Convective anisotropy of energetic particles can be used as a proxy to V (Takahashi et al., 2002). Figure 2 shows examples of toroidal waves detected by two representative magnetospheric spacecraft with low orbital inclination: Van Allen Probes (Radiation Belt Storm Probes, RBSP)-B and Time History of Events and Macroscale Interactions during Substorms (THEMIS)-D. Figure 2A shows the equatorial projection of the selected orbits of these spacecraft in solar magnetic (SM) coordinates. RBSP-B, with apogee at ∼6 R _E, covers the inner magnetosphere. THEMIS-D, with apogee at ∼12 R _E, covers the outer magnetosphere. Figure 2B shows the coordinate systems used for spacecraft position and measured vector quantities.

Figures 2C,D were generated from the toroidal components, δE _ν and δB _ϕ, measured by RBSP-B on the orbit shown in Figure 2A. The cross-power spectra (Figure 2C) show several toroidal harmonics, the most obvious being the fundamental (T1), second (T2), and third (T3) harmonics. The cross-phase spectra (Figure 2D), displayed only when the δE _ν-δB _ϕ coherence is higher than 0.5, also show several bands corresponding to the band structure in the cross-power spectra. The cross-phase spectra are the key to identifying the harmonic modes when many harmonics coexist or when the spectral intensity changes with time (Takahashi and Denton, 2021).

Figure 2E was generated from the δV _ϕ component measured by THEMIS-D on the orbit shown in Figure 2A. For an equatorially orbiting spacecraft, the velocity is a sensitive indicator of odd mode waves, which have antinodes at or near the equator. The velocity often is the best quantity for toroidal wave analysis when the electric fields measured on the same spacecraft are contaminated by spacecraft wake or charging. In this example, a strong T1 spectral line is visible. A caveat with the velocity data is that the δ V spectral intensity diminishes at L < 7, where the δ V is weak because of strong background B.

Included in Figure 2A are the orbit of geostationary satellites (e.g., Geostationary Operational Environmental Satellite, GOES) and an orbit of Active Magnetospheric Particle Tracer Explorers (AMPTE)/Charge Composition Explorer (CCE). The GOES satellites provide continuous B field data (not shown) at L ∼ 7. Harmonic mode identification is relatively easy with the GOES data because the magnetic latitude (MLAT) of the spacecraft does not change. Denton et al. (2016) used 12 years of data from multiple GOES satellites to develop a number of models of varying complexity for ρ at geostationary orbit. The most complicated models could determine ρ within a factor of 1.6, accounting for about two-thirds of the variance. Some GOES spacecraft carry detectors for energetic (>80 keV) protons (Rodriguez, 2014), and data from the detectors can be used to determine the frequency of oscillatory convective anisotropy induced by standing Alfvén waves (e.g., Takahashi et al., 2002). This capability remains to be utilized.

CCE had an orbital configuration intermediate between THEMIS and RBSP and was operational between 1984 and 1989. Min et al. (2013) used magnetometer data from this spacecraft to construct mass density models covering L = 4–9 and MLT = 0300–1900. The study also determined toroidal wave frequencies using GOES magnetometer data and found the frequencies to be consistent with those at CCE for L ∼ 7.

Magnetoseismology works when the driver fast mode waves for exciting toroidal waves have a wide spectral band to excite multiharmonic toroidal waves over a wide range of L. If the fast mode waves have a narrow spectrum, toroidal waves will be excited in a narrow L range and we will not be able to determine the L profile of ρ. Monochromatic fast mode waves such as wave guide modes could be excited in the magnetosphere and could produce ground magnetic pulsations with L-independent frequencies (Marin et al., 2014) while exciting magnetospheric toroidal waves on an isolated L shell. However, as Figure 2 indicates, toroidal waves (especially in the dayside magnetosphere) are usually excited over a wide frequency range in response to broadband fast mode waves generated either by dynamic pressure variations intrinsic to the solar wind (Sarris et al., 2010) or by compressional ULF waves generated in the ion foreshock (Takahashi et al., 1984). We believe that broadband fast mode waves are always present in the magnetosphere in addition to possible waveguide modes.

A number of studies used toroidal wave frequencies (f _Tn , n being the harmonic mode number) to find an optimal value of the α parameter appearing in Eq. 5. These studies found α values in the range 0–2, which is closer to α = 0–1 for the electron diffusive equilibrium expected in the plasmasphere rather than a collisionless distribution (α = 3–4) expected in the plasmatrough (Angerami and Carpenter, 1966; Takahashi et al., 2004). For example, Takahashi et al. (2004) obtained α ∼ 0.5 by a statistical analysis of the f _Tn /f _T1 ratio of toroidal waves detected by the Combined Release and Radiation Effects Satellite (CRRES) spacecraft at L = 4–6 in the postnoon sector. Takahashi et al. (2015a) obtained α ∼ 0 from a detailed analysis of multiharmonic toroidal waves (n = 1–11) detected by the RBSP spacecraft during a plasmaspheric pass in the dawn sector. A statistical analysis of the f _Tn /f _T3 ratios at RBSP in the noon sector (Takahashi and Denton, 2021) found α ∼ 2 at L = 4–6 in both the plasmasphere and the plasmatrough. Note that the α value does not need be the same between the electron density (n _e) and ρ because multiple ion species with different masses and charge states, which in general have different pitch angle distributions, contribute to the latter.

In magnetoseismology, multiharmonic toroidal waves are interpreted to be superposition of independent linear waves. If the waveform is nonlinearly distorted, it will lead to regularly spaced spectral peaks and will affect α estimation. It is known that nonlinearly distorted poloidal waves produce regularly spaced spectral peaks (Higuchi et al., 1986; Takahashi et al., 2011). It is not clear whether similar distortion occurs during toroidal wave events. However, statistically determined frequency spacing between toroidal harmonics is not even, and we believe that the distortion is rare. Note that the theoretical frequencies of linear toroidal waves are evenly spaced in a dipole magnetosphere if we set α = 6 in Eq. 5 (Cummings et al., 1969; Schulz, 1996). The statistical results favoring α < 2 are an indication that nonlinear distortion is negligible.

A statistical analysis of GOES magnetometer data (Takahashi and Denton, 2007) determined the f _Tn /f _T3 ratios at geostationary orbit (L ∼ 7) for n = 1–5 as shown in Figures 3A,B. This analysis indicated that the power-law model is only a rough approximation and that the ratios change with MLT. This finding led the authors to adopt the model given by Eq. 7. The results (Figure 3C) indicate that the mass density is peaked at the equator with the peak more pronounced at the later local times. The cause of the peak remains to be determined.

A follow-up study (Denton et al., 2015) using the same data as those of Denton et al. (2016) developed a model for the α index at geostationary orbit, α 3 Dmodel = 2.2 + 1.3 ⋅ cos M L T ⋅ 15 ° + 0.0026 ⋅ A E ⋅ cos ( M L T − 0.8 ) ⋅ 15 ° + 2.1 ⋅ 1 0 − 5 ⋅ A E ⋅ F 10.7 − 0.010 ⋅ F 10.7 , (8) where F _10.7 is the solar extreme ultraviolet (EUV) flux index, AE is the auroral electrojet index, and MLT is in hours. Eq. 8 modeled binned values of α within a standard deviation of 0.3.

A recent study used observationally determined MLAT of the nodes of toroidal waves to select α values (Takahashi and Denton, 2021). The results shown in Figure 4 were obtained by a statistical analysis of the MLAT dependence of the amplitude and the phase of the δE _ν and δB _ϕ components measured by RBSP over a 6-month period during which the spacecraft were located on the dayside. The panels in the top and middle rows show the results for T3–T5 waves. The panels at the bottom indicate the relationship between α and the node latitudes assuming a dipole field, and the intersects of the vertical dashed lines (observed node latitudes) and the theoretical curves give the α values. This analysis indicated α ∼ 1.7 (horizontal dashed line), not far from α ∼ 2 derived in the same study using the frequencies.

A useful variable in magnetoseismology is the average ion mass M given by M = ρ / n e , (9) where ρ is derived from toroidal wave frequencies and n _e is sometimes available from examination of plasma wave spectra observed at spacecraft such as CRRES (Anderson et al., 1992) and RBSP (Kurth et al., 2015). The magnetospheric plasma is mostly composed of H⁺, He⁺, and O⁺, which means that the mass density is expressed as ρ ∼ n H + m H + + n He + m He + + n O + m O + , (10) with the constraint of charge neutrality n e ∼ n H + + n He + + n O + , (11) where n _i and m _i (i = H⁺, He⁺, or O⁺) are the number density and the mass of the ion species, respectively. Although we cannot determine all of the three ion number densities from the two variables ρ and n _e, we can use M to infer the ion composition. The value of M should be between 1 amu (all-H⁺ plasma) and 16 amu (all-O⁺ plasma). The n He + / n H + ratio being relatively stable (Craven et al., 1997; Krall et al., 2008), M is a good indicator of the variability of n O + .

Figure 5 shows the statistical properties of M samples derived using CRRES data for 1991 (solar maximum) and setting α = 0.5 (Takahashi et al., 2006). The occurrence distribution (Figure 5A) is mostly confined within 1–16 amu as expected with a median value of 3 amu. Because He⁺ cannot raise M to >4 amu and H⁺ carries the highest number density in general, it is concluded that there were substantial amounts of O⁺. In addition, M differs between the plasmasphere and the plasmatrough. Figure 5B shows that M ∼ 1.5 in the plasmasphere (n _e > 100 cm⁻³) and M ∼ 3 amu in the plasmatrough (n _e < 20 cm⁻³). Geomagnetic activity also controls M, as shown in Figure 5C, with higher values occurring when the ring current index Dst has larger magnitudes. O⁺ ions originate from the ionosphere, and the solar EUV intensity (F _10.7) controls the density, temperature, and scale height of the O⁺ ions that are transported to the magnetosphere. A study that combined ρ determined using toroidal wave frequencies and ions detected by particle experiments at geosynchronous orbit (Denton et al., 2011) showed that while ρ has maximum value at solar maximum, the electron density has minimum value, so that the n O + / n e ratio varied between ∼0.2 at solar maximum and ∼ 2 × 10^–3 at solar minimum.

This n O + variability leads to a solar cycle variation of ρ, as demonstrated in the statistical result (Takahashi et al., 2010) shown in Figure 6. There is a high degree of anticorrelation between F _10.7 and the T3 wave frequency (f _T3) at GOES, which means a positive correlation between F _10.7 (solar flux units, 10^{–22 m−2} Hz⁻¹) and ρ (amu cm⁻³) expressed as, log ⁡ ρ = 0.421 + 0.00390 F 10.7 . (12)

This equation indicates a factor of ∼4 variation of ρ over a solar cycle. Note that the GOES measurements were made at L ∼ 7 and the ρ samples were taken from the 0600–1200 MLT sector, which means that the spacecraft was mostly in the plasmatrough. A similar study using Geotail data (Takahashi et al., 2014) found that ρ at L ∼ 11 in the 0400–0800 MLT sector varied by a smaller factor of ∼2 over a solar cycle. This difference could be accounted for by the L localization of an O⁺-rich region, as described next.

Spatial localization of heavy ion concentration is one of the important magnetospheric phenomena that magnetoseismology can uniquely address. We show two examples. The first example (Figure 7) is taken from Takahashi et al. (2008) and shows the L profile of magnetoseismic variables for a drainage plume crossing by the CRRES spacecraft. Although the n _e profile (Figure 7B) clearly indicates the distinction between the plasmatrough and the drainage plume, there is no change in f _T1 (Figure 7A) at the trough-plume boundary. This difference is explained by a higher heavy ion concentration in the plasmatrough, which is evident in the L profiles of M (Figure 7C) and n O + / n e (Figure 7D). Figure 7E shows that O⁺ ions account for the majority of the mass density in the plasmatrough.

The second example (Nosé et al., 2015) is shown in Figure 8. Figure 8A shows that M is elevated as high as 8 amu over an L distance ∼1 just outside the electron plasmapause located at L ∼ 3. Another study, which combined observations of Arase and RBSP (Nosé et al., 2020), demonstrated that O⁺ enhancement is limited in local time as well.

The lower two panels of Figure 8 are included to show the difficulty in determining ρ using particle data. Figure 8B shows ion number densities calculated using data from the RBSP Helium, Oxygen, Proton, and Electron (HOPE) mass spectrometer (Funsten et al., 2013) in the energy range 30 eV—1 keV. Although HOPE has a lower energy limit at 1 eV, energies lower than 30 eV were excluded to avoid spacecraft charging effects. The HOPE-derived densities (<1 cm⁻³) are well below n _e (>50 cm⁻³) determined from the plasma wave spectra, meaning that the bulk of the mass density is carried by ions with energies lower than 30 eV. Figure 8C shows that there is no indication of heavy ion enhancement in the HOPE data at the location where M is elevated. This example demonstrates that magnetoseismology captures low-energy ions that contribute to the mass density but are not measured by particle instruments.

An exception to the limitation of particle experiments occurs when the cold ion population is embedded in a fast bulk flow so that the population can be detected by particle instruments with the lower energy cutoff well above the thermal energy of the cold ions. Such flow events occur during Pc5 wave events in the outer magnetosphere (Chen, 2004; Hirahara et al., 2004; Lee and Angelopoulos, 2014) and provide opportunities to validate results from magnetoseismology.

Figure 9 shows a comparison of the ion mass density derived from ion flux measurements during a Pc5 wave event reported by Hirahara et al. (2004) and the ρ value estimated from the frequency of the wave (Takahashi et al., 2014). A 12-s snapshot of the ion phase space density (Figure 9B) exhibits three peaks at negative velocities. These peaks are attributed to the cold O⁺, He⁺, and, H⁺ ions that are convected at the same δ E × B velocity. The velocity (>100 km/s peak to peak, Figure 9A) is higher than the background plasma convection velocity. The peaks are separated because the ion instrument (an electrostatic analyzer) does not distinguish ion species and the velocity is calculated assuming all detected ions are protons. By assigning a correct mass value to the ions contributing to each peak, it is possible to obtain the number density and mass density for each ion at each instrument duty cycle, as shown in Figure 9C. The mass density summed over the three ion species is 3.9 amu cm⁻³ when averaged over the time interval shown in Figure 9C. This is close to the value 3.1 amu cm⁻³ that is obtained from the toroidal Pc5 wave frequency. This comparison could be extended to a statistical study using many flow events detected by Geotail (e.g., Hirahara et al., 2004) and THEMIS (e.g., Lee and Angelopoulos, 2014).

A major goal of magnetoseismology is to develop a global model of ρ. Ideally, the model will reach a degree of maturity similar to that of existing models of the magnetic field (e.g., Sitnov et al., 2008), the electron density (e.g., Carpenter and Anderson, 1992; O’Brien and Moldwin, 2003; Archer et al., 2015; Liu et al., 2015), the He⁺ density (e.g., Gallagher et al., 2021), and the density of low energy (but excluding cold) ions (e.g., Kistler and Mouikis, 2016). Magnetoseismic studies using ground magnetometer data have made significant progress in this regard. For example, Del Corpo et al. (2020) generated a global equatorial ρ model covering L = 2.3–8 and MLT = 0600–1800 using measurements by ∼20 pairs of stations included in the European quasi-Meridional Magnetometer Array (EMMA) magnetometer network.

By contrast, there is much room for improvement in spacecraft data analysis. Spacecraft studies are invaluable because they provide information on the configuration of the background magnetic field and on n _e (for derivation of M), as we stated in sections 2.3 and 3.3. Although multiyear spacecraft observations cover the entire dayside magnetosphere (see Figure 2) as well as a large portion of the magnetotail, statistical analysis of the data has been limited. Notable exceptions are GOES (e.g., Takahashi et al., 2010) and Geotail (Takahashi et al., 2014) studies covering a solar cycle and an AMPTE/CCE study covering ∼4 years (Takahashi et al., 2002; Min et al., 2013).

Figure 10 illustrates the potential of spacecraft data for the global model. Figure 10A shows the rate of detection of fundamental toroidal waves obtained in a study (Takahashi et al., 2015b) that used ion bulk velocity but did not convert the wave frequency to ρ. The waves are detected at a high rate on the dayside from L = 6 to L = 12 (spacecraft apogee). The rate becomes low on the nightside, but toroidal waves are still detected there along with Pi2 pulsations, most often after substorm onsets (Takahashi et al., 1988; Takahashi et al., 2018) and when ULF waves generated in the ion foreshock penetrate deep into the magnetosphere (Takahashi et al., 2020). The presence of substorm-related toroidal waves is evident in Figure 10B as a region of large δV _ϕ amplitudes in the premidnight sector.

Figure 11 shows a magnetoseismic study using toroidal waves detected by the Arase spacecraft in the midnight sector away from the magnetic equator. The waves were detected after Pi2 onsets on the ground. An RBSP spacecraft located near the magnetic equator detected compressional oscillations, which can be cavity mode oscillations. Because the AE index had moderate values (<200 nT) during the wave event in this example, we expect that nightside Pi2 waves and toroidal waves are commonly excited and can be easily detected off the magnetic equator in association with small substorms or other minor disturbances in the magnetotail. Because it appears difficult to determine nightside toroidal wave frequencies with ground magnetometers (Takahashi et al., 2020), we expect that a model derived using spacecraft data will perform better on the nightside than models derived using only ground data.

To relate observed f _Tn to ρ, standing wave equations (e.g., Eq. 3) are solved usually assuming perfect reflection, corresponding to infinitely high height-integrated Pedersen conductivity Σ_P, at a fixed ionospheric altitude. In reality, the ionosphere has a finite thickness and the conductivity is finite. We discuss whether the assumption is appropriate.

The assumption of a thin ionosphere is justified because the thickness of the ionosphere (∼300 km) is much shorter than the hemispheric length (>15,000 km) of magnetic field lines at L > 2.5, where reliable measurements of wave frequency can be made by spacecraft on low-inclination elliptical orbits (Takahashi and Denton, 2021). On these field lines, the Alfvén wave velocity at the ionospheric altitude is usually higher than near the equator, making the Alfvén wave travel time through the ionosphere much smaller than that through the region above the ionosphere. This means that the details of wave propagation through the ionosphere do not affect f _Tn in any significant way.

There are questions about the Σ_P. This conductivity, which controls the damping rate of toroidal waves, depends on solar illumination and particle precipitation from the magnetosphere, both of which are a function of latitude and local time as well. According to numerical studies of the Σ_P dependence of f _Tn (e.g., Newton et al., 1978), the frequency is very close to that of perfectly reflected waves when Σ_P is higher than a critical value (denoted Σ_P0) corresponding to impedance matching between the ionosphere and the waves. At locations where solar illumination is low or zero, the conductivity may become lower than Σ_P0, leading to strong damping of the waves or transition of the waves to free-end modes with lower f _Tn values (Newton et al., 1978). If one end of a field line is anchored to the sunlit part of the ionosphere and the other is anchored to the dark part, theory predicts that the usual half-wave T1 modes turn to quarter-wave modes (Allan and Knox, 1979). Quarter-wave modes at L ∼ 3 have been detected at the dawn terminator by ground magnetometers (Obana et al., 2008). If the quarter-wave and half-wave modes are not distinguished, it will lead to a serious error in ρ. Investigation of the quarter-wave modes in space remains to be done.

Yet nightside multiharmonic toroidal waves are readily detected by spacecraft and exhibit properties consistent with high ionospheric reflection even when observed within the plasmasphere (Takahashi et al., 2020), where precipitation is not expected to be high enough to maintain high Σ_P according to empirical ionospheric density models (e.g., Wallis and Budzinski, 1981). This poses an interesting question of what elevates the ionospheric conductivity and whether high conductivity occurs commonly to make nightside magnetoseismology possible. One interesting possibility is that ULF waves themselves enhance ionospheric conductivity through modulation of electron precipitation (Jaynes et al., 2015). Wang et al. (2020) reported strong modulation of the conductivity by storm-time compressional Pc5 waves. Whether similar precipitation modulations occur during less geomagnetically active periods remains to be understood.

Magnetoseismology based on spacecraft data poses challenges that are not encountered with ground magnetometer data. First, wave frequency and amplitude seen from a moving spacecraft change continuously even if the waves do not have intrinsic temporal variations. The L dependence is particularly important. For example, f _T1 within the plasmasphere decreases from ∼20 mHz at L ∼ 2.5 to ∼ 5 mHz at L ∼ 4 (Takahashi and Anderson, 1992; Takahashi and Denton, 2021). As spacecraft such as THEMIS and RBSP move rapidly in the radial direction in the L < 4 region, it is necessary to choose a proper time window for spectral analysis so that the spatial variation of the frequency is resolved. The studies cited in this paper used data windows of a fixed length, although Min et al. (2013) introduced a variable data sampling rate to handle the spatial variation of f _Tn at L > 4.

The inferred ρ is proportional to the inverse square of the observed Alfvén frequency (Eq. 4), and we usually consider that frequency to be the largest source of error. The uncertainty of the frequency can be determined from the frequency spectrum (Denton et al., 2001). Other sources of error are the magnetic field at the spacecraft location and the field line dependence of the magnetic field and mass density. In most cases, the magnetic field at the spacecraft location is known. In the best case, the measurement is near the magnetic equator, which is the region that most greatly affects the Alfvén frequency. The field line dependence of the magnetic field and mass density does have some effect on the inferred value of the mass density, but that effect is usually smaller than the uncertainty associated with the frequency. If, on the other hand, the field line is mapped from low-Earth orbit or the ground, or if the spacecraft is in the outer magnetosphere (particularly at L greater than 8), there can be significant uncertainties for the magnetic field and/or field line mapping (Takahashi et al., 2006; Takahashi et al., 2010).

Another issue related to spacecraft motion is frequency shift. The shift occurs because toroidal waves have a finite L width (FLR width) within which the wave phase changes by 180°. In regions where the wave frequency decreases with L (e.g., the plasmasphere), the wave phase is delayed from lower L to higher L. This spatial phase structure leads to frequency downshift (upshift) at spacecraft moving to higher (lower) L by an amount given by Δ f = v L 2 π ε , (13) where v _L is the equatorial L-crossing speed and ɛ is the equatorial semiwidth of FLR (Vellante et al., 2004; Heilig et al., 2013). The frequency shift is found to be significant (20–25 mHz for T1 waves at L < 2.4) when observations are compared between the polar-orbiting Challenging Minisatellite Payload (CHAMP) spacecraft and ground magnetometers (Heilig et al., 2013).

At equatorial-orbiting spacecraft, Δf is smaller but may not be negligible. As an example, we evaluate Δf at RBSP. For the spacecraft, v _L has a peak value of ∼ 5 km/s at L = 1.5, decreases to ∼3 km/s at L = 3, and becomes zero at L ∼ 6 (apogee). If we assume ɛ = 200 km, which was the case for a 14 mHz toroidal wave at L ∼ 5 (Takahashi et al., 2015a), we get Δf ∼ 2 mHz at L = 3. This frequency shift is ∼ 10% of f _T1 (∼20 mHz) at L = 3 (Takahashi and Anderson, 1992) and translates to a ρ error of ∼ 20% (see Eq. 4). This error can explain why M derived from f _Tn observed on outbound RBSP passes is higher at L < 3 (large v _L ) than at L > 3 (small v _L ) (Vellante et al., 2021). Evaluation of ɛ at various radial distances and local times is necessary to improve our understanding of Δf.

Finally, fully automated methods are lacking to determine f _Tn in space. As a result, only a small fraction of satellite data that are potentially useful has been used in magnetoseismic studies. The procedure is the easiest for T1 waves detected in the outer magnetosphere using plasma bulk flow data (Takahashi et al., 2014; Takahashi et al.,2016). Regular oscillations in the azimuthal component of the velocity are almost always associated with T1 waves, and not much manual work is required to distinguish them from other waves. Processing magnetic field data from geostationary orbits is also relatively easy with the spacecraft staying at a fixed L leading to a stable appearance of the spectral intensity of each harmonic. Processing data from elliptically orbiting spacecraft is the most difficult because of the rapid change of spacecraft radial distance, crossing of the nodes, and the presence of waves other than toroidal waves.

It appears that at some stage we need to introduce a technique such as neural network analysis to automate the interpretation of the wave spectra. The main required capability of the technique is to reject spectral peaks that do not result from toroidal waves. In this regard, we note that the quality of E and B data from spacecraft depends on the mode of sensor operation, location and attitude of the spacecraft, and the plasma environment (for E). In addition, spacecraft spin and nutation of sensor booms introduce noise lines at fixed frequencies but with varying amplitudes. Some of the spectral peaks caused by these artifacts are predictable, but they can overlap f _Tn as spacecraft move in L. Also, there are unpredictable peaks that originate from ULF waves that are not toroidal waves. An automated method to determine electron density has been developed by applying a neural-network algorithm to plasma wave spectrograms (Zhelavskaya et al., 2016), and we may design a similar algorithm for f _Tn .

Although significant progress has been made modeling the mass density and field line dependence at geostationary orbit (e.g., Denton et al., 2015; Denton et al.,2016), further work needs to be done to develop an accurate radially dependent model. Neural network analysis might also be helpful here, as it has been for modeling electron density (Chu et al., 2017).

Another need is for event-specific mass density field line dependence. With the possible exception of the event studied by Denton et al. (2009), which had a particularly accurate set of frequencies, significant uncertainty in the observed Alfvén frequencies has precluded an accurate determination of event-specific field line dependence. So, most studies have been statistical (e.g., Denton et al., 2006, Denton et al., 2015; Takahashi and Denton, 2007). Perhaps the mode structure technique described in subsection 3.2 will enable event-specific determination of the field line dependence, but that is yet to be shown.

Finally, we note that magnetoseismology belongs to a family of techniques to probe the magnetospheric plasma structure without using in situ particle measurements. Other techniques include EUV imaging of the plasmasphere (Sandel et al., 2003), energetic neutral atom remote (ENA) sensing of energetic ions (Roelof et al., 1985; Brandt et al., 2005), estimation of n _e using whistler waves (Carpenter and Smith, 1964; Park, 1974), spacecraft potential (Pedersen et al., 1984; Archer et al., 2015), or plasma wave spectra (upper hybrid resonance) (Mosier et al., 1973; Moldwin et al., 2002; Thomas et al., 2021). These indirect techniques are complementary to each other. For example, the global imaging techniques are capable of taking snapshots of plasma structures, which cannot be obtained using the magnetoseismic or n _e techniques unless we have a large number of measurement points. The magnetoseismic and n _e techniques provide the total densities, whereas ENA images provide information on the density of energetic ions (>10 keV). Improvement of magnetoseismic techniques and associated datasets is much desired to enhance the synergy of different density-related techniques.