We solve Maxwell’s curl equations together with the equation of motion (Lorentz force equation) for a cold plasma fluid. The system of equations can be written as ∇ × B = ∑ i N p μ 0 J s i + J i + 1 c 2 ∂ E ∂ t , (1) ∇ × E = − ∂ B ∂ t , (2) ∂ μ 0 J i ∂ t = ω p i 2 c 2 E + μ 0 J i × Ω i . (3) Here we use the SI system of units, with c standing for the speed of light and μ ₀ for vacuum permeability. E and B represent the dynamic electric and magnetic fields, while the static dipole field B ₀ is implicitly included through the signed vector gyrofrequency Ω i = q i B 0 m i , (4) where the index i stands for particle species in the cold plasma: electrons and the three ions H+, He+, and O+ (number of particle species N _p = 4). The assumption | B ₀ |≫| B | was used in Eq. 3 to remove the nonlinear term J × B . The number density n _i of individual species is introduced through the plasma frequency ω p i = n i q i 2 ε 0 m i . (5) The symbols q _i , m _i , and ɛ ₀ stand for the particle charge, mass, and vacuum permittivity, respectively. J _i stands for current density associated with the fields E and B , and J s i is the external current source. For the purposes of EMIC wave propagation, the electron mass can be replaced by m _e,num > m _e to gain a major computational performance boost at the expense of minor inaccuracies in wave propagation properties near the resonance cone—see the Supplementary Material for a short discussion of these inaccuracies.

The equations are solved by implementing the finite-difference time-domain (FDTD) method described by Pokhrel et al. (2018), where the Boris method is used to evolve the current J _i . In the case of ion cyclotron waves, the time step Δt in the numerical solution of the Maxwell equations is typically much larger than the time step Δt _c required by the Boris algorithm to advance the Lorentz equation for current with good precision. The singular update method from Pokhrel et al. (2018) provides matrix equations which advance J by MΔt _c = Δt in a single step, where M is an integer. A 2D version of the staggered Yee grid is used (Yee, 1966), with the current being placed symmetrically at the center of the Yee cell.

The minimum required size of each cell is determined by wavelength λ _w of the modes supported by the plasma fluid confined in the chosen simulations box and must be generally tested before each run. The time step Δt must not violate the Courant-Friedrichs-Lewy (CFL) condition Δ t < Δ t C F L = 1 / c 1 / Δ z 2 + 1 / Δ x 2 , (6) where Δz and Δx are the dimensions of a grid cell (Gedney, 2011). The CFL stability condition is necessary but not sufficient, so an even lower time step might be required, depending on the choice of Δz and Δx. The grid is located in the meridional plane and represented in Cartesian coordinates, with z aligned with the dipole axis and x pointing away from the Earth. The origin (z, x) = (0, 0) sits in the center of the source region.

Similarly to Hosseini et al. (2021) and Hanzelka and Santolík (2022), we use a 1D current density source located along z = 0 with a finite halfwidth w _J. The current amplitude distribution is described by the shape function G x = cos 2 π x 2 w J i f | x | < w J (7) G x = 0 o t h e r w i s e . (8) Time-dependence of the current amplitude is incorporated in the function T t = T t t T m t i f 0 < t < t max (9) T t = 0 o t h e r w i s e , (10) where T _t is a tapering function T t = cos 2 π t − t r a m p 2 t r a m p i f t < t r a m p (11) T t = 1 i f t r a m p ≤ t ≤ t max − t r a m p (12) T t = cos 2 π t − t max + t r a m p 2 t r a m p i f t > t max − t r a m p , (13) and T _m represents subpacket modulations. Subpackets are defined by N _s − 1 local modulation minima and N _s local maxima, with the adjacent minima and maxima being connected by cos² functions, similarly to the tapering function above.

The components of the source current must be obtained from the cold plasma dispersion relation. Considering that the EMIC waves are generated through the first-order cyclotron resonance with ions that propagate along the field line, it is reasonable to assume that the current should be near-circularly polarized, producing quasiparallel electromagnetic emissions. Nevertheless, the simulation code is not limited to any specific wave mode and allows for an arbitrary value of θ _k, so the full 3D field description is given in the following. We start by defining time-dependent frequency ω(t) and wave normal angle θ _k(t). As the next step, the refractive index μ and the complex conductivity tensor σ _i are obtained (e.g., Gurnett and Bhattacharjee (2017), Chapter 4). We continue by calculating the complex electric field with E _x normalized to −1j, where j is the imaginary unit. The corresponding magnetic field is obtained from Faraday’s law in the Fourier space, B x = − k z E y / ω (14) B y = k z E x − k x E z / ω (15) B z = − k x E y / ω , (16) where the azimuthal angle was set to zero, and thus k _y = 0. Since the calculations are performed at the equator, the solar magnetic coordinates used for the Cartesian grid coincide with field-aligned coordinates. The fields are then renormalized so that |B| = 1, and after that multiplied by G(x)T(t)B _w0, where B _w0 is the desired equatorial peak amplitude of EMIC wave field. Having the properly normalized electric field, we employ the generalized Ohm’s law and obtain the current J i = σ E . (17) As a final step, we need to consider the effects of discretization. Because the current source behaves as a delta function in the z direction, but the grid cells are finite, an additional numerical factor must enter into the initialization of our simulation. This factor g(ℓ) depends purely on ℓ = λ _w/Δz and has to be determined numerically. We ran simulations with ℓ going from 8 to 24 and found the least-squares power law fit g ℓ = 2.201 ℓ − 0.891 (18) with a coefficient of determination R ² = 0.9991.

The source current vector can be formally written as (dropping the particle index i) J s = J x ⁡ cos ⁡ ψ x , t , − J y ⁡ sin ⁡ ψ x , t , J z ⁡ cos ⁡ ψ x , t , (19) where the J _x , J _y , and J _z components are obtained from Eq. 17 with the correction introduced in Eq. 18, and ψ x , t = ψ 0 − ∫ − w J x d x ̃ k x x ̃ , t + ∫ 0 t d t ̃ ω ′ t ̃ (20) is the harmonic phase (ψ ₀ chosen such that ψ(0, 0) = 0). Assuming a constant θ _k and constant frequency, we get ∂ ψ ∂ x = − k x (21) ∂ ψ ∂ t = ω . (22) However, in a general case, the frequencies appearing in (20) and 22 are different. With a constant wave normal angle θ _k (i.e., time-independent k _x ) an a chirp rate r _c = (ω ₁ − ω ₀)/t _max, setting ω′(t) = ω ₀ + r _c t will result in ω = ω′ in Eq. 22. When k _x (t) includes explicit time dependence, it will enter the time derivative of phase and introduce a spatial gradient of frequency: ω = ω′ − x∂k _x (t)/∂t. For simplicity, we will further consider only sources with a fixed k -vector direction.

Wave reflection at the simulation box boundaries is mitigated by implementing damping regions. Following Umeda et al. (2001), we introduce a masking factor f _m ≤ 1 and multiply the electromagnetic fields by this factor in each time step. At distances from the boundary larger than d _damp, we set f _m = 1. In the damping region near box boundaries, f _m decreases parabolically. The optimal rate at which the masking factor should decrease depends on the time step and the phase velocity. The efficiency of damping improves with increasing d _damp, which is usually chosen as a small integer multiple of a typical wavelength. Due to the variability of wave vectors and frequencies in our simulations, there is no simple method that would prescribe the best choice of damping parameters. We therefore numerically tested the efficiency of damping to ensure that the power of reflected waves in results presented in Section 3 is less than 10^−3.5 of the incident waves.

The cold electron density distribution is based on the empirical model from Denton et al. (2002) adapted for dipole. The latitudinal dependence is described by n e = n e 0 cos ⁡ λ 2 α , (23) where n _e0 is the equatorial profile, and the formula for the exponential factor in SI base units reads α = 36.08 − 6.36 ⁡ log 10 n e 0 + 0.28 log 10 2 n e 0 − 0.43 L . (24) The equatorial profile uses the best-fit power law for the plasmatrough (Denton et al., 2004) n e 0 = 10.13 − 4.1 ⁡ log 10 L (25) rescaled to coincide with the prescribed ω _pe0 = ω _pe (z = 0, x = 0). The rescaled model does not represent the best fit, but the radial density gradient remains within the range predicted by the various density models discussed by Denton et al. (2004).

The cold electron density in our model can be further modified by density crests (increase) and troughs (decrease) with a Gaussian radial profile, which can be used to simulate ducted propagation. Symbolically, n _e,tot = n _e n _d, where n d = 1 + ∑ j = 1 N d δ n j e − L − L d j 2 2 σ L j . (26) The duct parameters are: number of ducts N _d, relative density change δn, central L-shell of the duct L _d, and characteristic width σ _L. As noted by Hanzelka and Santolík (2022), these ducts are different from the commonly assumed 3D tubes (Angerami, 1970; Koons, 1989); instead, they behave as density slabs, infinite in the y-direction. Dispersive properties of 2D and 3D ducted wave modes are principally different (Zudin et al., 2019), and therefore, the wave propagation properties in the 2D model may differ from in situ observations, especially where the azimuthal wave vector angle ϕ _k is concerned. Unfortunately, full-wave numerical investigations of 3D EMIC wave propagation are beyond our computational possibilities; furthermore, there is no suitable in-situ data we could use for comparison with the hypothetical 3D numerical results.

The dispersive properties of hydrogen and helium band EMIC waves strongly depend on the concentration of heavy ions (Lee et al. (2021) and references therein). The concentrations of the three ions in our simulations (η _p, η _He, and η _O) are set to be constant across the whole computational domain. Their impact can be seen in the behavior of characteristic frequencies below the proton gyrofrequency Ω_p. Specifically, the following three frequencies have an impact on the propagation and mode conversion of hydrogen band waves: the crossover frequency ω _cr, where the L-mode and R-mode become coupled and have linear polarization; the L-mode cutoff frequency ω _lc, which creates a stopband and leads to reflection of incident waves; and the ion hybrid frequency ω _ih, which can cause reflection of oblique R-mode waves. A similar triplet of frequencies exists in the helium band (but not in the oxygen band). Approximate formulas from Chen et al. (2014) are used in this paper to evaluate the characteristic frequencies.

The simulation code provides all six electromagnetic components and 3N _p cold current components with sampling rate f _s = 1/Δt on a fine spatial grid with cell size Δz ×Δx. In most practical cases, this data is too large to be stored. Instead, we save a small number of snapshots of the fine grid and a continuous time evolution on a coarse grid Δ_g z ×Δ_g x with sampling rate f _sg = 1/Δ_g t ≪ 1/Δt. Wave propagation properties are calculated in the time domain based on the coarse-grid data.

The data analysis process can be divided into four steps. First, the field components are transformed to the field-aligned system, (x, y, z) → (x′, y′, z′), where z′‖ B ₀ , x′ is perpendicular to the field line and lies in the x-z plane, and y′ = y. In the second step, the fields are converted to analytic signals with Hilbert transform. As the third step, we construct spectral matrices and average them over a short time interval (typically a small integer multiple of the average wave period T _avg given on the input). As the fourth and last step, we use the SVD methods (Santolík et al., 2003) and get normalized wave vectors, from which we can obtain various wave propagation properties: wave normal angle θ _k, azimuthal angle ϕ _k, B-ellipticity, and B-planarity. The angles are defined so that θ _k = 0° represents parallel propagation, θ _k = 180° anti-parallel propagation, and waves with ϕ _k = 0° propagate outward in the x′ direction. The Poynting vector and its polar and azimuthal angles θ _S and ϕ _S can be obtained directly from the cross-spectral components (Santolík et al., 2010). The instantaneous frequency of each B-field component is obtained by a simple forward difference of the analytic signal’s phase.

We also pick several grid cells (probes) and save the field data from those cells at a higher sampling rate. These high-resolution time series are used to construct spectrograms with the STFT (Short-Time Fourier Transform) method. Hann window with a 15/16 overlap is applied. Wave propagation properties in each time-frequency bin are obtained by the same SVD methods as described above, without any additional time averaging.

We use four different sets of input parameters to simulate the propagation of

Set 1: Constant-frequency hydrogen band wave,

A number of shared input parameters have no direct impact on wave propagation. Those are: Peak source amplitude B _w0 = 0.01, coarse grid sampling time T _avg/12, coarse grid spatial sampling λ _avg/4, probe data time step Δt/16, probe positions x′ = {−1,500 km, −750 km, 0 km, 750 km, 1,500 km} and λ = {0.5°, 2.5°, 5.0°, 10.0°, 15.0°, 20.0°} (5 ⋅ 6 = 30 probes in total), and 128 Boris time steps per electron gyroperiod (the smallest gyroperiod over the whole simulation box is taken).

All the other input parameters are listed in Table 1, with two exceptions. The density structure in the ducted case (Set 3) is composed of N _d = 3 Gaussians with characteristic equatorial widths 750 km ∼ σ _L = 0.118, amplitudes δn = 1 and centers L = L ₀ + {0, − 2σ _L, − 4σ _L}. Due to the properties of the Gaussian function, the sum of these three ducts creates a near constant elevation between L ₀ and L ₀ − 4σ _L with an increase of ∼ 1.2 δ n . Such structure can represent the plasmaspheric plume (Darrouzet et al., 2009). The rising-tone emission (Set 4) exhibits amplitude modulations of the source, with maxima of {0.5, 1.0, 1.0, 0.5}B _w0 at {0.125, 0.375, 0.625, 0.875}t _max and minima of {0.25, 0.50, 0.25}B _w0 at {0.25, 0.5, 0.75}t _max.

We first analyze the wave propagation properties of a simple hydrogen band EMIC wave packet with no amplitude modulations and a constant frequency. The corresponding input parameters can be found in Table 1, Set 1. In Figure 1, we plot the total magnetic and electric fields at time t = 25 s along with selected field components. Because the source is symmetric, we will discuss only the waves propagating towards z > 0 (northward); the fields propagating towards z < 0 (southward) quickly enter the damping region near box boundaries and dissipate.

Several important features of unducted EMIC wave propagation can be discerned from these snapshots. The magnetic and electric fields in Figures 1A,B display spatial oscillations with a period of λ _w/2, hinting at a rapid increase in ellipticity away from the source. Highly elliptical polarization suggests oblique wave vector direction, which can be confirmed by observing the angle between wave crests and field lines. When the wave frequency approaches the ion hybrid frequency (represented by the red curves in Figure 1), the magnetic field diminishes, while the electric field remains strong, confirming that the oblique EMIC wave is becoming electrostatic at the hybrid resonance.

The perpendicular component B _x′ plotted in Figure 1C disappears near the crossover frequency (magenta curves) because the wave becomes near linearly polarized with most of its magnetic amplitude in the y′ direction. Because the source width 2w _J is about 4.4 equatorial wavelengths, the plane wave approximation is not heavily violated, and E ⋅ B = 0 implies that the B _y′ component will have a spatial distribution similar to E _x′ (Figure 1D).

Figure 1E shows the magnetic field component B _z′, which reveals weak parts of the wave field moving across field lines, deviating from the expected quasi-parallel propagation of energy. Most of these can be shown to be right-handed and are related to mode conversions near characteristic frequencies—they will be described with the help of polarization analysis in the following paragraphs. However, some of these weak R-mode waves originate directly in the source. This observation may seem surprising since the source current density was obtained based on the assumption of a left-handed circular polarization. However, those calculations relied on the plane wave approximation in a homogeneous plasma, which is not exactly satisfied when the source is finite and the electron density in the radial direction changes by about 10% per wavelength. How the finite 1D source affects wave properties can be shown by constructing the magnetic field wave equation from Eqs 1, 2 (dropping the particle species index i) ∇ × ∇ × B = μ 0 ∇ × J s − 1 c 2 ∂ 2 B ∂ t 2 . (27) The circularly polarized current can be represented by J s = J s 0 x cos ⁡ ψ , sin ⁡ ψ , 0 , (28) where J _s0(x) represents the amplitude profile. If we take the curl of J _s , we get ∇ × J s = J s 0 k z ⁡ cos ⁡ ψ , J s 0 k z ⁡ sin ⁡ ψ , ∂ J s 0 / ∂ x sin ⁡ ψ . (29) The z-component will propagate into the calculation of B and cause deviations from the circular polarization. This behavior is confirmed in Figure 1E, where we show B _z′ to be nonzero in the source for z = 0, x ≠ 0. The resulting wave field can be decomposed into left-hand and right-hand polarized components, and thus the finite source supports both the L-mode and the R-mode. However, because the plane wave approximation is violated only weakly, the corresponding R-mode wave field is also weak. With decreasing w _J, the source is becoming increasingly point-like, supporting radiation in directions far away from θ _k = 0°. For the sake of completeness, we may also construct the electric field wave equation − ∇ × ∇ × E = μ 0 ∂ J s ∂ t + 1 c 2 ∂ 2 E ∂ t 2 . (30)

The spatial derivatives of current density are not present here; therefore, the electric field component E _z′ remains zero in the source, as documented in Figure 1F.

To better understand the wave propagation and polarization properties, we run SVD analysis on the coarse grid, following the methods described in Section 2.3. The results are shown in Figure 2, using the same snapshot as in Figure 1. Figure 2A presents the Poynting flux amplitude |μ ₀ S | in units of nT mV⁻¹. The lower threshold for all data is set to 10^−3.5 of the maximum Poynting flux in the chosen snapshot. As expected, the peak energy flux follows the starting field line L = 5.5, with only a slight deviation towards high L-shells, which is further confirmed by the low values of polar angle of the Poynting vector θ _S. The only regions where θ _S becomes large are the reflection region near the hybrid resonance, and the southern hemisphere (which will not be further discussed). Figure 2C shows that the wave starts as near left-hand circularly polarized, then becomes linearly polarized when crossing f = f _cr, and in the reflection region, a mixture of left-hand and right-hand polarization appears. The wave normal angle θ _k in Figure 2D shows a steady increase from near zero in the source up to 90° during reflection. However, the WNA values are somewhat noisy, especially near f _cr, near f _ih, and at the edges of the wave packet. Similar features are displayed by the azimuthal angle ϕ _k in Figure 2E, with the values jumping from 0° to 180°. At the two above-mentioned characteristic frequencies, the B-field polarization ellipse is degenerate (linear polarization), and thus the direction of the wave vector cannot be determined. At the wave packet edges, the variations come from the mixing of L-mode with the very weak R-mode. The planarity stays above 0.7 (Figure 2F), confirming that the use of the SVD analysis is meaningful.

A peculiar behavior can be seen near the source, where the azimuthal angle shows large eastward and westward deviations (Figure 2E). An explanation can be provided through the same calculations that led to Eq. 29. The additional part of the B _z component arising from the ∂/∂x gradient has a phase shift of 90° with respect to B _x, differing from the 180° shift expected in an oblique EMIC plane wave propagating in the meridional plane. This phase difference is demonstrated by Figures 1C,E where B _x′ is almost zero in the source while B _z′ attains its maximum or minimum at the same time. However, as long as the wave normal angle is small, the deviations in ϕ _k have little impact on wave propagation away from the source.

Apart from the main field-aligned packet, weaker fields with an oblique energy propagation direction appear in Figure 2, labeled W1 and W2. The W1 field has both θ _S and θ _k moderately or highly oblique ( > 50 ° ) , is right-hand elliptically polarized, and has origin in the weak R-mode radiated from the source. This origin can be confirmed by looking at B _z′ (Figure 1E) and tracing it back in time. A portion of the W1 field first propagates through the f _cr, experiences polarization reversal, reflects at the L-cutoff, and then goes through a second polarization reversal before arriving at the point of observation. The W2 field has similar properties but shows a smaller obliquity of θ _S and represents waves that were initially quasiparallel, did not experience reflection at the ion hybrid frequency and became unguided.

Neglecting the weak W1 and W2 fields, we may conclude that in a cold plasma with a high He + concentration (17%), an unducted hydrogen band wave that started with low values of WNA near the source will be entirely reflected back to the equator. To further confirm this conclusion, we show the evolution of the wave packet in three snapshots plotted in Figure 3. Figures 3A–F cover the initial stage before reflection, where we can clearly see the W1 field escaping away from L = 5.5 and the main field going through polarization reversal. We skip over the intermediate stage already discussed in Figure 2 and go to the final stage at t = 45 s, where the majority of the wave packet has been reflected. The reflected wave goes through a second polarization reversal and reaches the source region as a very oblique, left-hand highly elliptically polarized EMIC emission. Due to the slight deviation of S from the field-aligned direction, the reflected wave passes above the center of the wave source when it returns to the equatorial plane. A space probe flying through the source region can, therefore, easily miss the reflected wave, depending on the probe’s precise position and velocity. Our results demonstrate that more data from multi-spacecraft observations at close separations are needed to evaluate the occurrence and physical properties of reflected EMIC waves.

In a cold plasma with the three ions H+, He+, and O+, the dispersive properties of the helium band EMIC wave are generally similar to those of the hydrogen band wave, with the most significant difference being in wavenumber. However, because He-band is associated with higher densities (Meredith et al., 2014), even the wavenumbers can be close in value. One clear distinction appears once waves pass through the L-mode stopband: while the H-band waves will encounter another stopband at higher latitudes, the He-band can propagate unimpeded down to ionospheric altitudes. Moreover, low oxygen concentrations in the plasmatrough can push the characteristic frequencies close to Ω_O.

In the simulation run with input parameter Set 2, we choose ω ₀/Ω_He0 = 0.6, which is the same value as ω ₀/Ω_p0 in Set 1. The plasma frequency is doubled, which still makes the wavelength larger than in the hydrogen band case. Furthermore, the source extent is decreased from ±1,500 km to ±1,075 km. The full source width 2w _J is now only about two equatorial wavelengths, which leads to enhanced radiation in oblique directions. This change in the radiation pattern is demonstrated in Figures 4A,B, where the magnetic field components B _x′ and B _z′ show oblique wave crests emanating from the source. The power of these weak oblique fields is about two orders of magnitude below the peak of the main packet (Figure 4C), making them more significant than in the wide-source case from Figure 2. The polar angle θ _S reveals that some of these weaker waves stay at latitudes below 10°, displaying very oblique propagation of energy, while some are more parallel and quickly (i.e., at time t = 27.5 s) reach latitudes up to above 20°. The first group is marked W1 in Figure 4E and has two components. The outward propagating component is right-hand polarized and quickly reaches the damping boundary—its properties are similar to the W1 field from Figure 2, hence the shared label. The second component propagates inward, experiences polarization reversal, and reflects at the L-cutoff. The second group, labeled W2, propagates in a quasiparallel fashion down to the lower box boundary without experiencing any notable changes in propagation properties, becoming unguided after passing below f _ih. We may notice in Figure 4F that at the boundaries between L-mode-dominated and R-mode-dominated regions, the wave normal angle appears to be near 90°.

The narrow source not only affects the amplitude of the weak R-mode, but also changes the propagation properties of the main L-mode packet. This is easier to observe during later times of the wave field evolution, as demonstrated by the three snapshots in Figure 5. At t = 17.5 s (Figures 5A–C), the separation into the main packet and secondary R-mode packets is already clear. At a later time, t = 32.4 s (Figures 5D–F), we notice that a part of the main packet near the bottom edge, labeled W3 in Figure 5E, has a near-circular left-hand polarization and quasiparallel wave vector. In Figures 5G–I (snapshot t = 40.0 s), this W3 field is shown to propagate through the crossover frequency without any significant loss of wave power and remain mostly left-hand polarized, while the more oblique portion of the main wave packet undergoes polarization reversal. However, a small amount of wave energy is transferred into a reflected R-mode component near f = f _cr, labeled W4 in Figure 5H. The above-discussed weak fields W1 and W2 are also marked for comparison with Figure 4.

After inspecting the propagation of unducted EMIC waves in Sections 3.1 and Sections 3.2, we turn to ducted propagation on steep density gradients (input parameter Set 3). In this ducted H-band simulation, the source width is smaller than in Set 1, w _J = 1,125 km, but the wavelengths are shorter due to the increased density in the ducting structure. Because of the large radial density gradients, the wavelength in the source cannot be well represented by a single value. Nonetheless, in the chosen setup, the source is still wide enough to emit only negligible power into the oblique directions.

The B _x′ component plotted in Figure 6A shows stark differences from the unducted picture in Figure 1C. When reaching f _cr, the wave crests follow the field line, and the L-value of this field line matches well with the region of steep density drop off (see Figure 6B for a plot of 2D ω _pe/Ω_e0 distribution). As shown in Figures 6C,D, most of the wave power reaches the hybrid resonance and passes towards higher latitudes. However, a significant amount of wave power becomes reflected near the f = f _cr and f = f _lc surfaces. In Figure 6E, we can see that the quasiparallel ducted EMIC wave turns into a mixture of L-mode and R-mode after passing through the crossover. The left-handed part reflects to higher L-shells and goes through a polarization reversal—we label it S3, as it has a similar propagation path to W3, except for being stronger due to ducting. S2’ is an R-mode wave, which has similar properties to the weak field W2 but originates in the polarization reversal instead of being generated directly in the source. Additionally, a weak part of the main packet gets reflected at the crossover and becomes right-hand polarized (W4 in Figure 6E). Due to the mixture of R-mode and L-mode waves with similar power, the wave normal angles plotted in Figure 6F are difficult to interpret, but the quasiparallel component at the center of the wave packet can still be traced.

Values of θ _k below the f = f _ih surface suggest moderate obliquity, which means that these waves are decoupled from the R-mode at ω < Ω _O , and most of their energy will be reflected before reaching low altitudes. However, the simulation would give nearly identical results for He-band with a similar wavenumber (requiring a high-density background), up to minor differences related to the ion composition. The S2’ wave field would then propagate through down to the ionosphere.

In the previous sections, the current density source was not dynamic, except for the slow changes in amplitude that formed the single-peaked wave packet envelope. We now modulate the wave field into four subpackets and introduce a constant frequency chirp from the initial frequency ω ₀ = 0.45Ω_p0 to the final frequency ω ₁ = 0.65Ω_p0. As listed in Table 1, the concentration of He+ and O+ has been decreased to prevent the f ₀ = f _cr surface from crossing the source. The source width is very slightly (by ∼ 7 % ) smaller than in Set 1 and the initial wavelength is larger, making the source relatively more narrow. This can be noticed in the B _y′ plot in Figure 7A, where a weak quasiparallel field appears at the inner edge of the wave field. In Figure 7B, we show the instantaneous frequency at time t = 15 s, shortly before the reflection of the first subpacket at the hybrid resonance. At the local minimum between the first and second subpacket, which is located near 10° of latitude, the frequency can be seen to have a localized decrease, despite the linearly rising frequency of the source. These are the effects of group velocity dispersion consistent with ∂V _g/∂ω < 0, which have been previously discussed in the context of rising-tone chorus elements by Hanzelka and Santolík (2022).

To provide another view on the evolution of wave frequency and the effects of reflected waves, we process data from two selected simulation probes, P1 and P2, whose position is shown in Figure 7B (P1 near the source, P2 at the f _cr crossing of the first subpacket). Unlike in the presentation of the constant-frequency wave propagation (Section 3.1), we can show the whole time evolution in one plot, but we are limited to a single point in space. Figure 7C presents the amplitude envelopes of all three magnetic field components B _x′, B _y′, and B _z′, as well as the total magnetic field | B |. Before t ∼ 30 s, the four northward propagating subpackets have a negligible parallel component, and the peak amplitude is near equal to B ₀/100, as dictated by the source properties. The two perpendicular components have similar magnitudes, suggesting circular polarization, which is further confirmed by the negligible oscillations in | B |. After t ∼ 30 s, the reflected wave packet passes over probe P1, but its amplitude is diminished, with the second subpacket having less than 15% of its original amplitude. The B _y′ component dominates, and | B | exhibits strong oscillations at two times the wave frequency, which are signs of linear polarization and high obliquity in the hydrogen band L-mode. The frequencies near the source (Figure 7D) are linearly growing, except for minor oscillations in the frequency derived from B _z′. The chirp rate within subpackets of the reflected wave is higher than the initial value, which is again the effect of group velocity dispersion that we already noted when describing Figure 7B.

Moving to probe P2, we observe the behavior of two linearly polarized waves propagating in opposite directions. This overlap creates multiple very short subpackets associated with large variations in the instantaneous frequency. When such a structure is observed, we must separate the two modes, either by the Hilbert-Huang transform (HHT; for an application on EMIC waves, see (Ojha et al., 2021)) or by inspection of time-frequency spectrograms. We choose the latter method and plot spectrograms of power spectral density (PSD) and propagation properties in Figure 8. The PSD near equator, as shown in Figure 8A, confirms that the reflected emission has a considerably lower power (down by almost two orders of magnitude) than the forward propagating wave. Depending on the signal-to-noise ratio of the original wave packet, only parts of the reflected riser may be visible, or none at all. It is of note that the rising-tone element has a considerable spectral width, which is the consequence of the Fourier uncertainty principle for short wave packets. Even with the 93.75% overlap of the STFT time windows, the frequency resolution is too low to determine whether the individual subpackets are chirping—the line plots obtained from Hilbert transform are better suited for this purpose but might require some form of mode decomposition like the one included within HHT. Figure 8 supports our previous assessment of the reflected wave’s propagation properties, clearly showing the near-linear polarization and highly oblique wave normal angle.

Finally, Figures 8E–H show spectrograms constructed from the probe P2 data. Due to the close frequencies of the two wave packets, it is not immediately clear from the power spectrum in Figure 8E that we are indeed observing two risers propagating in opposite directions. Fortunately, the polar angle θ _S in Figure 8H shows a clear division into two northward and southward propagating elements. Figures 8F,G further show that one element has θ _k about 65° and ellipticty of −0.2, while the other has θ _k ≈ 100° and positive ellipticity < 0.1 (these values slightly differ between subpackets). However, in practical applications with the presence of noise, these small differences in WNA and ellipticity could not serve as reliable criteria, and the Poynting flux data would be needed. Before we conclude this section, we must note that the probe P2 was shifted to a slightly higher L-shell than P1 to ensure similar wave powers of the forward-propagating and reflected element—a property that is not necessarily satisfied within real spacecraft data.