Complex Whistler‐Mode Wave Features Created by a High Density Plasma Duct in the Magnetosphere

A Van Allen Probes observation of a high‐density duct alongside whistler‐mode wave activity shows several distinctive characteristics: (a)—within the duct, the wave normal angles (WNA) are close to zero and the waves have relatively large amplitudes, this is expected from the classic conceptualization of ducts. (b)—at L‐shells higher than the duct's location a large “shadow” is present over an extended region that is larger than the duct itself, and (c)—the WNA on the earthward edge of the duct is considerably higher than expected. Using ray‐tracing simulations it is shown that rays fall into three categories: (a) ducted (trapped and amplified), (b) reflected (scattered to resonance cone and damped), and (c) free (non‐ducted). The combined macroscopic effect of all these ray trajectories reproduce the aforementioned features in the spacecraft observation.


Introduction
Whistler mode waves in the magnetosphere are well known to have a significant role in the near-Earth space physics and space weather processes (Helliwell, 1969;Millan & Thorne, 2007;Shprits et al., 2008).Naturally occurring whistler mode waves such as chorus, hiss and lightning-whistlers play a major role in energetic particle losses and acceleration within the Earth's radiation belts (Horne, 2007;Millan & Thorne, 2007;Tsurutani et al., 2009).Global radiation belt models require knowledge of the spatial distribution and properties of these waves for accurate predictions of particle fluxes (Albert & Young, 2005;Glauert & Horne, 2005).However, the propagation of waves in the magnetosphere is complex due to the inherent inhomogeneity and anisotropy (Maxworth & Gołkowski, 2017;Streltsov & Goyal, 2021).As such, understanding relevant features of whistlermode wave propagation remains an active area of space physics research.
Of particular scientific interest are the impacts of cold plasma density structures known as "ducts" on the propagation of whistler-mode waves.Ducts are enhancements (or depletions) of background cold plasma density that conform to the Earth's magnetic field lines; they can effectively guide waves over thousands of kilometers in a manner akin to an optical fiber (Karpman & Kaufman, 1982;Walker, 1976).Ducts were first hypothesized in Storey (1953) as means of explaining ground observations of lightning-induced whistlers.Since then, a large body of research has studied the theoretical guiding properties of ducts (Karpman & Kaufman, 1982;Ke et al., 2021;Scarabucci & Smith, 1971;Smith et al., 1960;Sonwalkar, 2006;Streltsov et al., 2006).Although the theoretical conditions required for ducted propagation are now well understood, the majority of prior evidence for their existence has been indirect inference from observed wave features (Angerami, 1970;Carpenter, 1988;Hayakawa & Iwai, 1975;Ondoh, 1976).Although a handful of spacecraft observations were published in the 1990s (such as Sonwalkar et al., 1994), direct in-situ observations of magnetospheric ducts and associated waves with high resolution has only been feasible within the last two decades of spacecraft missions (VAP, MMS, etc.).
In recent years, researchers have directly observed several examples of ducts and cold density structures and have examined properties of wave guiding and focusing (R. Chen et al., 2021;Hosseini et al., 2021;Streltsov & Goyal, 2021;Williams & Streltsov, 2021).For instance (Streltsov & Bengtson, 2020), used measurements of high density ducts from VAP alongside simulations to confirm the guiding ability of ducts with observed magnetospheric parameters.Their analysis, however, neglected the impact of curvature of the geomagnetic field and primarily focused on waves within the duct.Another recent study by Ke et al. (2021) also utilized an Van Allen Probes observation of a duct, however, the density profile was coupled to a 2D particle-in-cell simulation with curvature included.The results showed that ducted waves remained quasi-field aligned while undergoing an increase field strength via focusing and gyro-resonant amplification.Such recent studies have successfully confirmed the long-expected guiding properties of realistic magnetospheric ducts.However, the distribution of waves external to ducts have mostly been ignored in prior work.It has thus been somewhat unclear as to the combined impact of ducting and external scattering on the wave distribution outside the duct.
Using an in-situ observation on the Van Allen Probe A (VAP-A), we analyze a ∼400 km sized high density duct and the associated wave fields.The observations are consistent with a large equatorial source of lower band whistler mode waves.It is found that, in addition to typical ducting behavior (i.e., parallel propagation and high amplitude inside the duct), several new features are present.The wave fields are shown to have unexpectedly high wave normal angles (WNA) at the left edge (earthward) of the duct and a "shadow" region of diminished wave power on the right (non-Earthward) side that is wider than the duct itself.Using raytracing simulations, with inclusion of cyclotron growth and Landau damping, all the salient features of the observations are readily explained as a combination of propagation and hot plasma effects.Section 2 describes the relevant features of the observations, and Section 3 discusses raytracing simulation results that match observations.

Spacecraft (Van Allen Probes) Observation of a Large High Density Structure
We examine a time interval from 12:00 to 14:00 UT on 28 May 2017 from the Radiation Belt Storm Probes A (VAP-A, operated from 2012 to 2019) measurements.The data is collected during the early recovery phase of the geomagnetic storm of May 28.The storm main onset was around 07:00 UT on 28 May with Dst index 140 nT and the average AE of ∼1,000 nT from 00:00 to 07:00 and maximum value of approximately 1,600 nT at ∼02:00 UT.VAP-A was located in the local afternoon sector MLT = 14.33, at L ≈ 4.5 4.6, at magnetic latitude λ ≈ 10.6°.The plasmapause was at L ≈ 2 (more details about this geomagnetic storm can be found in Da Silva et al. ( 2022)).
Of particular interest in this study is the sharp density enhancement from 25-30 to 200-250 cm 3 which was observed from 13:10:25 to 13:12:35 when VAP-A crossed the region with the plasmasphere plasma parameters.The spatial scale of this structure was found to be ∼300-400 km across the background magnetic field direction; density profile is shown in Figure 1a.Presumably this structure was a duct-like (following the latitudinal structure of the background magnetic field) product of the plasmapause erosion during the geomagnetic storm.Hiss-like whistler-mode wave activity is also simultaneously observed in the 0.2-2 kHz frequency range during this time interval.
The observation is rich with several features, however, we focus on specific characteristics for this study: 1. Duct Interior: The wave magnetic field amplitude inside the duct is shown to be up to 4× higher than the average field strength outside (∼13:10:30-13:12:30). Furthermore, the WNA tends to stay under ∼10°which implies approximately field aligned propagation.This is shown on Figures 1b and 1c. Figure 1d shows the azimuthal angle of the wavevector relative to the geomagnetic field.The variability shows some threedimensional effects within the duct, however, this is not directly investigated for this study.2. Shadow Region: The wave magnetic field amplitude on the right side of the duct (∼13:12:30-13:15) is significantly depleted (essentially at the noise floor).This "shadow" region is larger than the size of the duct itself.The WNA in this region is seemingly high at approximately 40°, however, since the wave fields are at noise levels, this feature is not considered in this study and is assumed to be due to quasi-electrostatic noise.This is shown on Figure 1b.
Unlike the high WNA in the shadow region, the wave fields on the earthward side have relatively high amplitudes.This is shown on Figure 1c.The relatively high value of planarity (above ∼0.5 for amplitudes that are large enough) shown in Figure 1e justifies the validity of evaluating the WNA even at these high levels (at least near 1 kHz which is the frequency utilized in the simulations of Section 3). 4. Free Region: The WNA far from the duct's edge (13:07:40-13:10) is in the range of 20-40°.This is expected from a ray-optics interpretation (∼4°W NA increase per degree latitude) for waves that do not interact with a duct, and is thus "free" of the duct's effects.This is shown on Figure 1c.
The features in the interior of the duct and the free region are expected from conventional ducting theory.However, the large shadow region and oblique WNA at the duct's edge do not have a straightforward interpretation.Given the relative size of these effects relative to the duct (particularly the shadow region), a deeper physical explanation is warranted.As such, we utilize a numerical approach to model the observed event; a description of the code and simulated results is provided in Section 3.

"DenRay" Model Description
In order to model the event presented in Section 2, we utilize a raytracing simulation (DenRay code) that is described herein.The geomagnetic field, B 0 , is modeled as an Earth-centered dipole.The background "cold" plasma is described by a static electron density, N c (r), that is a function of position, r.For latitudes (<12°) and frequencies of consideration (200 Hz-2 kHz), the dynamics of ions are ignored and simply play the role of a neutralizing positively charged background.Lastly, the energetic (source) electron population is described by a phase-space distribution function, f(v ∥ , v ⊥ ), where v ∥ and v ⊥ are the components of velocities parallel and perpendicular to the geomagnetic field.The hot plasma density, N h , can then be determined by integrating the distribution over all velocities.
In the magnetosphere, the cold electron density is typically orders of magnitude higher than the energetic electron density.With the aforementioned assumption (i.e., N c ≫ N h ), the propagation of whistler mode waves are exclusively governed by the cold plasma electrons.Energetic electrons play the role of growth or damping of waves via gyro-resonant instabilities, however, the cold population plays a concurrent role in defining the resonance conditions (Brinca, 1972).
In order to the model the propagation of waves, DenRay utilizes a raytracing formalism.Raytracing uses an Eikonal equation (Haselgrove, 1954) approach to track the propagation of wavefronts (rays) in a computationally efficient manner (R. B. Horne, 1989;Xie et al., 2022).This methodology is justified for a background density and magnetic field that does not vary appreciably over a wavelength and is a common approach in whistler-mode wave modeling in the magnetosphere (Abel & Thorne, 1998;Bortnik et al., 2007;L. Chen et al., 2009;Inan & Bell, 1977).The code tracks the position of a ray, r, and its corresponding wavevector, k, through time, t, via the raytracing equations, The quantity ω represents the wave's angular frequency and v g is the group velocity.D(ω, k, r) is the dispersion function of the cold plasma which is given by, The quantity μ 2 = |k| 2 c 2 / ω 2 is the index of refraction.A, B, and C are cold plasma parameters (Stix, 1992) and are reproduced in Appendix A. Given any initial condition for k and r (that satisfies D(ω, k, r) = 0), Equations 1 and 2 can be integrated forward in time to determine the temporal evolution of rays in the background cold magnetized plasma.
The relevance of energetic electrons becomes apparent when determining the growth or damping of rays via kinetic gyro-resonant instabilities (Brinca, 1972).For this study, the energetic electron distribution function is taken to be a Bi-Maxwellian, with parallel and perpendicular thermal velocities, v ∥ T and v ⊥ , respectively.The expression for the distribution function is shown in Appendix A.
The growth or damping waves from can be determined by linearizing the Vlasov equation and enforcing the condition, N h ≪ N c (Kennel, 1966).The temporal growth rate, γ L , of whistler mode waves is then given by, The quantity θ is the wave normal angle and represents the angle between the geomagnetic field and the wavevector.The expressions for I 1 , I 2 , and I 3 are given by in Appendix A.
In this work, the nth ray is modeled as a localized "power packet" with power P n (r(t)) in Watts.Thus, the impact of growth and damping is accounted for using, Finally, the magnitude of the power density W m 2 ) is estimated numerically using a large collection of rays.A spatial grid is first constructed over the simulation domain and the initial power of each ray, P n (r(t = 0)), is assigned based on the initial power density at the source.Specifically, each ray's power is given by the power density multiplied by the surface area of the cell and then divided by the number of rays in the source cell.Rays are then traced continuously through time and as each ray enters a cell in the grid, the total power entering the cell is accordingly accumulated, where the modified power in each ray is computed using Equation 5.The power density in each cell is then determined by dividing the total power per cell by the surface area of the cell.Furthermore, the local WNA in each cell is computed using nearest neighbor interpolation as long as the power density in that cell in nonzero.It is worth noting that this methodology is different from most other raytracing codes that typically assume that the power density along each ray can only change via damping or growth.However, the approach in DenRay also includes changes in power density due to the focusing or spreading of rays, which is important in the case of ducts.
The DenRay code thus allows efficient calculation of power density and WNA as a function of position and time.The results are expected to be reasonable as long as the growth of waves is not drastic enough to warrant nonlinear effects.As such, Section 3.2 shows simulation results of the event described in Section 2 using DenRay.

Simulation Results
In order to model the duct observation described in Section 2, a few conditions must be assumed.First, the source of waves is assumed to occur on the equatorial plane, which is well-justified by observations (LeDocq et al., 1998).Second, the wavevectors at the source plane are assumed to be field aligned, which is also expected for some classes of whistler-mode waves (Goldstein & Tsurutani, 1984).Third, the spatial extent of the wave source region is assumed to be much larger than the duct itself.Although a non-uniform distribution of wave power is possible (and may also explain some of the observed features, such as the shadow region), it has not been considered as part of this study.Fourth, the density profile of the duct is assumed to be field aligned (i.e., constant along constant an L-shell) which is a typical assumption of magnetospheric ducts.Fifth, the source of waves is assumed to emit continuously from the equator, as such DenRay is run in "continuous wave" mode so only the spatial distribution (and not temporal effects) is considered for this study.Finally, the simulation is run in only two dimensions (x, z) for computational efficiency and ease of visualization.The 2D nature of the simulation implicitly assumes that there is no variation in density in the y-direction.As such, any azimuthal effects of density variation or propagation are ignored in this simulation study.
It is also worth noting that the density gradients are relatively high within the duct, however, the high values of electron density (150-200 el/ cm 3 ) still ensure that the small-scale variations are several wavelengths in size which justifies the use of raytracing.However, the edges of the duct have spatial scales that are on the order of a wavelength or larger.As such, some edge effects (such as diffraction and boundary reflection) may not be adequately captured by raytracing and can be investigated with full-wave modeling in future work.
The simulation of the event assumes that the duct is located at 4.83 < L < 4.88.The density profile at the observation latitude is directly taken from the data and the field-aligned assumption is then enforced.This result of this process can be seen in Figures 2a and 2b.10,000 rays are launched from the equatorial plane (z = 0) and are uniformly distributed in the range 4.5 ≤ L ≤ 5.The wave frequency is set to 1 kHz which is approximately in the middle of the band seen in the observation.The energetic electron population has parameters given by v ∥T = 0.05c,v ⊥T = 0.055c, and N h = 0.25el/ cm 3 , which is a weakly anisotropic distribution.These values are chosen based on the range of "typical" parameters of source electrons in the magnetosphere.
Prior to examining the spatial distribution of power density and WNA, it is first useful to understand the general types of ray trajectories in the system.Figure 2a shows 20 ray trajectories in white (out of the 10,000 in the full simulations) superimposed on the background density profile.Upon visual examination, it can be seen that the rays each fall into three basic categories.As such, representative rays in each category are shown in color in Figure 2a.Additionally, the relative power along each of these representative rays along with the WNA as a function of magnetic latitude is shown on Figures 2c and 2d respectively.The three categories are as follows: 1. Free Rays: These rays are far from the left edge of the duct (green) or on the right side of the duct (cyan).The free rays do not experience any change in trajectory due to the presence of the duct and propagate as they would in a smooth background plasma.As shown in Figure 2d the WNA of these rays increases approximately linearly with latitude.This in turn leads to increased Landau damping which is apparent in the reduction of power with latitude as shown in Figure 2c.2. Reflected Rays: These rays are launched close enough to the left edge of the duct that they "collide" with the duct within the simulation domain and reflect from it (red); the "reflection" occurs since the WNA exceeds the Gendrin angle upon encountering the duct and thus switches direction.Upon reflection, the rays abruptly increase in WNA and saturate at the resonance cone angle (≈82°); this phenomena can be clearly seen in Figure 2d.Due to the sudden change in WNA, the reflected rays experience significant Landau damping and rapidly decay once reflection occurs, as seen in Figure 2c.Although these rays do quickly lose power from the reflection process, they also tend to "bunch-up" at the left edge of the duct at higher latitudes as seen in Figure 2a.Thus the power density of reflected rays is determined by a balance of strong Landau damping and strong bunching (focusing) or rays.3. Ducted Rays: These rays are guided within the duct and stay trapped over the entire extent of the simulation (black ray); this is the classic conceptualization of ducting that has been discussed in the literature for over 60 years (Smith et al., 1960).As seen in Figure 2d, the WNA of the ducted ray stays close to zero which enforces approximately field aligned propagation during the entire simulation.Furthermore, the low range of WNA allows for growth via gyro-resonance to outcompete Landau damping (Hanzelka & Santolík, 2019).This is evident in Figure 2c where the ducted ray is shown to increase in power as it propagates away from the equator.
In addition to the general types of ray trajectories described above, Figure 2a also shows the emergence of the "shadow" region.Specifically, rays that are trapped (ducted) by the duct would have otherwise followed the free ray trajectories.The region where these rays would have been is now devoid of wave power.It is worth noting that this only the case because the duct follows the topology of the curved geomagnetic field and would not be the case in straight duct.This is a defining feature of curve ducts that has not been noted in prior (to the best of the authors' knowledge).
The ray trajectories shown in Figure 2 provide a useful framework for interpretation, however, a quantitative comparison to data requires computing the spatial distribution of power density and WNA.Figures 3a and 3b show the power density (normalized to the source value) and WNA distribution using 10,000 rays and the method described in Section 3.1.The dashed black line in Figures 3a and 3b represents the spacecraft trajectory (λ ≈ 10°).
Figure 3c shows both the power density and WNA as seen by the spacecraft (with the normalized electron density superimposed) along with blue dotted lines representing the approximate duct boundaries.Each of the major regions described in Section 2 are also labeled on the simulation results shown in Figure 3c.It is first worth discussing the presence of expected features both inside and outside the duct.Within the duct itself, the power density increases with latitude while the WNA tends to stay under 20°(with some exceptions).This is consistent with discussion of the ducted ray category described previously and the typical conceptualization of the ducting phenomenon.That is, the ducted waves tend to have lower WNA and hence get amplified via the gyro-resonant instability.In contrast, far enough from the duct, the power density is approximately 40% of the value inside the duct.Interestingly, the power density stays approximately constant everywhere outside the duct even though Landau damping does occur.This is because the rays tend to converge as latitude increases and this focusing effect counterbalances the small loss from Landau damping.Additionally, the WNA of these waves increase approximately linearly with latitude as expected.These aforementioned features are what is to be expected from ducted and unducted waves from basic theory (Bortnik et al., 2007(Bortnik et al., , 2011;;Breuillard et al., 2012) and known from spacecraft statistics (Agapitov et al., 2012(Agapitov et al., , 2013(Agapitov et al., , 2017)).
The first new finding shown in Figure 3a is the large shadow region to the right of the duct.This shadow region is also present in the VAP observation shown in Figure 1b and is correspondingly reproduced in the simulation in Figure 3c.As previously discussed, the shadow is caused by the trapping of rays by the duct which leaves behind a large depletion of power.A noteworthy consideration is that the area of shadow region is considerably larger than the duct itself, and this area increases rapidly with latitude.This phenomena is only present due to the duct's curvature and would not be the case for a straight duct that is often considered in theoretical studies (such as A. V. Streltsov et al. (2006)).As such, this is an important finding that should be considered when modeling the distribution of wave power in the presence of ducts.
The second new finding is the increase in WNA at the left edge of the duct with a spatial extent that is on the order of the size of the duct itself.This feature is also present in the VAP observation shown in Figures 1c and 3c (green curve representing the average between 500 Hz and 1 kHz) and is also seen in the simulation (red) in Figure 3c.
The origin of this phenomena is associated with the reflected rays discussed previously.The rays that reflect off the left edge of the duct are pushed to the resonance cone angle, which explains the acute increase at the edge of the duct.Furthermore, since all the reflected rays propagate at approximately the same angle, they tend to bunch together and cause an increase in amplitude by almost an order of magnitude.However, these rays also undergo significant Landau damping due to their high WNA which counteracts the focusing effect.This complex extraducting effect has not previously been seen in the literature.It is worth noting that although the observations and model are qualitatively similar, the exact value of WNA differs between the two due to the generic choice of hot plasma parameters in the model (not directly from the observed event).As such, the weighting due to wave power effects in computing the WNA is not one-to-one between observations and simulations.Regardless, this local and rapid increase in WNA near the duct's edge is present in both observations and simulations and can have a considerable impact on subsequent wave-particle-interactions.As such, this change in WNA may need to be included in radiation belt models for conditions where ducts are present.

Summary and Discussion
We have analyzed a VAP-A observation of a ∼400 km wide duct alongside the corresponding whistler mode waves just outside the plasmasphere.Several features are observed that are explained using raytracing simulations: 1. Traditional Ducting: Waves within the duct are approximately field aligned and have a wave amplitude that is higher than outside the duct.Raytracing simulations show that the duct traps rays in a snake-like pattern which enforces a low WNA.At the same time, the low WNA is conducive for growth via first order gyroresonance.Accordingly, the waves remain mostly field aligned and are amplified in the process.This result is expected from the general understanding of ducting physics and corroborates substantial prior research.The accumulating of wave power inside the ducts presumably explains the enhanced scattering at high latitudes leading to bursty precipitations of subrelativistic electrons observed as microbursts (Breneman et al., 2017;Mozer et al., 2018, p. 201) accordingly providing a spatial constrain to the microburst scales (Shumko et al., 2020), which are generally smaller than the chorus source region (Agapitov et al., 2017(Agapitov et al., , 2021)).Also, the systems of the dense duct related to the plume region are favorable for wave power penetration into the plasmasphere (Agapitov et al., 2018).2. Shadow Region: A large shadow region of depleted wave power larger than duct itself is seen on the right side of the duct.Using raytracing, it is shown that the shadow is formed by the trapping of rays in a curved duct.Ducted rays follow the magnetic field line of the duct and hence leave behind a large gap in wave power.This phenomena has not been reported in prior work, and demonstrates a relatively drastic feature of a realistic duct.
Journal of Geophysical Research: Space Physics 10.1029/2023JA032047 HARID ET AL.

Resonance Cone Reflection:
The WNA on the left side of the duct is highly oblique for a few hundred kilometers away from the duct.Using raytracing, it is shown that a large family of rays outside the duct can collide and reflect from the duct's left outer boundary.The reflection process forces these rays to be close to the resonance cone.The high WNA force these rays to quickly damp (via Landau damping) as they propagate away, this process forces the high WNA region to be localized around the edge of the duct.
The aforementioned properties demonstrate the complexity of whistler wave dynamics in and around a realistic high density duct in the magnetosphere.Although several other features remain to be explored (frequency dependence, subwavelength scattering, nonlinearity, ion effects, etc.), this work presents a valuable step forward in the understanding of whistler-mode wave propagation in the presence of magnetospheric density structures.

Appendix A: Expressions for Plasma Quantities
The quantities A, B, and C (Stix, 1992) is given by, The quantities P, R, L, and S in Equation A1 are given by, where ω c = ⃒ ⃒ ⃒ q e B 0 m e ⃒ ⃒ ⃒ and ω p = ̅̅̅̅̅̅̅ N c q 2 e ϵ 0 m e √ represent the angular gyrofrequency and plasma frequency respectively.The quantities q e , m e , and ϵ 0 are the electron charge, electron mass, and free space permittivity respectively.
The Bi-Maxwellian distribution is given by, The quantities v ∥ T and v ⊥ T represent the thermal parallel and perpendicular velocities.
The quantities I 1 ,I 2 and I 3 are given by Brinca (1972),

Figure 1 .
Figure 1.(a) The electron density from Van Allen Probe A EMFISIS measurements on 28 May 2017 (time in UTC).(b) The magnetic field dynamic spectra.(c) The wave normal angle distribution.(d) Azimuthal angle around the background magnetic field.(e) Planarity of the wave field.

Figure 2 .
Figure 2. (a) Twenty sample ray trajectories (white) and four representative trajectories (colored) superimposed on the cold plasma density background.(b) Equatorial profile of the duct based on VAP data.(c) Relative power along each ray for the four representative trajectories, and (d) corresponding wave normal angles for the four trajectories.

Figure 3 .
Figure 3. (a) Relative power density distribution, (b) wave normal angles (WNA) distribution, and (c) modeled WNA and power density (with normalized electron density and WNA from observations superimposed) as seen by VAP.

Journal of Geophysical Research: Space Physics
HARID ET AL.