Nodal Structure of Toroidal Standing Alfvén Waves and Its Implication for Field Line Mass Density Distribution

We have conducted a statistical study of toroidal mode standing Alfvén waves detected by the Van Allen Probes spacecraft in the dayside inner magnetosphere, with an emphasis on the nodal structure of the fundamental through fifth harmonics. We developed a technique to accurately assign harmonic mode numbers to peaks in the power spectra of the electric (Eν) and magnetic (Bϕ) field components of toroidal waves and then determine the spectral intensities of Eν and Bϕ and the coherence and cross‐phase between these field components for each harmonic. The magnetic latitude (MLAT) dependence of these quantities was statistically examined to determine the location of the nodes. In addition to the equatorial nodes located close to the equator (MLAT = 0), we identified several nodes away from the equator within the MLAT range from −20° to +20°. We found that the Eν‐Bϕ cross‐phase is very close to ±90° except near the nodes, indicating that the fixed‐end approximation is appropriate in modeling dayside toroidal waves. Noting that the node latitudes depend on the distribution of the mass density (ρ) along the background magnetic field, we inferred the distribution from the nodes observed at L = 4–6. If we adopt a model field line mass density (ρ) distribution of the form ρ ∝ (1/r)α, where r is geocentric distance to the field line and α is a free parameter, the statistically determined node latitudes indicate that α∼1.5 is appropriate for both the plasmasphere and the plasmatrough.

. The plasma bulk velocity (V) also oscillates because it is related to E and B through the magnetohydrodynamics equation E = −V × B (e.g., Junginger et al., 1984).
Toroidal waves are excited by driver fast mode waves through steady field line resonance (Chen & Hasegawa, 1974;Southwood, 1974) or impulsive triggering (Takahashi & Heilig, 2019). The driver waves may originate from solar wind dynamic pressure pulses (Sarris et al., 2010), foreshock ULF waves (Clausen et al., 2009), or nightside disturbances associated with substorms (Takahashi et al., 1996). Because these source disturbances often have a broad frequency spectrum, they are capable of exciting toroidal waves at multiple harmonics (Cahill et al., 1986;Engebretson et al., 1986;Takahashi & McPherron, 1982;Waters et al., 1994). The presence of harmonics is central to the present study. We refer to the nth harmonic of a toroidal wave as the Tn wave and denote its frequency as f Tn , with n = 1 for the fundamental mode.
The motivation for our mode structure analysis is twofold. First, we are interested in confirming theoretical prediction of the mode structure of toroidal waves in the dayside magnetosphere. Although dayside toroidal waves are usually assumed to be close to the idealized fixed-end modes (perfect reflection at the ionosphere) because of the high ionospheric conductivity (Newton et al., 1978), few studies have examined the structure in a quantitative manner. We compare spacecraft observations with theoretical models.
Second, we are interested in using the mode structure in constraining mass density models. Observed frequencies of toroidal waves have been used to estimate the mass density using a technique known as magnetoseismology (Menk & Waters, 2013). In using this technique, one needs to assume a functional form of the mass density (ρ) variation along the background magnetic field line because the frequencies are related to the integral of the Alfvén velocity over the field line. The field line distribution has important implications for the way the background plasma establishes an equilibrium and how various plasma waves behave (Denton et al., 2006). Previous studies used the frequency of toroidal harmonics to constrain the model (Takahashi & Denton, 2007;Takahashi et al., 2004). Observationally determined nodal structure of the waves is a new approach to the field line mass density distribution. We discuss the results of the mode structure analysis from this point of view.
The remainder of this paper is as follows. Section 2 presents theoretical models of toroidal waves. Section 3 describes experiments. Section 4 presents data analysis in the frequency domain. Section 5 presents statistical results on the frequency and mode structure of toroidal harmonics. Section 6 presents discussion, and section 7 presents the conclusions.

Model Toroidal Waves
To guide the data analysis, we review key properties of theoretical fixed-end toroidal waves. We obtain the frequencies and mode structures of the waves using the toroidal wave equation of Cummings et al. (1969) with the mass density specified by the power law model (Radoski & Carovillano, 1966) where ρ eq is the equatorial mass density, L is the magnetic shell parameter, 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. If α = 0, ρ is constant along the field line. A positive (negative) α means that ρ increases (decreases) from the equator toward the ionosphere. This mass density model has been used widely (Dent et al., 2006;Kabin et al., 2007;Menk et al., 1999;Orr & Matthew, 1971;Price et al., 1999;Waters et al., 1996) not only for its simplicity but also for its relevance to theory (Angerami & Carpenter, 1966). We solve the equation for L = 5, which is representative of observations by Van Allen Probes. Although α is typically assumed to be in the range 0-4 (Orr & Matthew, 1971;Poulter et al., 1984;Takahashi et al., 2004;Vellante & Förster, 2006), we vary α from −6 to +6. We compute the mode frequencies and node latitudes for integer values of α. For non-integer α values, we obtain the frequencies and node latitudes by linear interpolation. One could adopt other field line density models with more free parameters ), but we do not consider such models in this study.
TAKAHASHI AND DENTON 10.1029/2020JA028981 2 of 28 Figure 1 shows the α dependence of f T1 through f T5 . In Figure 1a, the frequencies are normalized to f T1 for α = 0. The normalized frequencies decrease as α increases. This occurs because the integral of ρ over the field line is larger for larger α when ρ eq is held constant, that is, the frequencies are lower when the field lines are "heavier." Although Figure 1a is helpful in understanding how α is physically related to the frequencies, normalization by f T1 is not a good choice in practice when we try to constrain α using the observed frequencies. The reason is that T1 waves are not the easiest to detect, as indicated in Figure 9, which is presented in Section 2. This leads us to use f T3 as the reference for frequency normalization, as was done previously (Takahashi & Denton, 2007). Figure 1b shows the frequencies normalized to f T3 for each value of α. As α increases, f T1 /f T3 and f T2 /f T3 increase while f T4 /f T3 and f T5 /f T3 decrease. The rate of the change of the ratio with α is greatest for T1 waves: a 13.0% change from 0.238 for α = 0 to 0.269 for α = 3. For T2, T4, and T5 waves, the changes are much smaller: 2.4% (from 0.620 to 0.635), −1.1% (from 1.378 to 1.363), and −1.7% (1. 756-1.726). This means that f Tn /f T3 for n ≥ 2 is less sensitive to changes in α than for n = 1.

Mode Frequencies
For a fixed value of α, the variation of the frequency ratios with L is small.
If we take f T1 /f T3 as an example, the ratio decreases from L = 4 to L = 6 by 2.4% for α = 0 and by 2.0% for α = 3. This justifies the comparison made in Section 5.5 between the theoretical frequency ratios for L = 5 and the frequency ratios statistically determined using satellite observations made at L = 4-6. Figure 2 shows the mode structure of T1-T5 waves obtained for α = 1. Each panel shows the normalized amplitudes of the E ν and B ϕ components as a function of magnetic latitude (MLAT in the northern hemisphere regardless of the magnitude of MLAT. As the harmonic number increases, the number of nodes increases, and the MLAT spacing between the nodes decreases. As a consequence, even a spacecraft with a low orbital inclination encounters off-equatorial nodes. Consider Van Allen Probes as an example. The spacecraft had an orbital inclination of 10° and covered the MLAT range from −20° to +20°, as indicated by the black horizontal bar at the bottom of Figure 2. If we take the T5 wave ( Figure 2e) as an example, the spacecraft have the possibility of encountering up to five nodes, excluding those located at MLAT = ±20°. Figure 3 shows the α dependence of the theoretical node latitude (denoted MLAT node ). For both E ν and B ϕ , the non-equatorial nodes move away from the equator as α increases. For example, Figure 3b indicates that the E ν node moves from 10.0° for α = 0 to 17.3° for α = 6. The implication is that we can infer the α value with an accuracy of ∼1, if we can determine MLAT node with an accuracy of ∼1°. Of course, this is contingent upon having a good model for the background magnetic field, which also controls the mode structure. The real magnetic field can be very different from a dipole field, depending on the location and geomagnetic TAKAHASHI AND DENTON 10.1029/2020JA028981 3 of 28 activity. Nevertheless, Figure 3 gives us good motivation to search for the node locations as a means to determine α.

Mode Structures
To evaluate whether the mode structures derived using a dipole field are good enough for comparison with observations, we examine the mode structures for realistic model magnetic fields. Figure 4 shows a dipole and T89c (Tsyganenko, 1989) comparison of the mode structures of T3 waves excited on field lines that pass the dipole equator at geocentric distance 5 R E and at 12 h magnetic local time (MLT). The structures were obtained by numerically solving the wave equation derived by Singer et al. (1981). For simplicity, we considered the case of zero dipole tilt angle. Two T89c fields are considered, one for Kp = 3 (IOPT = 4) and the other for Kp ≥ 6 (IOPT = 7), where IOPT is the input parameter required for the T89c model. We refer to the dipole and the two T89c fields as model 1, model 2, and model 3, respectively. The mass density is given by Equation 1   (a-e) Theoretical mode structure of the E ν (solid line) and B ϕ (dashed line) components of (a) the fundamental (T1) through (e) fifth (T5) harmonics of toroidal waves obtained by solving the toroidal wave equation of Cummings et al. (1969) for L = 5 and α = 1. The black and white dots indicate the nodes of E ν and B ϕ , respectively. The amplitude of each component is normalized to the maximum value. The orange (blue) shading indicates the magnetic latitude (MLAT) domains in which E ν leads (lags) B ϕ in phase by 90°. The nodes excluding those at the ionospheric foot points are numbered starting from the south. The black horizontal bar at the bottom indicates the MLAT range covered by Van Allen Probes. The white vertical dotted line indicates the MLAT location (10.4°) of Van Allen Probe B at 1410 UT on April 22, 2014. the amplitudes normalized as in Figure 2. The mode structures are symmetric about the magnetic equator, so the figure covers only the northern hemisphere (MLAT ≥ 0).
The mode structure differs very little between model 1 and model 2. With model 1, the nodes are found at 10.6° for E ν and at 26.3° for B ϕ . With model 2, the nodes are found at 10.4° and 25.8°. Within the MLAT range covered by Van Allen Probes, only the E ν node will be detected, and the node latitudes differ only by 0.2° between model 1 and model 2. We also examined the mode structures on other field lines (not shown) and found very small changes from those shown. Specifically, on field lines passing the equator at 4 R E , the E ν node is located at 10.5° with model 1 and at 10.3° with model 2. On field lines passing the equator at 6 R E , which is close to the apogee of the spacecraft, the E ν node is located at 10.6° with model 1 and at 10.9° with model 2.
We need to be careful when we consider the mode structures during geomagnetically active periods. The model 3 results shown in Figure 4 indicate a notable departure of the node latitudes from those obtained using model 1 or model 2. With model 3, we find the nodes at lower latitudes, 9.1° (E ν ) and 22.8° (B ϕ ). We also need to pay attention to the local time dependence of the magnetic field, which requires use of realistic models such as T89c. However, we believe that the dipole model is adequate in the present study. The reason is that we discuss the mode structures by statistically processing satellite observations that were made mostly in the noon sector, at L = 4-6, and during quiet to moderately disturbed times (Kp < 3).
We assume that the node locations do not change unless the background magnetic field or plasma mass distribution changes. Compressional Pc5 waves are known to exhibit frequency doubling and waveform deformation near the magnetic equator, which can be explained by a latitudinal oscillation of the equatorial node that is phase-locked to the waves (Takahashi et al., 1987). Frequency doubling and waveform deformation have also been reported for fundamental poloidal waves . We have not seen frequency doubling for toroidal waves and believe that these waves do not induce detectable node oscillations.

Experiments and Data
Data used in this study were acquired by Van Allen Probes A and B. The twin spacecraft made scientific measurements from 2012 to 2019 on elliptical orbits with an inclination of ∼10°, apogee of ∼5.8 R E (maximum dipole L value of 6.6), and orbital period of ∼9 h. The spacecraft were spin stabilized with the spin periods maintained at ∼11 s. The spin axes approximately pointed to the Sun with offset angles lower than 27°. The Van Allen Probe experiments relevant to the present study are the electric field spectrum analyzers (for electron density) and the fluxgate magnetometers included in the Electric and Magnetic Field Instrument Suite and Integrated Science (EMFISIS) experiments (Kletzing et al., 2013) and DC and low-frequency electric field experiments included in the electric field and waves experiments . The E-field data used in the present study were obtained by the spinfit method applied to the two orthogonal components measured in the spacecraft spin plane and rotated into modified geocentric solar ecliptic (MGSE) coordinates . The data have a time resolution of the spacecraft spin period (∼11 s), placing the upper (Nyquist) limit of our wave frequency analysis (f Nyq ) at ∼45 mHz. When the angle λ B between the measured B-field and the spin plane is sufficiently large (>10°), we can use the B⋅E = 0 assumption to get a reliable estimate of the electric field component parallel to the spacecraft spin axis (E xMGSE ). According to Ali et al. (2016), the error in E xMGSE is 0.5 mV/m when λ B = 5.75° and decreases quickly as λ B increases. Our threshold λ B of 10° is lower than the 15° that is normally recommended (e.g., Dai et al., 2013). We use this relaxed threshold because we are mainly interested in finding wave frequencies and need not be very strict about the accuracy of the amplitude of E.
To extract perturbations associated with ULF waves and examine their spectral and polarization properties, we rotate the measured B-field vectors into locally defined magnetic field aligned (MFA) coordinates. The MFA coordinate axes are defined using a reference magnetic field B ref and the spacecraft position vector R relative to the center of the Earth. We obtain B ref by fitting a function to the three components of the observed B-field vector time series given in geocentric solar ecliptic (GSE) coordinates. The function is a polynomial of the form c 0 + c 1 τ + c 2 τ 2 + c 3 τ 3 + c 4 τ 4 , where τ is UT rescaled to the range from −0.5 to 0.5 for the data segment selected for analysis. The segment length is 15 min in the routine moving data window analysis of toroidal waves described below. The coefficients c 0 -c 4 are determined using the least-squares method. In the MFA coordinate system, the parallel or compressional ( By using this method, we can effectively remove slow variations of the magnetic field, either spatial or temporal, without using magnetic field models. We can define B ref by taking running averages of the measured B-field, but we prefer the polynomial fitting technique. In taking running averages one usually selects a data window with a fixed length (T w ). To define perturbations in a time segment [t 0 , t 1 ], the running average must be calculated using B-field data starting at t 0 − T w /2 and ending at t 1 + T w /2. This means that perturbations near t 0 and t 1 are affected by B-field data outside [t 0 , t 1 ]. This can be a problem when we analyze spacecraft data taken near perigee, where the background B-field changes rapidly in a nonlinear manner. By contrast, polynomial fitting does not require data outside [t 0 , t 1 ].
We rotated the spinfit E-field vector samples into an MFA system that uses the spin (∼11 s) averages of the observed B-field to define the field-aligned component. This MFA system was used to produce E-field data files for the entire Van Allen Probes mission period, before the start of this study. As far as wave E-fields are concerned, the difference between this system and the MFA system defined above for wave B-fields is negligible because the E ν or E ϕ values are essentially the same regardless of whether we use the observed B or B ref in applying the B⋅E = 0 assumption in coordinate transformation. To make the spectral properties of the E-field vectors consistent with those of the B-field vectors in the MFA system, we remove the trend, also a polynomial of fourth degree, from E ν and E ϕ .
To compare observations with the theoretical mode structures (see Figure 2), we use dipole coordinates for the spacecraft position. We define L, MLAT, and MLT using the Gauss coefficients relevant to the centered dipole that are specified by Thébault et al. (2015) for the International Geomagnetic Reference Field model.

Frequency-Domain Analysis
To conduct mode structure analysis, we first need to separate harmonics of toroidal waves in the frequency domain. This section describes how we determine f Tn . With the exception of geostationary satellites, spacecraft move in L and MLAT. L shell crossing means continuous changes in f Tn , and changing MLAT means that spacecraft cross E ν and B ϕ nodes (Cahill et al., 1986). In this situation, correctly assigning harmonic mode numbers to multiple spectral peaks is often challenging. In addition, ULF waves TAKAHASHI AND DENTON 10.1029/2020JA028981 6 of 28 other than toroidal waves are present in the magnetosphere. Spectral peaks caused by these waves need to be excluded.
In previous studies using spacecraft data, toroidal wave frequencies were determined by visually inspecting the dynamic spectra of B ϕ (Min et al., 2013;Nosé et al, 2011Nosé et al, , 2015Nosé et al, , 2018Takahashi et al., 2010), E ν , or both (Takahashi et al., 2018). In some studies, the plasma bulk velocity (Takahashi et al., 2014) or ion flux anisotropy (Takahashi et al., 2002) in the azimuthal direction was used as a proxy to E ν . Our approach is similar. However, to improve the accuracy of mode identification, we also compute the cross-phase between E ν and B ϕ and compare it with a theoretical model.

Dynamic Spectra
These are derived from the elements of the 2 × 2 spectral matrix constructed from the Fourier transform of input E ν and B ϕ time series as described by Bendat and Piersol (1971). We compute these spectral parameters in a 15 min data window, which is shifted in 5 min steps. The window size gives a frequency resolution of 1.1 mHz for the raw Fourier components. A three-point (3.3 mHz) boxcar smoothing in frequency is applied to the elements of the spectral matrix before deriving the spectral parameters displayed. Figures 5f and 5g show the λ B angle defined in section 3 and the electron density n e derived using plasma wave spectra . The condition λ B > 10° is satisfied except after 1700 UT, which means that the E ν data are mostly reliable. The rapid changes of n e occurring at ∼1000 UT (L ∼ 3.8) and ∼1530 UT (L ∼ 4.6) indicate plasmapause crossings. The n e variation is very smooth both inside and outside the plasmapause.
Multiharmonic toroidal waves are evident in all spectrograms. The waves produce multiple spectral lines, which change smoothly in frequency but vary significantly in intensity. The most persistent spectral line starts at ∼1040 UT (L ∼ 4.6) at 30 mHz, reaches a minimum frequency of 18 mHz at ∼1300 UT (L ∼ 5.9), and lasts until ∼1520 UT (L ∼ 4.9). This spectral line is attributed to T3 waves, as described below. T3 waves reappear after ∼1600 UT and indicate a rapid frequency change as the spacecraft moves inward within the plasmasphere (L < 4.1). T1 and T2 waves follow the same trend. The toroidal waves are similar to those observed by other elliptically orbiting spacecraft (Cahill et al., 1986;Clausen et al., 2009;Sarris et al., 2009;Takahashi et al., 2004). Features to be ignored in this example are the elevated    Takahashi et al. (2015), the most persistent noise appears between 20 and 25 mHz, and the intensity of the noise varies from orbit to orbit and also during an orbit. This noise limits our ability to determine toroidal wave frequencies using E ν . Fortunately, this is not a very serious problem in our statistical analyses of toroidal waves because the noise occupies a small fraction of the frequency range surveyed. Also, there are many orbits on which the noise is very weak or nonexistent, and the B ϕ component is free from this noise line. We

Sample Data Segment
We explain how we determine the harmonic modes using a sample 20 min data segment. The segment is marked by the black bar at the bottom of Figure 5e. The spacecraft was located at L∼5.7 and MLAT∼10°. Figures 6a and 6b show that E ν and B ϕ oscillate with different waveforms. The E ν oscillation consists of a TAKAHASHI AND DENTON 10.1029/2020JA028981 7 of 28 strong (2 mV/m peak-to-peak) long period (∼3 min) component and weaker shorter-period components. The B ϕ oscillation is dominated by a 1 min component with peak-to-peak amplitudes of ∼1 nT, but the oscillation is not purely monochromatic.  Figure 2c shows that an E ν node is located at MLAT = 10.6°, very close to the spacecraft, making a prediction of    E B for the T3 wave difficult. However, the observed    E B of ∼ 90° ( Figure 6g) implies that the spacecraft was located slightly south of the E ν node. The proximity of the node to the spacecraft can be inferred from the missing T3 peak in   E E S (Figure 6c).

Interactive Identification of Mode Frequencies
We have developed a routine procedure to identify toroidal wave frequencies and harmonic mode numbers using Van Allen Probes data. Automated frequency detection procedures have been developed for ground magnetometer data (Berube et al., 2003;Chi et al., 2013;Sandhu et al., 2018;Wharton et al., 2018), incorporating the cross-spectral analysis technique (Waters et al., 1991). This approach cannot be used for spacecraft data because no spacecraft pair maintains the short radial distance that is required for the technique to work properly. Fortunately, we can determine toroidal wave frequencies using data from individual spacecraft. As Figure 5 shows, dynamic spectra generated from spacecraft data exhibit clear spectral peaks arising from toroidal waves. A major difficulty is assigning a harmonic mode number to each peak. To reduce this difficulty, we pay attention to the nodal structures of toroidal waves.
Figure 7 illustrates our approach using Van Allen Probe B data for orbit 1596. Figure 7a shows that the spacecraft was mostly in the prenoon sector, where toroidal waves in the Pc3-4 band (7-100 mHz) are routinely detected (Anderson et al., 1990). Our first step is to find spectral peaks that are potentially associated with toroidal waves. We do this by simply searching for peaks in S ij computed using the same 15 min moving data window as that used in generating Figure 5 . If a spectral peak is found at frequency f p , we examine whether we can define the full-width-at-half-maximum (FWHM, delimited by f 0 and f 1 ) around the peak, within the entire band covered by the spinfit data, 0-f Nyq . If we can, and there are no other spectral peaks between f 0 and f 1 , we compute the weighted frequency f w , given by This frequency is a good representation of the wave frequency when a wave has a bandwidth spanning several Fourier components (Takahashi & McPherron, 1982). We then integrate   from f 0 to f 1 and take the square root of the integral to obtain the root-mean-square amplitudes δE ν , δE ϕ , δB ν , δB ϕ , and δB μ . We similarly define the band integral of the E ν -B ϕ coherence and cross-phase, denoted     

E B
and      E B , respectively. If there are multiple peaks between f 0 and f 1 , we discard the spectral peak at f p to avoid confusion in the subsequent analyses. The large colored dots in 7b are the f w data points that resulted from the peak search in    Figure 5. The small black dots in Figure 7b indicate the f w samples that survived the downselection procedure.
In the third step, we generate a plot of model wave frequencies with information on the E ν -B ϕ cross-phase included to guide manual determination of the harmonic mode numbers. In the example shown in Figure 7c, the frequencies and cross-phase are obtained by solving the dipole field toroidal wave equation (Cummings et al., 1969). Here we adopt α = 1 for the mass density model and an average ion mass (M, given by ρ/n e ) of 4 amu. Once this setup is done, the n e data shown in Figure 7d give the Alfvén velocity along the field line that is necessary to solve the wave equation. To reduce computation, the equation is solved at discrete L shells (2, 3, …, 7) for a fixed equatorial mass density (1 amu cm −3 ), and the wave frequencies and the node latitudes are saved in a lookup table. The theoretical frequencies at the spacecraft at any UT epoch are obtained by L interpolation of the frequencies saved in the lookup table and then adjusting the frequencies by the mass density at the spacecraft given by Mn e . The node latitudes are also interpolated. From the MLAT of the spacecraft (Figure 7e) and the theoretical node latitudes, we can predict the sign of    E B at the spacecraft, which is color coded in Figure 7b. We caution not to take the model frequencies too literally, because M varies with solar activity, geomagnetic activity, and L (Denton et al., 2011;Nosé et al., 2015;Takahashi et al., 2004). Also, the real magnetic field will be very different from the dipole model when the geomagnetic activity level is high, resulting in large differences between the observed and model frequencies.
In the fourth step, we visually compare Figures 7b and 7c and use an interactive computer tool to select f w samples that belong to a target harmonic mode. The tool allows one to use a computer mouse to draw a rectangular area called "marquee" in an orbital f w plot displayed on a computer screen in the format of Figure 7b and to save all data points in the area into a temporary data array by choosing the "save" option from a pulldown menu. The saved data points appear on the screen with a different symbol. If one selects wrong data points by mistake, they can be removed from the temporary array by invoking the "remove" option. We typically repeat the marquee procedure several times to cover all possible f w samples belonging to a target harmonic mode. When the selected data points appear satisfactory for the entire orbit, the data in the temporary array are written out into an output file with the file name specifying the spacecraft, orbit number, and the harmonic mode.
In the present example, the target is the third harmonic (T3), which is highlighted by the thick line in Figure 7c. The model predicts a switch of the sign of    E B at 1418 UT from positive to negative, in association with the satellite crossing of an E ν node MLAT = 10.6°. In Figure 7b, we find a sequence of f w samples that fit the model, in terms of both the frequency trend and the sign of the cross-phase. Because the real field line mass density distribution is unknown, the model prediction of the node latitudes based on a single α value should have some error relative to the observation. However, we expect the error to be of the order of a few degrees at most and find that the predicted cross-phase plotted in the format of Figure 7c is quite helpful to determine harmonic modes. The procedure is repeated for the T1 and T2 modes. The downselected frequencies are denoted f wT1 , f wT2 , and f wT3 .
Our frequency selection procedure involves human judgments, and the results certainly depend on the individual (operator) who runs the procedure. To make proper judgments, the operator needs to be familiar with the properties (e.g., instrument noise) of the spacecraft experiments and magnetospheric ULF waves, and the outcome depends on the experience of the operator. However, we believe that our results on the mode frequencies and node latitudes depend little on the operator because the judgment errors are likely random and averaged out in statistics. In case studies, we can afford paying attention to far more details of the wave properties to determine toroidal wave frequencies with high confidence (e.g., Takahashi et al., 2015). A possible approach in the future to remove the operator dependency is to generate high-confidence frequency data files for a subset of Van Allen Probes orbits and use the files to develop a neural network technique to automatically determine the frequencies. This approach has been taken in the Neural-network-based Upper hybrid Resonance Determination technique (Zhelavskaya et al., 2016) developed for automated determination of electron density using plasma wave spectra.

Statistics
We determined f wT1 , f wT2 , and f wT3 , using Van Allen Probe A and B data for January-June 2014. We use the resulting data set to statistically evaluate the frequencies and mode structures of toroidal harmonics.

Overview
Figures 8a-8c show the locations of wave detection. Each data point is the spacecraft position at the center of the 15 min data window mapped to the equatorial plane of solar magnetic (SM) coordinates along the dipole field line. The samples for each harmonic mode are distributed in the noon sector and in the L range from ∼3 to 6.6 (geostationary orbit). Figures 8d-8f show time series plots of all f wT1 , f wT2 , and f wT3 samples. The vertical spread of the data points results mostly from satellite motion in L. On each orbit, the frequencies change with L as shown in Figure 7.
The f wT1 samples are found mostly between 4 and 20 mHz (Figure 8d). The lower cutoff occurs because we use a 15-min data window for spectral analysis. The frequency resolution of the Fourier transform with this data window is 1.1 mHz. To define a spectral peak after taking three-point smoothing, we need to have a few Fourier components on either side of the peak, which is the reason we do not have f w values lower than ∼4 mHz. The absence of f wT1 samples above ∼30 mHz results from our inability to detect T1 waves in the 15 min data window when a spacecraft is located at L < 3 where the frequency changes rapidly as the spacecraft moves in L.
The f wT2 samples are more evenly spread out in frequency (Figure 8e). There is no obvious low frequency cutoff originating from the finite data window length. The variation of the lower cutoff of f wT2 , seen most clearly in April and May, is caused by changes in the space environment. As the geomagnetic activity increases (see Dst in Figure 8g), the plasmapause shrinks, and the spacecraft are in the plasmatrough when located outside of L∼4. At a given L, the frequencies are higher in the plasmatrough than in the plasmasphere because the mass density is lower in the plasmatrough.
The f wT3 samples are distributed similarly to f wT2 (Figure 8f). A notable difference is that f wT3 reaches the upper cutoff (f Nyq , 45 mHz) of the 11 s data. An example of the cutoff is found in Figure 5b at 1610 UT, when the spacecraft was in the plasmasphere and moving inward.

Plasma Domains
We are interested in separating the plasmasphere and plasmatrough because ions in these domains come from different source regions and go through different transport processes. Toroidal waves could provide valuable information on the difference between the domains, for example, regarding the ion composition and field line mass density distribution. We could use n e data to identify the domain. However, the two plasma domains can be distinguished using a statistical approach without relying on the n e data. Figure 9 shows the number (given in color code) of the f wT1 (left column), f wT2 (middle column), and f wT3 (right column) samples in L versus frequency bins, sorted by the source spectrum (rows). The f wT1 distributions form a group that resides in the L range 2.0-4.5 with the frequency falling from ∼20 to ∼6 mHz and another group that resides in the L range 4.0-6.5 with the frequency falling from ∼15 to ∼6 mHz. The f wT2 and f wT3 distributions show similar patterns. The distinction between the two groups is most clearly seen in Figure 9f  The interpretation of the two groups is easy. The low-L group comes from measurements made in the plasmasphere, and the high-L group comes from measurements made in the plasmatrough. The two groups overlap in L because the plasmapause distance depends on the solar wind and geomagnetic conditions. The large vertical (frequency) spread of the samples in the second group can be explained by larger variability of the mass density outside the plasmasphere than inside.
To distinguish the two f wT3 groups, we introduce a model boundary frequency f bT3 : which is shown by the dotted line in Figure 9f. We assumed a power-law L dependence following previous studies of toroidal waves detected by spacecraft (e.g., Takahashi et al., 2014). This particular equation was obtained from a least-squares fit of a power-law function to several data points taken along the trough in the 2-D distribution shown in Figure 9f. In the following analyses, measurements are judged to be made in the plasmasphere (plasmatrough) if f wT3 < f bT3 (f wT3 ≥ f bT3 ).
We find that the largest number of frequency samples is obtained using   B B S for T3 waves (N = 15,565).
This was also the case in studies that used magnetic field data from geostationary satellites ( Takahashi &   TAKAHASHI (Fairfield, 1971) and a circle (L = 6.6) representing geostationary orbits are included.
(d-f) All f wT1 , f wT2 , and f wT3 samples. (g) Dst index. Denton, 2007;Takahashi & McPherron, 1982). This leads us to use f wT3 as the key parameter in our statistical data analysis.

Frequency Ratios and Field Line Mass Density Distribution
Frequency ratios among toroidal harmonics depend on the field line mass density distribution, as illustrated in Figure 1. This dependence has been used to infer the distribution from observed toroidal wave frequencies (Takahashi & Denton, 2007;Takahashi et al., 2004). We took the same approach to the Van Allen Probes data and show the results in Figure 10. frequencies are shown for the plasmasphere and plasmatrough as time series (Figures 10a and 10c) and as histograms (Figures 10b and 10d).
We defined the f wTn /f wT3 domains (bands) occupied by T1-T5 waves by paying attention to notable features of the total f wTn /f wT3 histogram for the plasmatrough (Figure 10d). The bands are 0.10-0.36, 0.36-0.82, 0.82-1.19, 1.19-1.55, and 1.55-1.92. The T1 band was defined based on the lower cutoff (0.10) and the first minimum (0.36) of the histogram. The limits for the remaining bands are defined to be the midpoints between the neighboring peaks of the histogram. We apply the same band definition to the plasmasphere because the total f wTn /f wT3 histogram for the plasmasphere are peaked very near the center of the bands. In each f wTn / f wT3 band, we calculated the mean (denoted f wTn /f wT3 ) and standard deviation of the mean (denoted σ mean ). The results are summarized in Table 1 along with the number of samples (N). The mean values are also shown to the right of Figures 10b and 10d. We do not define f wT4 /f wT3  in the plasmasphere because the total histogram does not exhibit a clear peak.
The frequency ratios are similar to those obtained in a previous study     Table 1. The filled red circles mark the intersections of the red lines with the black curves, which give α estimates labeled "α freq_obs ." As an indicator of the error, we evaluated the α values corresponding to f wTn /f wT3  ± σ mean and show them using short red horizontal bars. The errors are small compared with the difference in α freq_obs between different harmonics. In both plasma domains, the estimated α freq_obs varies among the harmonic modes, which means that the mass density model given by Equation 1 should be taken only as first approximation.
Given the harmonic dependence of α freq_obs , we define an optimal α value by minimizing is the theoretical ratio given as a function of α. With this approach, which follows Takahashi and McPherron (1982), we obtain α = 1.8 for the plasmasphere (Figure 11a) and a nearly identical result, α = 1.9, for the plasmatrough (Figure 11b).
Regarding the differences of α freq_obs among the harmonics, we examine the accuracy of f wT1 /f wT3 . In Figures 9a, 9d, and 9g, we find a 4-mHz lower cutoff to the lower-frequency (plasmasphere) group of f wT1 data points for L > 4. This is an artifact of the finite length of the data window for Fourier transform, as described in Section 5.1. The cutoff results in a sharp lower cutoff of the distribution of the normalized frequency for the plasmasphere at 0.16 (Figure 10b). A cutoff occurs in the plasmatrough distribution as well but at a lower value of 0.10 ( Figure 10d). This difference can be explained by the fact that toroidal wave frequencies are higher in the plasmatrough; most of them are higher than 4 mHz. As a consequence of this skewing, f wT1 / f wT3  is higher for the plasmasphere. To avoid this skewing, we applied Equation 4 only to f wT2 /f wT3  and f wT5 /f wT3  and obtained α = 2.0 (green horizontal dashed line in Figure 11a). We also replaced the f wT1 /f wT3  value with that for the plasmatrough and obtained α = 1.2 (blue horizontal dashed line in Figure 11a). This exercise indicates that the skewing of the f wT1 /f wT3 distribution does not produce a large error in estimating plasmaspheric α. Figure 12 illustrates how we derive the spectral parameters and use them to detect the nodes of toroidal waves. We focus on T3 waves using data from Van Allen Probe B acquired on orbit 1596 (Figures 5 and 7) as an example. Our desire is to continuously track the behavior of the spectral parameters at the T3 frequency. If the waves were continuously detected, this would be easy. The reality is that toroidal waves often become undetectable either from node crossings or from disappearance of driver fast mode waves. In Figure 12a, f wT3 samples are found from 1055 to 1520 UT, with a gap at ∼1315 UT. By comparison, f wT1 and f wT2 samples are sparse but are present outside the time span of the f wT3 samples.

Example of a Node Crossing
When f wT3 is not directly determined from power spectra, we can still define a T3 frequency using the T1 or T2 frequencies, which are included in the frequency survey illustrated in Figure 8. This is done first by estimating f wT3 from f wT1 or f wT2 using the frequency ratios found in Figure 10d for the plasmatrough: We use the ratios for the plasmatrough because we have a larger number of f wTn samples from that region than from the plasmasphere. We have confirmed that adopting the plasmaspheric ratios does not change the outcome of the mode structure analysis described in Section 5.5. The plus symbols in Figure 12b indicate both the f wT3 values directly derived and those estimated using f wT1 or f wT2 . Given the possibility of multiple f wT3 estimates at a time step, we obtain a unique value simply by calculating the mean of the f wT3 values defined above. The mean is denoted f T3_mean and is plotted by black filled circles in Figure 12b. We use f T3_mean as the key to define the mean frequencies of other harmonics. Adopting the plasmatrough f wTn /f wT3  values listed in Table 1, we get the following mean frequencies: f T1_mean = 0.24f T3_mean , f T2_mean = 0.62f T3_mean , f T4_mean = 1.36f T3_mean , and f T5_mean = 1.73f T3_mean . These frequencies are used in Section 5.5 to statistically evaluate the mode structure of toroidal waves.  Figure 10. The red filled circles labeled "α freq_obs " indicate the α values corresponding to the observationally determined frequency ratios. The red horizontal bars indicate the α values that take into account the standard error of the mean listed in Table 1 This example indicates that node crossings can be used to gain information on field line mass distribution. The node latitude may be a superior variable to look at than the frequency ratios (Takahashi & McPherron, 1982) for the distribution because wave frequencies cannot be very accurately determined at each time step. In the Takahashi and McPherron (1982) study, the length of the data window was 1 h, 4 times that of the present study, but the estimated α still varied between 0 and 4, very likely because of the uncertainly in the wave frequencies. On the other hand, detection of a node is dictated by the spacecraft orbit. A clean node crossing like the one illustrated in Figure 12 occurs only a few times per orbit at most. This make a statistical approach essential to understand the mode structure of toroidal waves.

Mode Structure Statistical Analysis
We have applied the spectral analysis method illustrated in Figure 12 to T1-T5 waves at Van Allen Probes A and B during January-June 2014. This resulted in orbital data files containing   that are evaluated at the five frequencies f T1_mean , …, f T5_mean defined in Section 5.4. We demonstrate in Figure 13 that we can routinely determine the most important parameter for the mode structure analysis,    E B , using the T3 waves observed in a 4 day period in January 2014. Figure 13a shows all f T3_mean samples from both spacecraft. Figure 13b shows    E B for cases with high coherence (     0.5

E B
). The data points in this figure are clustered around ±90°, with the sign being consistent with the MLAT of the spacecraft shown in Figure 13c. That is,    E B is ∼−90° (∼90°) when the spacecraft is in the MLAT domain shaded blue (orange).
We used the full 6 months data to generate plots of the MLAT dependence of the spectral parameters, which are shown in Figure 14 for the plasmasphere and in Figure 15 for the plasmatrough. Note that we limited the L range to 4-6 for comparison with the theoretical models for L = 5 (Figure 3).  (Figure 2a). Accordingly, we judge that a B ϕ node is on average located within the MLAT bin centered on 0 (MLAT node_obs = 0) and indicate the judgment by a blue vertical line.
A cautionary statement is warranted here regarding the behavior of   E E S . Going back to the theoretical T1 mode structure shown in Figure 2a, we find that the E ν amplitude is actually peaked at MLAT = ±42° and has a minimum value (83% of the peak value) at the equator. This is inconsistent with the equatorial peak of   E E S appearing in Figure 14a. This discrepancy could be explained by broadband equatorial E ν oscillations that are not associated with toroidal waves. We speculate that fast mode waves generated by irregular solar wind dynamic pressure variations are the source of such E ν oscillations. Mixture of toroidal and non-toroidal waves is inevitable in our analysis because we evaluate   In the results for T2 waves (Figures 14f-14j), we also find a generally good agreement with theory (Figure 2b). There is an equatorial minimum of     / E E B B S S (Figure 14h), which is colocated with the minimum of    E B (Figure 14i) and a positive-to-negative change of    E B (Figure 14j). Accordingly, we judge that an E ν node is located at the equator, that is, MLAT node_obs = 0. However, we note that   E E S does not exhibit an equatorial minimum expected from the E ν node predicted by theory (Figure 2b), possibly because of the non-toroidal waves we speculated in discussing the T1 wave results. The non-toroidal waves are also likely the cause of the gradual change of    E B at the equator.
The equatorial nodes found for T1 and T2 waves are valuable confirmation of the basic concept of standing Alfvén waves and their symmetry about the magnetic equator (Sugiura & Wilson, 1964). Equatorial nodes are expected regardless of the harmonic modes provided the field line mass distribution and the background magnetic field are symmetric about the equator. Therefore, the equatorial nodes do not provide information on how the mass density changes with distance from the equator. The mode structure analysis yields more interesting results as the mode number increases because we begin to detect nodes away from the magnetic equator that can be used to infer the mass density distribution.
For T3 waves (Figures 14k-14o  The results for T4 waves (Figures 14p-14t) are puzzling. For example,    E B changes from ∼−90° to ∼90° across MLAT = 0, which is opposite to what we expect from the theoretical model (Figure 2d). In addition, we find no evidence for the B ϕ nodes that theory predicts at MLAT ∼±8°. We can explain these findings by low intensity of T4 waves. Figure 10b indicates that the spectral peak detection rate is low and that a clear histogram peak is missing in the T4 band. Weak T4 waves can be masked by non-toroidal waves. Also, it is TAKAHASHI AND DENTON 10.1029/2020JA028981 20 of 28 possible that signals from T3 and T5 waves leak into the T4 band masking T4 waves. This latter explanation is plausible because the negative-to-positive    E B change across MLAT = 0 seen in Figure 14t is qualitatively the same as those found for T3 ( Figure 14o) and T5 (Figure 14y) waves.
The results for T5 waves (Figures 14u-14y In the plasmatrough (Figure 15), T1, T2, T3, and T5 waves show node latitudes that are generally similar to those found in the plasmasphere. The    E B plots for these waves indicate an equatorial node (Figures 15e, 15j, 15o, and 15y). In addition, we find E ν nodes of T3 waves, MLAT node_obs = −11° and MLAT-node_obs = 12°, the latter of which matches the example shown in Figure 12, and B ϕ nodes of T5 waves, MLAT node_obs = ±13°. A difference from the plasmasphere results is that B ϕ nodes of T4 waves are visible at MLAT node_obs = ±8° (Figures 15r-15t), which is not surprising because Figure 10d shows clear presence of T4 waves in the plasmatrough. Another difference is that no E ν nodes are selected for T5 waves. Figures 14 and 15 both indicate small N-S asymmetries in the mode structure. In the plasmasphere, the detected E ν nodes of T3 and T5 waves are located at latitudes that are 1° or 2° higher in the northern hemisphere than in the southern hemisphere. The same is true in the plasmatrough for T3 waves. This asymmetry could be explained by assuming that the effective magnetic equator is located slightly (∼0.5°) north of the dipole equator (MLAT = 0). Note that the equatorial nodes identified in Figures 14 and 15 come with an uncertainty of 0.5° because the MLAT bins for the statistics have a width of 1°. Alternatively, the asymmetry could be attributed to an asymmetry in the background magnetic field, mass density distribution, or the ionospheric boundary condition. It is also possible that the asymmetry is related to the MLAT bias in sampling the northern and southern hemispheres, combined with the MLT dependence of the background magnetic field or field line mass distribution. The dipole tilt angle might also induce a small north-south (N-S) asymmetry to the background magnetic field. In a study of compressional Pc5 waves, mostly detected at L∼8 in the dawn and dusk sectors, Takahashi et al. (1990) suggested that the equatorial node of the waves is shifted from the dipole equator by as much as 9° when the tilt angle is large. A similar tilt angle effect might cause a node shift of dayside ULF waves. We hope to conduct a detailed analysis of the asymmetry in the future when we have processed Van Allen Probes data for more orbits. tionship for T3, T4, and T5 waves repeated from Figure 3 in a modified format. The filled circles indicate MLAT node_obs and the corresponding α, which is denoted α node_obs . The crosses indicate MLAT node_obs ± 1°, in consideration of the MLAT bin width used in the statistics, and the corresponding α, which serve as a measure of the uncertainties in α node_obs . In this case, the uncertainty of α node_obs is ∼±1.5.

Field Line Mass Density Distribution Inferred From Node Latitudes
Ignoring the uncertainties, we find that α node_obs varies from one node to another and also between the plasmasphere and the plasmatrough. This situation is similar to α estimation from the frequency ratios ( Figure 11) and implies that Equation 1 does not completely describe the variation of density along the field line. Nevertheless, we can define an optimal α value by minimizing This equation gives α = 1.4 when applied to data for the plasmasphere (Figures 16a-16c, thick black horizontal dashed lines) and α = 1.7 for the plasmatrough (Figures 16g-16i). Noting the N-S difference of the node latitudes, we shifted the magnetic equator northward by 0.5° to generate additional figures for the plasmasphere (Figures 16d-16f) and the plasmatrough (Figures 16j-16l). This adjustment does not change the α values from Equation 5 for either plasma domain.
The α values obtained from the node latitudes are close to those obtained from the frequency ratios (Figure 11). The various α values from the two methods are in the narrow range 1.2-2.0 and do not show any significant differences between the plasmasphere and plasmatrough. Based on these results, it appears that α∼1.5 describes the field line mass distribution for both plasma domains.
While the good agreement between the two methods implies that either method can be used for mass density modeling, we note that there are advantages in using the node latitude method. First, a clear node crossing such as that shown in Figure 12 can be used to model the density more accurately than is possible using the frequency ratios. Second, the node latitudes provide information on the possible N-S asymmetry of the mass density distribution. Third, the number of parameters for constraining mass density models can be larger in the node latitude method. Using Van Allen Probes data, we are able to identify six off-equatorial nodes, equal to the number of harmonic frequencies identified using Geostationary Operational Environmental Satellite (GOES) magnetometer data (Takahashi & Denton, 2007). If we use data from satellites reaching higher MLAT than Van Allen Probes, we will be able to determine more node latitudes and use that information to construct mass density models that are more structured than Equation 1.

Significance of Mode Structure Analysis
Numerous previous studies have shown the validity of the concept of standing Alfvén waves in terms of the N-S conjugacy on the ground (Sugiura & Wilson, 1964), a±90°phase delay between E and B fields measured in space (Cahill et al., 1986), and the diminishing amplitude of a field component near the magnetic equator (Singer & Kivelson, 1979).
Our results are new in that they include more quantitative details on the mode structure for multiple harmonics. For example, the equatorial nodes are confirmed to be very close (within 1°) to the dipole equator. Nodes away from the equator are detected, and their locations indicate high degree of symmetry (also within 1°) about the magnetic equator. Another important finding is that the cross-phase between E ν and B ϕ is close to ±90° except near the nodes. This means that the time-averaged Poynting flux along the background magnetic field is small and confirms that the fixed-end approximation (Newton et al., 1978) is valid for describing the mode structures observed by spacecraft. Therefore, our statistical analysis confirms that the existing theoretical models of toroidal waves in a dipole magnetosphere (Cummings et al., 1969;Radoski & Carovillano, 1966) fairly accurately describe the mode structure of toroidal waves excited at L ∼ 5 under average geomagnetic conditions.

Uncertainty of α Estimated Using Node Latitudes
The information we gain on the mode structure allows us to take a new approach to infer the field line mass density distribution. Nodes located off the equator but at MLATs accessible by spacecraft are of particular importance. In this regard, T1 waves are irrelevant because the waves have only one node in space (at the magnetic equator) regardless of α. The same is true for T2 waves because Van Allen Probes did not TAKAHASHI AND DENTON 10.1029/2020JA028981 24 of 28 go beyond MLAT = ±20° and thus were not able to fully straddle the off-equatorial B ϕ nodes expected at |MLAT| > 19° (for α ≥ 0) according to Figure 3a. The third or higher harmonics are valuable because they have off-equatorial nodes within the MLAT range of Van Allen Probes.
Not all off-equatorial nodes carry equal significance for estimating α. This has an impact on our ability to estimate α from observed node latitudes. To understand this, let us take a look at the MLAT node -α relationship for T5 waves shown in Figure 16c. In the case of the lower-latitude E ν node located in the northern hemisphere (node number 6 in Figure 2e), MLAT node changes from 5.6° for α = 0 to 10.3° for α = 6, an increase of 4.7°. In the case of the B ϕ node located in the northern hemisphere (node number 7 in Figure 2e), MLAT node changes from 11.8° for α = 0 to 20.9° for α = 6, an increase of 9.1°. Accordingly, for a given error in MLAT node_obs , the corresponding α node_obs has a smaller error with nodes located at the higher latitudes. Therefore, it will be advantageous to use spacecraft that reach high MLATs to expand on the α estimation technique introduced in this study. For example, Figure 2 indicates that with a hypothetical satellite reaching MLAT = ±35° (orbital inclination of 25°), we would be able to detect nodes 1 and 3 of T2 waves, nodes 1 and 5 of T3 waves, nodes 1 and 7 of T4 waves, and nodes 1, 2, 8, and 9 of T5 waves, in addition to those detected by Van Allen Probes.
On the other hand, it is also true that the wave amplitude becomes smaller as the harmonic mode number increases. This makes it difficult to accurately determine the node latitudes for higher-order harmonics. Therefore, there is a trade-off between using higher-order harmonics for α sensitivity and using lower-harmonics for better statistics of node latitudes. We have not discussed rigorous statistical approaches to address this issue. Nonetheless, the simple approaches described above demonstrate that we can obtain valuable information on the field line mass density distribution using the node latitudes detected by spacecraft.

Comparison with Previous α Estimation
Because our use of node latitudes for mass density modeling is new, we compare our results with previous results on α obtained using frequency analysis. Denton et al. (2006) statistically analyzed toroidal wave frequencies observed by CRRES and found that α∼2 is appropriate for L = 4-5 and α∼1 is appropriate for L = 5-6. This is in agreement with our results within a margin of ±0.5.
In the L > 6 region, the GOES satellites (L ∼ 7, usually the plasmatrough) have been a major data source for magnetoseismic studies Takahashi & Denton, 2007;Takahashi et al., 2010). Because the satellites provide a large volume of frequency data extending even to the nightside, it is possible to examine the MLT dependence of α. Denton et al. (2015) found that α changes from ∼0 at noon to ∼2 at dawn. Such MLT dependence in the L > 6 can be addressed by expanding the Van Allen Probes f Tn survey to locations outside the noon local time sector.
The power-law function (Equation 1) is also used to model field line distributions of electron density n e . A series of papers used n e data from the Polar spacecraft to estimate α for electrons Denton, Goldstein, Menietti, Young, 2002;Goldstein et al., 2001). Because the spacecraft was on a polar orbit, n e measurements were made over a wide range of MLAT, and the field line distribution was determined by binning data in L and MLAT. These studies reported α in the range 0-1 in the plasmasphere and in the range 1.6-2.1 in the plasmatrough. These results are remarkably close to our results for ρ derived using the node latitudes. Because different ion species (e.g., H + and O + ) could have different density variation along field lines, α does not need to be the same for ρ and n e . Thus, the Polar results might indicate that H + and O + do have similar field line dependence for L < 6.

Conclusions
Van Allen Probes routinely detected several harmonics of toroidal waves within the frequency range covered by spinfit data, which are well suited for studying magnetohydrodynamic waves. We developed a technique to identify the frequency of the harmonics and to derive spectral parameters for each harmonic that are related to the latitudinal mode structure of the harmonics. We conclude the following: 1. Wave frequencies can be classified into plasmasphere and plasmatrough groups 2. Spectral parameters for the harmonics exhibit clear MLAT dependence, from which we can determine the location of the nodes (node latitudes) of the harmonics 3. If we assume a r −α mass density variation along field lines, the node latitudes indicate that α∼1.5 is appropriate for both the plasmasphere and the plasmatrough in the L = 4-6 region where the majority of toroidal waves were detected

Data Availability Statement
Data used in this study are available from the following sources: NASA Goddard Space Flight Center Space Physics Data Facility Coordinated Data Analysis Web (https://cdaweb.gsfc.nasa.gov) for Van Allen Probes; World Data Center for Geomagnetism, Kyoto (http://wdc.kugi.kyoto-u.ac.jp) for the Dst index.