Electromagnetic waves emitted by ground-based very low frequency (VLF) transmitters can penetrate through the ionosphere into the plasmasphere, which becomes one important type of whistler mode waves in the inner magnetosphere. In the ionosphere, VLF transmitter emissions can cause ionospheric perturbations, particle precipitation and heating (Inan et al., 2007; Parrot et al., 2007; Parrot et al., 2009; Gamble et al., 2008; Sauvaud et al., 2008; Graf et al., 2009; Marshall et al., 2010; Bell et al., 2011; Li et al., 2012; Němec et al., 2020). In the inner magnetosphere, VLF transmitter waves can resonate with energetic electrons and contribute to the loss of energetic electrons in radiation belts (Gamble et al., 2008; Sauvaud et al., 2008; Graf et al., 2009; Rodger et al., 2010; Li et al., 2012; Ma et al., 2017; Ross et al., 2019; Claudepierre et al., 2020; Cunningham et al., 2020; Hua et al., 2020; Hua et al., 2021; Hua et al., 2022; Liu et al., 2022). The propagation of VLF transmitter signals in the magnetosphere can be divided into two main categories: ducted mode (Helliwell, 1966) and non-ducted mode (Cerisier, 1973). The ducted mode requires the presence of plasma density irregularities, which can confine the wave normals to nearly align with the background magnetic field line. The non-ducted mode occurs in smoothly varying density distribution, and the wave normal angles can become very oblique.

Different propagation modes of VLF transmitter signals result in different distributions of wave power and wave normal in the inner magnetosphere, thereby affecting the loss process of energetic electrons. Starks et al. (2008) built a VLF transmitter model combining a simulation of the fields in the Earth-ionosphere waveguide, ionospheric absorption estimates, and geomagnetic field and plasma density models with fully three-dimensional ray tracing. Comparing the modeled wave power with the observations from 5 satellites, the model systematically overestimates the median field strength in the plasmasphere by approximately 20 dB at night and at least 10 dB during the day. Clilverd et al. (2008) compared the distribution of transmitter signals from DEMETER and CRRES in ionospheric footprints and indicated that the non-ducted mode dominates in the low L region ( < 1.5), and transmitter propagation becomes highly ducted at higher L-shells. Using years of observations from Van Allen Probes and DEMETER satellites and simulation from a ray-tracing model, Zhang et al. (2018) revealed the propagation route for some important VLF transmitters and provided direct evidence of the dominant non-ducted propagation. The burst mode waveform measurement from Van Allen Probes can provide direct observation of the wave normal of the Alpha transmitter signals, and the corresponding statistical study has shown that the non-ducted mode dominates over ducted propagation in both the occurrence and wave intensity (Gu et al., 2020; Gu et al., 2021). Starks et al. (2020) used a full wave model to reproduce the distribution of VLF transmitter wave power over latitude and longitude at the altitude of DEMETER satellites (660 km) and applied this distribution to the ray tracing model to simulate the meridional distribution of the VLF transmitter wave power and wave normal angle in the inner magnetosphere for both ducted and non-ducted propagation. The simulation result shows the latitudinal dependence of transmitter wave power; the wave power is higher at high latitudes than near the equator.

However, due to the orbital limitation of these satellites, there exists a significant knowledge gap in the VLF transmitter emission propagation between the low latitudinal inner magnetosphere and the topside of the ionosphere, and the latitudinal dependence of the transmitter power has not been verified by in situ observations. In this study, we use approximately 3 years (2017-2020) of observations from the Exploration of energization and Radiation in Geospace (ERG/Arase) satellite and 7 years (2012-2019) of observations from Van Allen Probes to statistically study the meridional wave power distribution of NWC VLF transmitter emissions in the inner magnetosphere and analyze the latitudinal dependence of the wave power. Ray tracing simulation is also performed to reproduce the ray paths of the transmitter signal and the corresponding meridional wave power distribution. In Section 2, we give a brief introduction of the ERG/Arase and Van Allen Probes satellites and the data used in this study. In Section 3, we provide two events of NWC transmitter emissions observed by ERG/Arase and Van Allen Probes. In Section 4, the statistical meridional distribution of the NWC wave power is shown as well as that for other ground transmitters. Finally, we use the ray-tracing model to simulate the wave power distribution and compare it with the observational results in Section 5.

The Van Allen Probes (Mauk et al., 2013) consist of two satellites with identical instruments and nearly similar near-equatorial (with an orbital inclination of ∼ 10 ° ) highly elliptical orbits with a perigee of approximately 1.1 Earth radius (R _E ), an apogee of approximately 5.8 R _E and an orbital period of about 537 min. It was launched in August 2012, and its mission ended in October 2019, with an operating duration of approximately 7 years. The high-frequency receiver (HFR) of the Electric and Magnetic Field Instrument Suite and Integrated Science (EMFISIS) (Kletzing et al., 2013) can provide the measurement of wave power spectra of one component of the wave electric field in the plane perpendicular to the spacecraft spin axis with a time resolution of 6 s and 82 logarithmic spacing frequency bins from 10 to 400 kHz.

The Exploration of energization and Radiation in Geospace (ERG) (Miyoshi et al., 2018) satellite with the nickname “Arase” was launched on December 20, 2016, with apogee and perigee of ∼32,000 and 460 km, respectively, and an orbital period of ∼570 min. The High Frequency Analyzer (HFA) (Kumamoto et al., 2018) of the measurements from the Plasma Wave Experiment (PWE) instrument (Kasahara et al., 2018) provides two components electric spectra from a few kHz to 10 MHz (with 480 frequency bins with frequency width from ∼1.2 to ∼98 kHz) every 8 s (one spin period). The inclination of ERG/Arase orbits ( ∼ 31 ° ) is higher than that of the VAP, which can provide middle latitude coverage of VLF transmitter emissions.

To study the latitudinal distribution of transmitter wave power in the inner magnetosphere, we first check how the transmitter signals are captured by ERG/Arase and Van Allen Probes. Figure 1 shows two examples of NWC transmitter signals observed by ERG/Arase (1a) and Van Allen Probes (1b). Both events are selected at the nightside since the wave power of the NWC transmitter is stronger at the nightside in the inner magnetosphere (Zhang et al., 2018). For the ERG/Arase electric field wave power spectrum shown in Figure 1A, we can see strong wave power at frequencies near the frequency of the NWC transmitter signal (19.8 kHz). The frequency channel of 19.53 kHz shows the strongest wave power, and the wave power decreases sharply as the frequency increases or decreases away from 19.8 kHz. The orbit of the ERG/Arase satellite in this event is in the region with magnetic longitude from approximately 157° to 207° and magnetic local time from approximately 21 to 1 h. The magnetic latitude of the 100 km footprints of the satellite in the Southern Hemisphere is from −30° to −41°, which is close to the location of the NWC transmitter (geomagnetic longitude of 187° and geomagnetic latitude of −31°). Thus, we can confirm that the observed strong emission near 19.8 kHz is the signal of the NWC transmitter.

In Figure 1B, similar NWC signals are observed by Van Allen Probes, with a similar satellite footprint region and strong wave power near 19.8 kHz. The power spectra density value observed by Van Allen Probes is weaker than that by ERG/Arase due to the differences between the instruments onboard. The undermeasurement of the electric wave power of Van Allen Probes may be caused by plasma sheath effects (Hartley et al., 2016; Hartley et al., 2017). However, the actual value of the power spectra density does not affect our study of the latitudinal variation of the transmitter wave power since we only focus on the relative change of the wave power rather than the exact wave power value.

From the above cases, we are confident that the wave power of the NWC transmitter can be observed well above the background noise by both ERG/Arase and Van Allen Probes. We can use nearly 3 years of ERG/Arase observations and 7 years of Van Allen Probes observations to perform a statistical analysis of the NWC wave power distribution projected on a meridian plane. The following filters are applied to the data for the analysis: 1. The MLT should be at nightside for stronger emissions than dayside; 2. The wave power of the frequency bin closest to 19.8 kHz should be larger than the mean wave power of the two nearby frequency bins; 3. The magnetic longitude should be close to the longitude of the NWC transmitter (187°). The magnetic longitudinal filter range is selected to be [180°, 195°] (similar to the longitude range used in Zhang et al. (2018)), which is labeled by the two white vertical dashed lines in Figures 1A, B, and most of the strong NWC wave signals are confined within this longitude range.

The NWC signals are obtained by applying the filters mentioned in the previous section to the observations of ERG/Arase and Van Allen Probes satellites. For the ERG/Arase observation, the obtained data are sorted into bins with different geomagnetic latitudes (of latitudinal bin width 2°) and radial distances (of radial bin width 0.2 R _E ) from the center of the Earth. Figures 2A, B show the meridional distribution of the data number and mean value of the NWC wave power. The wave power is obtained by a product of the power spectra density and the width of the frequency bin containing 19.8 kHz. In Figure 2A, most of the data numbers for the bins are high enough to provide reliable statistical mean values. In Figure 2B, we can clearly see a region with strong wave power, illuminating the propagation route of the NWC signals from the Southern Hemisphere to the Northern Hemisphere. The wave power at the high latitude region is much higher (approximately by one order of magnitude) than that near the equator, especially in the Southern Hemisphere. This wave power distribution suggests that as the NWC emission propagates from the Southern Hemisphere to the Northern Hemisphere, the wave power decreases as the wave approaches the equator and then increases as the wave leaves the equator.

To establish the latitudinal dependence of the NWC wave power more directly, we perform further analysis on the meridional wave power distribution. For each latitude bin, we find the peak wave power A ₀ and corresponding peak location r ₀. We plot r ₀ versus latitude in Figure 2B as the black solid line, which indicates propagation path of NWC waves. The magenta solid line is the background magnetic field line calculated by the IGRF-12 model (Thébault et al., 2015) passing through the location of the power peak of the observed NWC signal from DEMETER from Zhang et al. (2018). The variation in A ₀ versus latitude is plotted in Figure 2C, in which we can clearly see that the NWC wave power decreases from nearly 4 × 10⁻⁶ V ²/m ² at approximately -25° latitude to approximately 0.4 × 10⁻⁶ V ²/m ² near the equator and then increases to nearly 1.2 × 10⁻⁶ V ²/m ² at approximately 25° latitude. The apparent sharp decrease near -30° latitude is due to no significant NWC wave power in the region near -30° latitude. The NWC emissions propagate from the ionosphere at approximately -25° into the inner magnetosphere and then propagate only northward.

To verify the latitudinal dependence of NWC wave power based on ERG/Arase measurement, we perform a similar analysis for the Van Allen Probes data (with 1° latitude bin and 0.1 R _E radial distance bin), with corresponding figures shown in Figures 2D–F. The data number in Figure 2D is even higher than in 2a due to the longer operation time and higher spectral data time resolution of the Van Allen Probes. Although the orbital coverage of the Van Allen Probes limits the latitudinal range of the results to only approximately -20°–0°, the latitudinal variation in NWC wave power from the Van Allen Probes observations is similar to that from ERG/Arase observations. In Figure 2F, the variation in NWC wave power decreases from approximately 1.4 × 10⁻⁷ V ²/m ² to approximately 4.6 × 10⁻⁸ V ²/m ², with a decreasing factor of approximately 1/3, which is comparable with the relative decrease from -20° to 0° of the NWC wave power measured by ERG/Arase (∼1/5) shown in Figure 2C. In summary, both ERG/Arase and Van Allen Probes observations confirm that the NWC wave power is minimized near the equator.

The statistical meridional distribution of the NWC wave power is mainly as a result of the propagation of the VLF transmitter emissions in the inner magnetosphere. Using the HOTRAY code (Horne, 1989), we perform ray-tracing simulations to construct the distribution of the NWC wave power and make a comparison against the observed wave power distribution. For the ray-tracing model used, the background magnetic field is the IGRF-12 field (Thébault et al., 2015), and the cold plasma density profile is identical to the typical nighttime profile in Bortnik et al. (2011). Landau damping during the propagation is also accounted for using a suprathermal electron model (Bortnik et al., 2011) with a moderate Kp value (=2). The rays are launched at an altitude of 700 km from −3° to −50° magnetic latitude with a spacing of 0.1°, and with the initial wave normal angle radially outward. In total, 471 rays are traced. The ray tracing simulation provides the paths of these rays on the meridional plane, the wave normal variation along the ray paths, and path-integrated wave damping.

To obtain wave power variation, we make use of the conservation of energy flux along the propagation. That is, when the wave damping is not considered, for the i-th ray with initial latitude λ _i , the dot product of the Poynting vector S and the cross-sectional area A should be conservative: S ⋅A = S ⋅ A = C _i , where S is the magnitude of the Poynting vector, A = d ⋅ r _⊥⋅Δϕ, where d is the cross-sectional distance between the (i − 1)-th and (i + 1)-th rays, r _⊥ is the distance from the Earth’s dipole axis, and Δϕ is an arbitrary fixed azimuthal angle width. Expressing in term of electric field power, the conservation equation can be written as: S P E P E ⋅ d ⋅ r ⊥ ⋅ Δ ϕ = C i (1) where P _E is the electric wave power, and the ratio S P E can be obtained from the dispersion relation along the ray path, and the geometric factor d ⋅ r _⊥ can be obtained by the ray tracing. Thus, the variation in the electric wave power can be expressed as: P E = C i S P E ⋅ d ⋅ r ⊥ Δ ϕ (2) where C i = P E0i ( S P E ) 0 i d 0 i r ⊥ 0 i Δ ϕ is the flux at the starting point of the ray path at λ _i . Introducing a damping factor e − 2 ∫ γ d t where γ _d is Landau damping rate due to the suprathermal electrons (Chen et al., 2012), we can model the electric wave power during wave propagation as: P E = S P E 0 S P E d 0 d r ⊥ 0 r ⊥ e − 2 ∫ γ d d t P E 0 (3) Where the three pairs of brackets on the right side of the equation underline the dispersion, geometric and damping factors, respectively.

The electric wave power at the starting point P _E0 is obtained from the statistical wave power distribution based on DEMETER observations (Zhang et al., 2018), and the entire meridional distribution of the electric wave power can be constructed. Figure 3A shows the simulated meridional electric wave power distribution for the NWC transmitter emission. For each latitude bin, we find the wave power peak A ₀ and the corresponding location r ₀ as we do in Section 4, and the trajectory of r ₀ is shown by the black dashed line. We use the minimum value of peak power A ₀ along the trajectory of r ₀ near the equator as the normalization factor (P _E, norm ∼ 6.9 × 10⁻⁸ V ²/m ²) and the simulated wave power is normalized by P _E, norm in Figure 3A. The magenta solid line is the IGRF magnetic field line starting from the maximum NWC transmitter wave power location observed by DEMETER. The simulated wave power distribution shows similar latitudinal dependence: the wave power decreases as the wave propagates to the equator and then increases as the waves leave the equator. This provides support for the presence of minimal wave power near the equator.

It should be noted that three sets of spectral intensity are used, which are measured by three different instruments from DEMETER, ERG/Arase, and Van Allen Probes. The DEMETER wave spectra that are used as source emission in ray tracing, whose results are compared with ERG/Arase and Van Allen Probes wave spectra. Therefore, the absolute wave powers of the simulation and the observation differ. However, of interest in our study is the relative change in the wave power over the magnetic latitude. For comparison with the simulated distribution, we replot the wave power distribution from the ERG/Arase observation normalized by the minimum value of peak power A ₀ along the trajectory of r ₀ (black solid line) near the equator (P _E, norm ∼ 3.2 × 10⁻⁷ V ²/m ²) so that the equatorial value of the observed wave power matches that of the simulation (Figure 3A). The comparison exhibits similar latitudinal dependence. The relative changes in the wave power between the equator and ± 25° are of the same order ( ∼ 1 0 0.5 to 10¹) from both the observation and simulation results (Figure 3C). In the Southern Hemisphere, the relative change from the ray tracing result is slightly smaller than that from the observation, but in the Northern Hemisphere, the relative change from the simulation is slightly larger. We also plot the variation of normalized wave power peak A ₀ from the Van Allen Probes observation as the green line in Figure 3C, which shows very good agreement with the simulation result. Overall, the relative wave power changes from the simulation and observation are comparable.

In addition to the NWC transmitter, we also analyze the wave power distribution for three other VLF transmitters, i.e., NAA, NLK and NML transmitters. Their respective wave frequencies are 24, 24.8 and 25.2 kHz, and their magnetic longitudinal ranges are [−5°, 16°], [294°, 304°] and [322°, 332°], respectively (similar to the longitude ranges used in Zhang et al. (2018)). The results for these three transmitters are plotted in Figure 4, in a format similar to Figure 3. The normalization factors P _E, norm for the NAA, NLK and NML transmitters are 5.2 × 10⁻⁹ V ²/m ², 4.8 × 10⁻⁹ V ²/m ² and 4.4 × 10⁻⁹ V ²/m ² respectively for the simulation result and are 4.4 × 10⁻⁸ V ²/m ², 1.1 × 10⁻⁸ V ²/m ², 1.5 × 10⁻⁸ V ²/m ² respectively for the ERG/Arase observation result. We do not show the results from Van Allen Probes observation because of their limited latitude range. Different from the NWC transmitter, the locations of these three transmitters are in the Northern Hemisphere and at relatively high latitudes ( ∼ 54 ° ) . Thus, the waves can propagate to higher L-shell regions. The meridional wave power distributions of these three transmitters all show a minimum near the equator: the wave power decreases first when propagating to the equator and then increases when the waves leaves the equator. The simulation results of these three transmitter wave power match very well with the observational results for both the normalized value of the wave power and the trajectories of the peak power. Additionally, the relative changes in the wave power from the high latitude region (±30°) to the equator minimum are all on the order of ∼ 1 0 1 to 10² for both the observation and simulation results.

To investigate the cause of the existence of VLF transmitter minimum wave power near the equator, we check the contribution of the three factors in Eq. 3, the dispersion, geometric, and damping factors. For the NWC transmitter, the minimum of the wave power near the equator is in part due to the variation of the dispersion factor but mainly determined by the geometric factor. The geometric factor, which is inversely proportional to the cross-sectional area, decreases when the wave propagates to the equator due to the spread of the rays, and then increases in the high latitude region due to the convergence of the rays (Buckley and Budden, 1982; Hanzelka and Santolík, 2019; Starks et al., 2020). The damping effect is rather weak at such a low L shell region due to the lack of the energetic electrons responsible for Landau damping (Li et al., 2010). For the three high latitude transmitters (NAA, NLK, and NML), both the dispersion and geometric factors play significant roles in forming the wave power minimum near the equator. Moreover, compared with NWC, the emission from these high latitude transmitters can propagate through the higher L shell region. As a result, the Landau damping factor becomes non-negligible and reduces the wave power along the propagation path, which results in the shift of the latitude of the wave power minimum farther away from the transmitter source region.

In this study, we use approximately 3 years of observations from the ERG/Arase satellite to statistically study the meridional distribution of the NWC transmitter wave power from the equator to the middle latitude region. The result shows significant latitudinal dependence of the wave power: the wave power is minimized near the equator and increases as the latitude increases. This latitudinal dependence is also confirmed by the statistical results of the 7-year observation from Van Allen Probes, although over a narrower latitude range (∼—20° to ∼ 0 ° ). Then, we use the HOTRAY code to perform ray-tracing simulation of the transmitter signals and reproduce the meridional distribution of the NWC transmitter wave power. The simulation results coincide with the observation very well, especially on the relative variation of the wave power over magnetic latitude and on the presence of transmitter wave power minimum near the equator. Similar statistical observational analysis and ray-tracing simulation are also performed for NAA, NLK and NML transmitters and the results also confirm the presence of wave power minimum near the equator. Although the transmitter wave amplitude minimum near the equator region is also noted in Starks et al. (2008) before, the altitudes of the satellites used in their study are nearly fixed and thus, we cannot confirm whether the observed minimum near the equator is caused by the variation in the transmitter power over different L-shells or by the variation in the power along the magnetic field line of the same L-shell. Our results demonstrate that the transmitter wave power has a strong dependence on the latitude and is minimized near the equator.

The comparison between the observed and simulated wave power distributions also provides some additional information about the propagation mode of the transmitter signal. For the three higher latitude transmitters, the simulated wave power distribution matches the observed distribution very well, which strongly suggests that the dominant propagation mode is the non-ducted mode in the region of L < 3. For the NWC transmitter, in addition to the strong wave power in the non-ducted propagating region, considerable wave power exists in the lower altitude region, which may associate with the ducted and pro-longitudinal (PL) mode (Thomson and Dowden, 1977; Starks et al., 2001). The PL mode usually occurs at lower altitude regions (for low latitude transmitters such as NWC) and is caused by the sharp latitudinal plasma density gradient in the topside ionosphere. As a final remark, the wave power distribution of different VLF transmitters on the meridian plane reproduced in our study not only reveal the physics of wave propagation, but also can be applied to estimate the energetic electron loss caused by VLF transmitter emissions in the inner magnetosphere. The latter will be left as future work.