Normal mode analysis is a technique that finds the resonant frequencies of any oscillating system. Before using this technique to describe waves in a plasma, it is useful to consider a simple example of two masses with springs on each side, as shown in Figure 2. The equations of motion for the position of each mass are found by applying Newton’s second law, m d 2 d t 2 x 1 = − 2 κ x 1 + κ x 2 , (1) m d 2 d t 2 x 2 = κ x 1 − 2 κ x 2 , (2) where κ is the spring constant from Hooke’s law.

This set of differential equations can be written in matrix form as m d d t 2 x 1 x 2 = − 2 κ κ κ − 2 κ · x 1 x 2 . (3)

To find the normal modes of this system, Equation 3 is Fourier transformed, so effectively d d t → i ω . The system of equations then becomes − m ω 2 x 1 x 2 = − 2 κ κ κ − 2 κ · x 1 x 2 . (4)

In this form, the goal of normal mode analysis becomes apparent: create and solve an eigenvalue equation. The left-hand side shows the eigenvalues m ω 2 , while the right-hand side is the linear transform corresponding to the forces on each mass. To complete this example, Equation 4 is rearranged: − m ω 2 + 2 κ − κ − κ − m ω 2 + 2 κ · x 1 x 2 = 0 . (5)

The solution of Equation 5 is only possible if either x 1 = x 2 = 0 (the trivial solution) or if the matrix is not invertible. Setting the determinant to 0 and solving the resulting polynomial gives the eigenvalues of ω : ω 2 = κ m 2 ± 1 . (6)

The two positive roots of ω 2 in Equation 6 are the resonant frequencies corresponding to the normal modes of the system. Since this system is linear, the general solution to the equations of motion (Equations 1, 2) is a superposition of these two normal modes. This technique can be generalized to any dimension system of linear differential equations, including the set of equations that describe the motions of a plasma.

The goal of normal mode analysis in a plasma is the same as the simple example above: to create an eigenvalue equation from a system of differential equations, then solve for the resonant frequencies of the plasma. The motions of a plasma are described at the kinetic level by the Vlasov equation for each species s: ∂ F s t , x → , v → ∂ t + v → · ∇ x F s t , x → , v → + q s m s E → t , x → + v → × B → · ∇ v F s t , x → , v → = 0 . (7)

The Vlasov equation is effectively a total time derivative of the distribution function, with the Lorentz force as the acceleration. To close this set of equations, the electric and magnetic fields need to be specified. While the following solutions and ideas work for general electromagnetic waves, we will restrict the analysis in this paper to electrostatic solutions since those are the wave modes measured by Thomson scatter radars. Therefore, Gauss’ law is appropriate to close the system. Additionally, we do not include a collision operator on the right-hand side of Equation 7 as collisions act to damp waves, but they do not significantly affect the resonant frequencies.

As defined in Equation 7, the Vlasov equation is nonlinear since the Lorentz force multiplies the velocity derivative of the distribution. For normal mode analysis to apply, the Vlasov equation and Gauss’ law need to be linearized and then Fourier–Laplace transformed. This is done using standard techniques, with each variable being decomposed as a zero-th order term and a first-order perturbation, such as F s = F 0 s + F 1 s , with any resulting second-order terms discarded. The full linearization process, including the justification for dropping second order terms, is detailed in Longley (2024). The linearized Vlasov equation is then − i ω − k → · v → − i γ F 1 s ω , k → , v → + q s m s E → 1 ω , k → · ∇ v F 0 s v → + q s m s v → × B → · ∇ v F 1 s ω , k → , v → = F 1 s t 0 , k → , v → (8) where k is the wavenumber from the spatial Fourier transform and ω is the frequency from the Laplace transform in time. Note that taking the Laplace transform in time creates the initial value term ( t 0 ) on the right-hand side. Now linearizing Gauss’ law, i k → · E → 1 ω , k → = 1 ϵ 0 ρ 1 ω , k → , (9) where ρ 1 = ∑ s q s n 1 s is the charge density. Since the number density is defined as n 1 s = ∫ d v → 3 F 1 s , this creates a closed system of equations for E 1 and n 1 s for each species.

The Thomson scattering spectra is defined as S ω = n 1 e ω 2 , so typically Equations 8, 9 are solved for the perturbed electron density n 1 e (Froula et al., 2011). However, to construct an eigenvalue equation for normal mode analysis, this system of equations is instead solved for the perturbed electric field. This is done by first integrating Equation 8 to obtain the perturbed densities, then substituting those into the charge density in Gauss’ law (Equation 9). This is easiest to demonstrate in the unmagnetized limit, where B = 0 . Integrating Equation 8 over all velocity space obtains n 1 s ω , k → = − i q s m s ∫ d v → 3 E → 1 ω , k → · ∇ v F 0 s v → ω − k → · v → − i γ + i ∫ d v → 3 F 1 s t 0 , k → , v → ω − k → · v → − i γ . (10)

Taking Equation 10 for the densities of electrons and a single ion species, the electric field in Equation 9 becomes i k → · E → 1 ω , k → = 1 ϵ 0 − i e 2 m i ∫ d v → 3 E → 1 ω , k → · ∇ v F 0 i v → ω − k → · v → − i γ + i e ∫ d v → 3 F 1 i t 0 , k → , v → ω − k → · v → − i γ − i e 2 m e ∫ d v → 3 E → 1 ω , k → · ∇ v F 0 e v → ω − k → · v → − i γ − i e ∫ d v → 3 F 1 e t 0 , k → , v → ω − k → · v → − i γ . (11)

the meaning of Equation 11 is clearer if the terms with E 1 are grouped together. Also recognizing that for electrostatic waves, k → and E → 1 are colinear ( k → · E → 1 = k E 1 ), we obtain i k E 1 ω , k → = − i e 2 k ϵ 0 E 1 ω , k →   1 m i ∫ d v → 3 k → · ∇ v F 0 i v → ω − k → · v → − i γ + 1 m e ∫ d v → 3 k → · ∇ v F 0 e v → ω − k → · v → − i γ + i e ϵ 0 ∫ d v → 3 F 1 i t 0 , k → , v → ω − k → · v → − i γ − ∫ d v → 3 F 1 e t 0 , k → , v → ω − k → · v → − i γ . (12)

The initial value terms ( F 1 s t 0 , k → , v → ) are not multiplied by E 1 , so Equation 12 is not a linear transformation of E 1 as it does not satisfy the additive property. These initial value terms are analogous to a source term for the driven oscillator and are the primary subject of Section 3. However, for normal mode analysis, the steady-state behavior of the system is of interest, so initial perturbations can be ignored. Dropping these initial value terms, we arrive at the desired eigenvalue equation: i k E 1 ω , k → = − i e 2 k ϵ 0 E 1 ω , k → 1 m i ∫ d v → 3 k → · ∇ v F 0 i v → ω − k → · v → − i γ + 1 m e ∫ d v → 3 k → · ∇ v F 0 e v → ω − k → · v → − i γ . (13)

Physically, this shows the left-hand side as eigenvalues for the density perturbations on the right-hand side.

The integrals in Equation 13 are defined as the susceptibility of each species: χ s ω , k → = ω p s 2 k 2 ∫ d v → 3   1 ω − k → · v → − i γ k → · ∂ f o s ∂ v → (14) where ω p s 2 = n 0 s e 2 / m s ϵ 0 is the plasma frequency for each species, and the distribution is now normalized using the notation F 0 s = n 0 s f 0 s , so that ∫ d v → 3 f 0 s = 1 . We can then simplify Equation 13 using Equation 14 E 1 ω , k → = − E 1 ω , k → · χ i ω , k → + χ e ω , k → . (15)

Solving for the perturbed electric field in Equation 15, 1 + χ i ω , k → + χ e ω , k → · E 1 ω , k → = 0 . (16)

The dielectric function of a plasma is defined as ϵ = 1 + χ i + χ e , so the final result is ϵ ω , k → · E 1 ω , k → = 0 . (17)

Equation 17 is the basis for the normal mode analysis of a plasma, and considerable literature exists on deriving equivalent forms (e.g., Bekefi, 1966). Solving Equation 17 means either E 1 = 0 (the trivial solution) or ϵ = 0 is required. Therefore, solving for roots of the dielectric function will obtain the resonant frequencies of a plasma. Note that in Equation 17, the dielectric function is a scalar because the electrostatic approximation was used. A general electromagnetic solution will lead to a 3 × 3 matrix for the dielectric function, and the roots are obtained by setting the determinant of the matrix equal to zero.

By writing Equations 16, 17 in terms of susceptibilities, the assumptions of electrostatic waves ( B 1 = 0 ), unmagnetized plasma ( B 0 = 0 ), and no collisions can be relaxed by using the appropriate susceptibilities. Standard plasma physics texts (Bellan, 2006; Froula et al., 2011) have derived these more general susceptibilities. Here, we make use of collisionless, electrostatic susceptibilities with Maxwellian velocity distributions. These are listed in Supplementary Appendix A, for both the unmagnetized and magnetized cases. Furthermore, Supplementary Appendix B shows how the Thomson scatter spectra are calculated from these susceptibilities.

The simplest solution for a root of the dielectric function is the Langmuir mode. This wave mode is anticipated to occur at a frequency near the plasma frequency, and therefore ion dynamics can be neglected (Longley et al., 2021). Plasma lines are well known to be enhanced by photoelectrons (Longley et al., 2021), but for simplicity we will only analyze the case of thermally driven plasma lines such as those detected in Vierinen et al. (2017). Using the unmagnetized electron susceptibility, the dielectric function at high frequencies is (Supplementary Appendix A) ϵ = 1 + χ e = 1 + α 2 1 − 2 x e Daw x e − i π x e e − x e 2 , (18) where the parameter α is the ratio of the wavelength to the Debye length, defined as α = 1 k λ D (19)

Figure 3 plots the real and imaginary parts of Equation 18 for different sets of plasma parameters. The resonant frequency is obtained by solving for the roots of ϵ = 0 , which correspond to eigenvalues of ϵ · E 1 = 0 . The parameters chosen in Figure 3 illustrate three different cases of the dielectric function for the Langmuir mode. For the T e = 1000   K curve, the parameter α = 3.5 means that the wavelength is significantly larger than the Debye length. With large α and small Landau damping due to the low temperature, there exists a single root to ϵ = 0 that is the Langmuir frequency. Raising the temperature in Figure 3 to 2500 K lowers α to 2.2 and creates appreciable Landau damping, so that Im ϵ ≠ 0 for frequencies near the root of Re ϵ = 0 . This makes a solution to ϵ = 0 impossible, but physically it represents a damped wave. A root to the equation Re ϵ = 0 when damping is present is called a “quasinormal mode” (Wikipedia, 2023). Note that Landau damping is mathematically described by the i π x s e − x s 2 term in Equation 18, and therefore the damping is never exactly 0 for finite x s . Physically, Landau damping is where particles with velocities near the phase velocity of the wave ( v = ω / k ) efficiently take energy from the wave, leading to a decrease in the wave’s amplitude (Chen, 2016). The amount of damping is dependent on the temperature of the distribution, which describes how many particles have velocities near the phase velocity of the wave. Even at the lower temperature of T e = 1000   K in Figure 3, some Landau damping is present, and this case could be strictly defined as a quasinormal mode. However, a practical distinction is applied where cases with Im ϵ ≈ 0 are considered normal modes.

Increasing the temperature once more to 5000 K in Figure 3 leads to α = 1.6 , and the wavelength is now the same order of magnitude as the Debye length. Physically, the Debye length is an exponential scale length describing distances where electrons can easily reconfigure to shield any charge imbalances. Waves propagate through a plasma by creating and sustaining charge imbalances. It is therefore no surprise that for the α = 1.6 case, there is no solution to Re ϵ = 0 . However, a distinct spectral peak in the scattering power still appears for this case.

The transition of the plasma line/Langmuir mode from a normal to a quasinormal mode, and to what will later be defined as a non-resonant mode, appears to depend on the α parameter. We can further investigate this by directly solving for roots of the real part of the dielectric function, 0 = Re ϵ = 1 + α 2 1 − 2 x e Daw x e . (20)

The standard approximation is to anticipate ω ≈ ω p e , which in most conditions means x e ≫ 1 . The Dawson function can then be Taylor-series expanded for large x e , so that 2 x e Daw x e = 1 + 1 2 x e 2 + 3 4 x e 4 + O x e − 6 . Putting this into Equation 20, 0 = 1 − α 2 1 2 x e 2 + 3 4 x e 4 . (21)

This creates the biquadratic equation in x e which is solved with the quadratic formula. With the assumption that α 2 ≫ 1 and substituting x e = ω / k v t h , e , Equation 21 simplifies to the well-known dispersion relation for Languir waves (Equation 22): ω 2 = ω p e 2 + 3 2 k 2 v t h , e 2 (22)

For the Langmuir mode, it is easy to solve for a root of the dielectric function if the approximations of x e ≫ 1 and α 2 ≫ 1 can be made. Without these assumptions, an analytical solution is not as simple, and in some cases not even possible. An easier and more interesting question is to make no assumptions on x e and ask what parameters are needed for a root of Re ϵ = 0 to exist. Looking at the real part of the dielectric function in Equation 20, we see that at x e = 0 , the dielectric is ϵ = 1 + α 2 , which is strictly positive. Therefore, a root is only possible if the right-hand side becomes negative for some value of x e . This gives the condition of 0 > 1 + α 2 1 − 2 x e Daw x e (23)

The only free plasma parameter in Equation 23 is α , so we solve for the values of α necessary for a root to exist. The inequality is easiest to satisfy when the term 1 − 2 x e Daw x e is at its minimum value. This occurs for x e ≈ 1.50198 . Substituting this for x e and rearranging, we obtain the condition for a root to exist: α > 1.874 (24)

For typical Thomson scatter experiments in the ionosphere, α ≫ 1 , and therefore a root to the dielectric function is expected near the plasma frequency. Note that the magnetized form of the dielectric function can slightly modify the condition in Equation 24. Figure 3 shows that for low densities and high temperatures, it is possible for the above condition to not be satisfied and therefore no root will exist. Observations of the plasma line at altitudes above 1500 km are reported in Hagen and Behnke (1976), with observed spectra at α ≈ 1.2 looking similar to that in Figure 3 with α = 1.6 .

The unmagnetized Langmuir mode and its relation to the plasma line are the easiest solutions to the dielectric function possible. Yet the solution still encounters several problems. The first is that if damping is present, we cannot solve for ϵ = 0 and instead need to solve Re ϵ = 0 . Figure 3 shows there are actually two roots to this equation, and we intuitively choose the root with less damping. Nonetheless, we have not justified the exclusion of the other roots nor defined criteria to assure a root-solving algorithm finds the correct root. Furthermore, while it is expected that the Langmuir mode usually exists in the ionosphere, the ion-acoustic and electrostatic whistler modes are often cutoff with no solutions to Re ϵ = 0 (see Section 4). However, Thomson scatter experiments still measure strong scatter in ion and gyro lines, necessitating a more robust characterization of what a plasma wave mode is.

In deriving the dielectric function in Section 2.2, the initial value terms from the Laplace transform in time were dropped so that a linear transformation of the electric field could be written and solved for eigenvalue frequencies. Keeping the initial value terms, the equation for the electric field becomes ϵ ω , k → E 1 ω , k → = i e k ϵ 0 ∫ d v → 3 F 1 i t 0 , k → , v → ω − k → · v → − i γ − ∫ d v → 3 F 1 e t 0 , k → , v → ω − k → · v → − i γ (25)

Without the initial value terms on the right-hand side, this is Equation 17.

Equation 25 is no longer solvable for eigenvalues of the dielectric function. Mathematically, this is because the initial value term on the right-hand side means this is no longer a linear transformation of the electric field (failing the additive property f x + y = f x + f y ). Physically, Equation 25 thus no longer describes the normal modes of the system that will naturally exist with small initial perturbations. Instead, the initial positions of particles will create an electric field with the strength being E 1 ω , k → = 1 ϵ ω , k → · i e k ϵ 0 ∫ d v → 3 F 1 i t 0 , k → , v → ω − k → · v → − i γ − ∫ d v → 3 F 1 e t 0 , k → , v → ω − k → · v → − i γ (26)

For normal modes of the system, ϵ → 0 at a resonant frequency and amplitude E 1 of the wave can be large for a given initial perturbation (RHS of Equation 26). Going one step further, at non-resonant frequencies, the amplitude of oscillation will be highest when the dielectric function is smallest. It is this behavior of the dielectric that is of interest.

Taking a step back, the initial value terms of every particle are unknowable for a plasma experiment. Instead, an ensemble average is taken of Equation 25 to find the average electric field strength, weighted by the likelihood the plasma started in a particular initial state. The ensemble average is defined in Froula et al. (2011) as X = ∫ d v   X v   P v ∫ d v   P v . (27) where the zero-th order distribution f 0 v is the probability P v of finding particles of species s at a given initial position. Formally, Equation 25 is squared and ensemble averaged using Equation 27 to obtain the average electric field (without squaring, the result is E 1 = 0 since there is an implied sinusoidal dependence when linearizing). Furthermore, this ensemble average is applied over all space and time, taking the limit as the volume and time going to infinity (note that γ = 1 / T ). Then ϵ ω , k → 2 E 1 ω , k → 2 = lim T , V → ∞ 1 T V e 2 k 2 ϵ 0 2 ∫ d v → 3 F 1 i t 0 , k → , v → ω − k → · v → − i γ − ∫ d v → 3 F 1 e t 0 , k → , v → ω − k → · v → − i γ 2 (28)

In squaring the left-hand side of Equation 28, there will be cross terms, but the standard treatment is to drop these by assuming that the initial positions of electrons and ions are uncorrelated (Froula et al., 2011). This assumption can be relaxed, though the resulting cross terms will only lead to an initial transient that decays as 1 / t (Froula et al., 2011). Carrying out the ensemble average on the right-hand side, ϵ ω , k → 2 E 1 ω , k → 2 = S e ω , k → + S i ω , k → , (29) where we define the source terms as S s ω , k → = lim T , V → ∞ 1 T V e 2 k 2 ϵ 0 2 ∫ d v → 3 F 1 s t 0 , k → , v → ω − k → · v → − i γ 2 . (30)

As an example, the Maxwellian distribution can be used for the initial distribution function F 1 s t 0 , and in the unmagnetized limit, the source term in Equation 30 evaluates to S s ω , k → = n s e 2 k 2 ϵ 0 2 π k v t h , s exp − ω 2 k 2 v t h , s 2 (31)

Comparing Equation 31 to Appendix Equation B3 shows the connection between the wave source terms S and the modified distributions M that describe Thomson scatter (Supplementary Appendix B). In general, these terms are proportional through the relation shown in Equation 32 S s ω , k → = n s e 2 k 2 ϵ 0 2 M s ω , k → (32)

Chapter 9 of Nicholson (1983) derives a similar expression to Equation 28 by neglecting ion dynamics and considering the electric potential of numerous moving test charges. Nicholson (1983) calls this result “fluctuations in equilibrium” but only applies the analysis to resonant Langmuir waves (i.e., ϵ = 0 condition). Furthermore, Bekefi (1966) develops the formalism for “non-resonant” waves driven by motions of charges but does not provide a treatment of the ensemble averaged system.

Equation 29 is the desired result for interpreting the existence of density fluctuations in a plasma. The normal modes can still be obtained by setting the source terms S equal to 0 and solving for roots of ϵ , but the strongest oscillations in a plasma are not necessarily at the resonant frequencies. The steady-state behavior of the plasma is obtained by taking an ensemble average, and therefore the source terms show how waves are continuously generated across all frequencies.

For physical intuition, we return to the mass-on-spring analogy of Section 2.1. The eigenvalue/eigenvector relation in Equation 5 can be written to include a sinusoidal driving force on each mass, giving the new relation − m ω 2 + 2 κ − κ − κ − m ω 2 + 2 κ · x 1 x 2 = A 1 ⁡ sin ω d t A 2 ⁡ sin ω d t (33)

In this example, the frequency of the driving force ω d will dictate the frequency each mass oscillates the same as the original eigenvalues of the matrix on the left-hand side of Equation 33. It is this balance between an external driving force (source term) and the system’s internal response (normal modes) that determines the full oscillating spectrum of the system.

The source term in its simplest form (Equation 31) is the velocity distribution evaluated at the condition v = ω / k (Equation 29). This is the Cherenkov condition for wave generation by particle motion, with more particles at a given velocity leading to stronger waves (Nicholson, 1983). This creates an analogy with the driven harmonic oscillator, where the source terms S e and S i act as a continuous driving force for waves at the frequency ω = k → · v → . In this interpretation, E 1 2 is the average amplitude of the oscillations, dictated by the value of the response function ϵ 2 .

Since ϵ = 1 + χ e + χ i , we can interpret the resonant ( ϵ = 0 versus non-resonant ( ϵ ≠ 0 ) response to driven oscillations in terms of the susceptibilities mean. If χ s is large at a given frequency, then species s is able to efficiently reconfigure and cancel out an applied electric field. An intuitive example of this is that for low frequency waves ω ≈ 0 , the susceptibility will be very high (e.g., Figure 3) since the long period of the wave allows plenty of time for electrons and ions to reconfigure and cancel out the wave’s electric field. The real part of χ s will act in phase with the wave, whereas the imaginary part of χ s acts out of phase with the wave and therefore will damp it out. If Re χ s ≈ 0 , then the plasma is not reconfiguring to cancel out the applied electric field but is instead moving with the electric field in a way that continues the wave’s propagation. Values of Re χ s that are small, but non-zero can be interpreted as the plasma trying to propagate the wave but being unable to fully do so in each cycle.

With this interpretation of the dielectric function, we are able to explain the concept of a non-resonant wave mode (e.g., Bekefi, 1966). Given a source term S ω that is constant at all frequencies, the amplitude of oscillation E 1 ω 2 is largest when the plasma is best able to propagate the wave. Therefore, when the dielectric function has no roots ( ϵ ω = 0 ), a local minimum of the dielectric function ϵ ω 2 will represent non-resonant waves where the source term is continuously driving the wave and the plasma is able to partially continue the oscillation but at a lower amplitude than if a normal mode resonance existed. Without the continuous driving of S ω , such waves would quickly decay and be unobservable. This interpretation also better characterizes finding roots of the real part of the dielectric when there is still an imaginary damping part— Re ϵ = 0 and Im ϵ ≠ 0 . Therefore, for heavily damped roots ( Re ϵ = 0 and Im ϵ ≫ 0 ), no wave will exist.

For a source term that varies with frequency, the largest oscillations will occur at a balance between the minima of the dielectric function and the maxima of the driving source. This balance can lead to the strongest scatter occurring at frequencies shifted away from the resonant or non-resonant frequencies, as the next sections will demonstrate.

For the ion-acoustic mode, both the electron and ion susceptibilities are important, and therefore the dielectric function is ϵ = 1 + χ e + χ i . It is expected from fluid theory that the ion-acoustic frequency will be approximately ω ≈ k v t h , i , which means x i ∼ 1 (Appendix Equation A2), and the Dawson function in the ion susceptibility cannot be Taylor expanded in either the large or small limits. The electron susceptibility, however, can be simplified. Since the Dawson functions are evaluated at normalized frequencies (Appendix Equation A2), the electron and ion arguments are related by x i x e = m i m e T e T i . For the ion-acoustic mode, x i ≈ 1 , so x e ≪ 1 in this frequency range so long as T e and T i are approximately within an order of magnitude. The electron susceptibility can then be Taylor-expanded to first order in x e as χ e ≈ α 2 1 − i π x e (34)

With Equation 34, the dielectric function at low frequencies is then ϵ ω = 1 + α 2 1 − i x i π m e T i m i T e + T e T i α 2 1 − 2 x i Daw x i − i π x i e − x i 2 (35)

The presence of the Dawson function requires a numerical solution to find any roots of Equation 35.

In the top panel of Figure 4, the real and imaginary parts of Equation 35 are plotted, showing that roots for Re ϵ = 0 only occur for specific ratios of the electron and ion temperatures. Exact criteria for when a root occurs can be derived from Equation 35. Firstly, note that at x i = 0 , the real part of the dielectric function is strictly positive— R e ϵ = 1 + α 2 1 + T e T i > 0 . Therefore, a root only exists if there is a value of x i where the real part of the dielectric is negative. This condition is 0 > 1 + α 2 + T e T i α 2 1 − 2 x i Daw x i (36)

Since 1 + α 2 > 0 , the last term in Equation 36 needs to be negative. This is easiest to satisfy if 1 − 2 x i Daw x i is at its minimum value, which happens at x i ≈ 1.50198 . Then T e T i > 1 + α − 2 min 1 − 2 x i Daw x i (37) T e T i > 1 + α − 2 0.28475 ≈ 3.51 1 + 1 α 2 (38)

Typically, α 2 ≫ 1 for Thomson scatter experiments in the ionosphere, so the ion acoustic mode only has a root for the dielectric function when T e T i ≳ 3.5 . In the ionosphere, the temperature ratio is rarely greater than ∼3 (Aponte et al., 2001), and therefore no root to the real part of the dielectric exists according to Equations 37, 38. Nonetheless, the ion line is always observed in the collective scatter regime where α ≫ 1 . Moreover, the ability to detect the ion line is arguably the single defining feature of the incoherent scatter radar (ISR) class of Thomson scatter radars. This highlights the problem of associating the resonant solution of the ion-acoustic mode with the ion line in Thomson scatter spectra.

The existence of an ion line in Thomson scatter experiments can be explained by the analogy with a driven oscillator described in Section 3.2. Waves will continuously be generated at low frequencies through Cherenkov radiation by particles moving at v = ω / k , with the source terms providing the strength of wave generation. The response of the plasma to these generated waves is quantified by the dielectric function, with local minima of ϵ ω 2 being defined as non-resonant wave frequencies. To test this idea, the middle panel of Figure 4 shows ϵ ω 2 , and the bottom panel shows the corresponding ion line spectra. Because the ion-acoustic mode is heavily Landau-damped (imaginary part of ϵ ), the minima values of the dielectric function do not correspond exactly to the roots of Re ϵ = 0 when they exist. Similarly, the peaks in the ion line spectra do not correspond to either the resonant or non-resonant frequencies of the ion-acoustic mode since the driving source term is also important.

The ion line is effectively unmagnetized for most aspect angles (Milla and Kudeki, 2011). Therefore, the driving source term for ion-acoustic waves is the Maxwellian distribution given by Equation 31. For electrons, the argument of the Maxwellian is ω k v t h e ≪ 1 , so the electrons drive waves of equal strength at all frequencies relevant to the ion line. However, the ion source term will drop off significantly on the range of frequencies relevant to the ion line. This means that while the plasma responds best at frequencies around ω k v t h i ≈ 2 (see Figure 4), the strongest driving force is at lower frequencies, ω k v t h i < 1 . The balance between the response of the plasma (dielectric) and the continuous generation and driving of waves leads to the characteristic double-hump shape of the ion line where the peak spectral power does not directly correspond to a resonant or non-resonant wave frequency.

Despite the peak ion line power having no relation to the ion-acoustic frequency, we can still define the ion-acoustic frequency as either a root to Re ϵ = 0 if it exists or the frequency where ϵ ω 2 is at a minimum. Figure 5 shows the calculation of the ion-acoustic frequency with this definition as a function of electron and ion temperature. The behavior of the ion-acoustic mode cleanly transitions from the resonant to non-resonant cases when T e = 3.5   T i .

Gyro lines in Thomson scatter experiments are typically associated with the electrostatic whistler mode. The whistler mode is inherently magnetized and propagates via the electrons’ gyro motion around the magnetic field. The relatively low power of the gyro line compared to the plasma and ion lines has led to few observations of it—mostly by the Arecibo Observatory (Bhatt et al., 2006; Janches and Nicolls, 2007; Hysell et al., 2017) and the European Incoherent Scatter (EISCAT) radar (Malnes et al., 1993). The gyro line has remained an enigma within the ionospheric radar community due to its limited observations and the complicated magnetized terms in the dielectric function. Hysell et al. (2017) provides a thorough examination of the resulting whistler mode dispersion relation, concluding that a simple formula for the gyro line frequency does not exist.

The standard theory for the gyro line frequency ω G L makes the following harsh assumptions (Hysell et al., 2017): k ∥ 2 v t h , e 2 ≪ k ⊥ 2 v t h , e 2 ≪ ω G L 2 ≪ Ω c e 2 ≪ ω p e 2 (39)

Using Appendix Equation A6 for the magnetized electron susceptibility and neglecting ions, the dielectric function is ϵ = 1 + α 2 1 − ω k ∥ v t h , e ∑ n e − k ⊥ 2 ρ ¯ e 2 I n k ⊥ 2 ρ ¯ e 2 2 Daw y e n + i π e − y e n 2 (40)

With the harsh assumptions of Equation 39, roots to the real part of Equation 40 can be obtained through the following steps: 1) assuming that k ⊥ 2 ρ ¯ e 2 ≪ 1 means I 0 k ⊥ 2 ρ ¯ e 2 ≈ 1 and I n ≠ 0 k ⊥ 2 ρ ¯ e 2 ≈ 0 , so the n ≠ 0 terms are dropped from the summation; 2) Taylor-expand the remaining Bessel function in the small argument limit; 3) Taylor-expand the Dawson function in the large argument limit; 4) retain only first order terms in both expansions; 5) solve for ω G L 2 using the quadratic equation. With these steps and a few minor approximations detailed in Hysell et al. (2017), the gyro line frequency is found to be ω G L 2 = Ω c e 2 ⁡ cos 2 ⁡ θ 1 + 2 k 2 v t h , e 2 Ω c e 2 Ω c e 2 ω p e 2 + 1 − 3 4 k 2 v t h , e 2 Ω c e 2 (41)

From Equation 39 we have assumed Ω c e 2 ω p e 2 ≪ 1 and k 2 v t h , e 2 Ω c e 2 ≪ 1 , so both those factors can be neglected to produce the often quoted gyro line frequency of ω G L ≈ Ω c e ⁡ cos ⁡ θ (42)

Note that in this study, the convention for the aspect angle is that θ = 0 ° corresponds to waves propagating parallel to the magnetic field. For radar observations, θ = 0 ° is obtained when the radar line of sight is parallel to the Earth’s magnetic field.

The assumptions in Equation 39 are required for a clean, simple solution for roots of the dielectric function. However, those assumptions are often not justified. At lower altitudes where gyro lines are often observed, both ω p e and Ω c e can be a similar order of magnitude (Bhatt et al., 2008). The constraint of k ∥ 2 v t h , e 2 ≪ k ⊥ 2 v t h , e 2 is not valid for any gyro lines at Arecibo, as it implies k ∥ ≪ k ⊥ , meaning that cos ⁡ θ ≪ sin ⁡ θ , but the aspect angles at Arecibo range from 30° to 60°. Furthermore, while k ⊥ 2 v t h , e 2 ≪ Ω c e 2 is typically justified and means that the argument of the Bessel functions is small, it is often not small enough to justify dropping the n ≠ 0 terms. A more robust assumption is to assume that k ⊥ 2 ρ ¯ e 2 is small enough that only the n = 1 term is comparable to the n = 0 term. This simplifies the dielectric function to ϵ = 1 + α 2 ( 1 − ω k ∥ v t h , e e − k ⊥ 2 ρ ¯ e 2 I 0 k ⊥ 2 ρ ¯ e 2 2 Daw y 0 + i π e − y 0 2 + I 1 k ⊥ 2 ρ ¯ e 2 2 Daw y 1 + i π e − y 1 2 ) , (43)

where y 0 = ω / k ∥ v t h , e and y 1 = ω − Ω c e / k ∥ v t h , e in Equation 43.

The imaginary parts and therefore the damping of the whistler mode are dominated by the terms exp − y 0 2 and exp − y 1 2 . The former describes Landau damping and is important at small frequencies, and the latter describes cyclotron damping at the first gyro-resonance and is maximized when ω ≈ Ω c e . Neglecting the damping components, the normal mode frequency of the gyro line can be obtained by solving for roots to the real part of the dielectric function. Noting that y 1 = y 0 − C e , which defines C e = Ω / k ∥ v t h , e as the normalized gyrofrequency, we can then solve for the roots of 0 = 1 + α 2 1 − 2 y 0 e − k ⊥ 2 ρ ¯ e 2 I 0 k ⊥ 2 ρ ¯ e 2 Daw y 0 + I 1 k ⊥ 2 ρ ¯ e 2 Daw y 0 − C e (44)

As with the plasma line, the dielectric function at ω = 0 is strictly positive, so a root will only exist if, for some non-zero frequency, the dielectric function is negative. We can write as the inequality 1 α 2 + 1 < 2 y 0 e − k ⊥ 2 ρ ¯ e 2 I 0 k ⊥ 2 ρ ¯ e 2 Daw y 0 + I 1 k ⊥ 2 ρ ¯ e 2 Daw y 0 − C e (45)

The problem in solving either Equations 44 or 45, is that the Bessel and Dawson functions are transcendental, and a general solution is not tractable unless the assumptions of Equation 39 are made to justify Taylor expansions. It is therefore not possible to obtain an analytical solution for the gyro line frequency or conditions for its existence unless the approximations in Equation 39 are used.

The existing gyro line theory in Equations 41 and 42 relies on a narrow set of assumptions needed to simplify the magnetized dielectric function. The primary difficulty in a general solution is the presence of the infinite summation of the Bessel functions, the argument of which is called the “finite Larmor radius parameter” and is defined in Equation 46 as: b ≡ k ⊥ 2 ρ ¯ e 2 = k ⊥ 2   v t h , e 2 2 Ω c s 2 (46)

For infinitesimal b, only the n = 0 term is needed for the dielectric function. However, the n = 1 term can become important even when b is as small as 0.03 (Figure 6). As the finite Larmor radius parameter increases, higher order terms in the summation in the dielectric function are needed. These higher order terms can either remove roots from the dielectric function (for intermediate aspect angles) or further complicate the problem by creating even more roots that correspond to the magnetized Berstein modes (aspect angles near 90°).

Figure 6 shows the gyro line’s dependence on electron temperature and therefore on b. At low temperature, the finite Larmor radius parameter is b = 0.035 , and a normal mode solution is clearly present even though the n = 1 term contributes to the dielectric function. Visually, the importance of the n = 1 term can be assessed by seeing the substantial increase in cyclotron damping ( Im ϵ ) when ω ≈ Ω c e . Increasing the temperature in Figure 6 shows that the normal mode resonance is lost when b = 0.071 , but a non-resonant mode is still obtained by finding the minima of the dielectric function. Further increasing the temperature leads to b = 0.106 , allowing the higher-order terms in the summation to wash out the root typically produced by the n = 0 term. At the higher temperature of 1500 K, the minimum of the dielectric function disappears, but a vestigial gyro line is still present in the scattering spectra. This interesting feature will be further discussed in Section 4.3.

The finite Larmor radius parameter can be minimized by either smaller temperatures, larger magnetic fields, or smaller wavenumbers. Note that changing the aspect angle will change b as well, but the tradeoff is that the argument ω / k ∥ v t h , e also changes and will modify the location of the roots and damping. While the magnetic field varies slightly with altitude in the ionosphere, both k ⊥ and Ω c e are primarily dictated by experimental setup. To investigate this dependence, the gyro line for a 230-MHz radar is calculated in Figure 7 (results are applicable to 224 MHz and 233 MHz EISCAT radars in Table 1). The Bragg scatter wavenumber for a 230-MHz radar is k = 9.64 , compared to Arecibo’s k = 18.02 at 430 MHz. Additionally, for EISCAT’s location in northern Scandinavia, the magnetic field at ∼200 km is 4.92 × 10 − 5   T , compared to 3.36 × 10 − 5   T at Arecibo. Both of these conditions lead to smaller finite Larmor radius parameters than at EISCAT compared to Arecibo for a given temperature. The smaller finite Larmor radius parameter leads to most of the assumptions in Equation 39 being valid, so the resulting whistler mode is a normal mode of the plasma with minimal damping. Figure 7 shows that the resulting gyro lines at EISCAT are considerably sharper and more powerful than the gyro lines at Arecibo.

The dependence of the gyro line frequency on plasma parameters is investigated in Figure 8 for Arecibo and Figure 9 for EISCAT. In both figures, the aspect angle is fixed at 45°, with plots at different aspect angles shown in the supporting information. In panel (a) of each figure, the gyro line frequency is obtained by solving for the minima of ϵ 2 . This frequency is compared to the frequency found from solving roots of Re ϵ = 0 in panel (b). For EISCAT, the roots are distinct and easy to obtain, so there is little difference between the two resulting frequencies. However, for Arecibo, the root is not present across a wide range of typical plasma parameters, and therefore the resulting gyro line is associated with non-resonant whistler waves. In panel (c) of Figures 8 and 9, the power at the gyro line frequency is calculated and compared to the ion line power. The ion line power is calculated analytically for ω = 0 as (Froula et al., 2011) S ω = 0 , k → = 2 π k ∥ v t h , i α 2 1 + α 2 1 + T e T i 2 (47)

Since the gyro line is not influenced by ion dynamics, it is assumed that T e = T i .

The estimate in panel (c) shows how easily the gyro line could detect relative to the ion line. However, as the electron temperature increases, the gyro line experiences more Landau and cyclotron damping, broadening the spectral peak. Eventually, for high enough temperatures, the whistler mode becomes non-resonant and decreases in power while broadening substantially (Figure 6). This could lead to experimental difficulties in detecting the gyro line peak within a noisy measurement of the scattering spectra. Panel (d) in Figures 8, 9 estimates the relative prominence of the gyro line peak by calculating the power at ω G L and 0.9   ω G L and plotting the ratio of the power in dB, 10 ⁡ log 10 S ω G L S 0.9 ω G L . While the choice of 0.9 is somewhat arbitrary, it does provide an indication of how prominent the gyro line peak will be compared to the broader noise-dominated spectrum.

Solving for the gyro line frequency shows that the minima finding technique has two significant advantages over the typical root solving in normal mode analysis. First, the gyro line frequency can be calculated across a wider range of plasma parameters, better aligning with them where gyro lines are observed. Second, when a root does exist, it is significantly easier to find it with a bracketing method that searches for it near the non-resonant frequency. The root solving in this paper bracketed the root between 0.8   ω G L , ⁡ min and 1.2   ω G L , ⁡ min , with ω G L , ⁡ min being found by minimizing ϵ 2 . This proved to be a robust root solving algorithm that always found the correct root with no convergence issues.

Figures 8, 9, along with the similar figures in the Supplementary Material, provide a full range of conditions needed for a radar to observe gyro lines at Arecibo and EISCAT VHF/3D. The data availability statement provides the code used to generate these figures and can readily create similar figures for different radars to predict the detectability of gyro lines.

Figure 6 plots the gyro lines at Arecibo for different electron temperatures. As the temperature rises, the root to the dielectric function disappears and then the minima of the dielectric disappear. Interestingly, a gyro-line-like feature remained present for each temperature. To examine this more closely, Figure 10 re-plots the same T e = 1500   K case at Arecibo from Figure 6 across a wider frequency range. While the plasma line at ω / k ∥ v t h , e ≈ 7 corresponds to a distinct minimum of ϵ 2 , there are no other minima of the dielectric. However, broad spectral peaks can still be observed at ω / k ∥ v t h , e ≈ 2 and ω / k ∥ v t h , e ≈ 5 . While the peak at ω / k ∥ v t h , e ≈ 2 was referred to as a “vestigial gyro line” in the previous section, there is no such transition from a gyro line for the ω / k ∥ v t h , e ≈ 5 peak. Therefore, these peaks need a more general interpretation.

In Figure 10, the only visible feature of the dielectric function at the vestigial gyro line is an inflection point. However, there is no obvious interpretation for what an inflection point in the dielectric function would physically mean, so we therefore attribute no significance to these inflection points. Furthermore, it has yet to be determined why there is scattering power at any of the other frequencies between the ion and plasma lines. Both the broad “vestigial gyro line” and the broader “shelf” feature between the ion and plasma lines can again be explained by the analogy of the driven oscillator. Previously, we focused on characterizing the plasma’s response to driven oscillations by looking for roots or minima of the dielectric function. The balancing part of this analogy is the source term that generates waves and drives fluctuations in the plasma. This source term is plotted in Figure 10. Again, there are inflection points at the peaks in the scattering power, but it does not appear to be fruitful or physically meaningful to try and characterize inflection points. However, it is clear that the non-constant source term balanced against the dielectric function leads to the bumps in the scattering spectra, as well as the general filling-in of the spectra.

The vestigial gyro lines and shelf features in Figures 3, 6, 7, and 10 are not dictated by the plasma’s response but by the driving of the system by the equilibrium distribution, and therefore we call these features “distribution driven fluctuations”. This choice of terminology reflects the dominant role of the source term in driving the fluctuations and creating a possibly measurable scattering power. For a distribution driven fluctuation to exist, the driving source term must be substantially large and continuously maintained in equilibrium in order to survive the ensemble average. In contrast, the normal modes in a plasma can be driven by an infinitesimal perturbation and still result in high scattering power.

The transition of gyro lines into broad distribution driven fluctuations is one example of non-standard Thomson scatter spectra. Other exotic spectra include perpendicular-to-B ion lines driven by Coulomb collisions (Kudeki and Milla, 2011; Milla and Kudeki, 2011), ion lines distorted by non-Maxwellian distribution functions (Goodwin et al., 2018), and plasma line splitting (Bhatt et al., 2008). In this section, we provide an example of interpreting these types of exotic spectra by examining the roots and minima of the dielectric function for plasma line splitting.

Plasma line splitting is a phenomenon first observed at Arecibo by Bhatt et al. (2008), where two distinct spectral peaks occur near the plasma frequency. This phenomenon was originally proposed in Salpeter (1961), predicting that two roots will appear in the dielectric function when the plasma frequency is near the second harmonic of the gyro frequency ( ω p e ≈ 2 Ω c e ). This splitting is shown in Figure 11, plotting the dielectric function and scattering spectra for several densities. At the lower density ( n e = 2 × 10 10   m − 3 ), the plasma line is sharp and corresponds to a normal mode of the plasma; a similarly sharp plasma line occurs at higher density ( n e = 6 × 10 10   m − 3 ). However, at the chosen intermediate density ( n e = 4 × 10 10   m − 3 ,) the plasma line has two distinct spectral peaks, one of which corresponds to a normal mode with a root to Re ϵ = 0 , and the other peak corresponds to a non-resonant wave where the dielectric function is at a minimum but has no root.

The parameters in Figure 11 show the plasma line occurring at a lower frequency (∼1.5 MHz), then jumping to a higher frequency (∼2.5 MHz), with the plasma line splitting occurring as an intermediate step. To understand this transition further, Figure 12 plots the dielectric function for a narrower set of density values, with the inset showing where roots to Re ϵ = 0 occur. The predicted behavior from Salpeter (1961) occurs, where the single root occurs at lower frequencies and then jumps to higher frequencies. Interestingly, the double root reported in Salpeter (1961) is actually a triple root to the real part of the dielectric ( n e = 3.6 × 10 10   m − 3 curve, corresponding to ω p e / Ω c e = 1.81 ). However, the scattering spectrum only has two peaks because cyclotron damping is strongest at the harmonics of the gyro frequency, and therefore the middle root has no effect on the wave behavior.

For each of the densities shown in Figure 12, the plasma line spectrum has two distinct peaks that sit on top of a broader spectral enhancement. This broader spectral enhancement is another example of distribution-driven fluctuations, where waves are continually excited near the second gyroharmonic. This is seen from evaluating the magnetized source term (Appendix Equation B4) with ω = 2 Ω c e , and seeing that exp ω − 2 Ω c e k ∥ v t h , e is maximized at this condition. The plasma line spectrum in this regime is therefore a balance between the driving of waves by the electron’s gyro motion and the plasma’s response in this frequency range. The imaginary part of the dielectric shows strong cyclotron damping at the second gyroharmonic, whereas the real part of the dielectric is close to 0 for a broad frequency range. The measurements of plasma line splitting in Bhatt et al. (2008) showed some filling in of the spectrum between the two spectral peaks, but a careful reanalysis of those experiments would need to be done to rule out instrumental or signal processing effects.

This study has separately examined the dielectric function for the plasma line, the ion line, and the gyro line. These are common names for the spectral features observed in ionospheric Thomson scatter experiments, but as we have shown, the underlying wave mode or fluctuation may have a different physical origin depending on the plasma and radar parameters. Table 2 consolidates the terminology used to describe these different types of waves and fluctuations, the required conditions for that type of fluctuation to be present, and the corresponding spectral features. The usage of this terminology for Thomson scatter experiments is demonstrated in Figure 13, which revisits the sample Arecibo spectra plotted in Figure 1.

Table 2 also highlights a major problem within the ionospheric radar community: every measurement is erroneously called “incoherent scatter.” The original idea of ionospheric radar was posited in Gordon (1958) and assumed that electrons in the ionosphere would be randomly distributed, and therefore the phases of scattered waves would be random and the total backscatter would be incoherent. The terminology of “incoherent scatter” has persisted despite its well-known inaccuracy. Colloquially, an incoherent scatter radar is any ionospheric radar capable of making routine ion-line measurements with enough sensitivity to fit the ion line for plasma parameters. Formally, these are high-power and large-aperture Thomson scatter radars that operate in the collective scattering regime where α ≫ 1 (Equation 19). When α ≫ 1 , the incident wavelength is larger than the Debye length, and the resulting scatter is off plasma waves. These wave fronts provide enough structure for the Bragg scatter condition to be met, where constructive interference occurs from scatter off successive wavefronts and creates coherence in the backscattered wave.

The distribution-driven fluctuations shown in Figure 10 are an interesting transition case between coherent and incoherent scatter. In terms of α , these fluctuations are well within the collective scatter regime. However, the scatter is weak and largely dictated by the equilibrium distribution. True incoherent scatter ( α < 1 ) is also weak, and the spectra exactly follow the electron distribution. The physical distinction is that true incoherent scatter is Doppler broadening of an incident electromagnetic wave, whereas distribution-driven fluctuations physically represent a forced oscillation at a non-resonant frequency, and therefore the scattered wave will have some degree of coherence.

The goal of this study has been to explain the presence of strong spectral features in Thomson scatter experiments when normal wave modes are not present. The ubiquitous measurements of ion lines in the ionosphere were a motivating puzzle which are now explained as non-resonant ion acoustic waves. Non-resonant waves are defined as frequencies where the magnitude of the dielectric function is at a local minimum. This holds a physical analogy to a driven oscillator, where waves are continuously created by Cherenkov radiation (source term) and the dielectric function characterizes the plasma’s response to continuously driven oscillations. Normal wave modes such as the Langmuir mode are also continuously driven, and their amplitudes are the result of a balance between the damping of the wave (dielectric function) and the driving source.

Our analysis used a specific framework (Froula et al., 2011) for calculating the dielectric function and source terms in a thermal plasma. This framework is ideal for this study as it is based on the plasma kinetic equations, but it suffers deficiencies in modeling collisions with the BGK operator. The more accurate Coulomb collision operators in Kudeki and Milla (2011) and Milla and Kudeki (2011) are required for accurate computations of the ion line at aspect angles within ∼10° of perpendicular to the magnetic field, and possibly the gyro line in the same regime. The ideas developed here can be generalized to this perpendicular-to-B regime by analyzing the dielectric functions from Kudeki and Milla (2011) and Milla and Kudeki (2011). For example, the ion line exactly perpendicular to B is created by collisional diffusion across magnetic field lines (Milla and Kudeki, 2011) and is best classified as a distribution driven fluctuation.

For extant radars, EISCAT-3D and EISCAT-VHF are best equipped to observe gyro lines and further explore the transition from normal modes at lower temperatures to quasi-normal or non-resonant wave modes at higher temperatures. Nonetheless, the highest resolution gyro line observations were made at Arecibo (Bhatt et al., 2006; Hysell et al., 2017). Future research will examine archived Arecibo experiments to look for gyro lines that transition from sharp spectral peaks to distribution driven fluctuations which would appear as a broad shelf feature between the ion and plasma lines. The −20 dB or lower power of the shelf feature places it at the edge of Arecibo’s sensitivity, although experiments such as Hagen and Behnke (1976) showed Arecibo to be capable of measuring spectra in the very weak, non-collective regime.