The Generation of 150 km Echoes Through Nonlinear Wave Mode Coupling

A fundamental problem in plasma turbulence is understanding how energy cascades across multiple scales. In this paper, a new weak turbulence theory is developed to explain how energy can be transferred from Langmuir and Upper‐Hybrid waves to ion‐acoustic waves. A kinetic approach is used where the Boltzmann equation is Fourier‐Laplace transformed, and the nonlinear term is retained. A unique feature of this approach is the ability to calculate power spectra at low frequencies, for any wavelength or magnetic aspect angle. The results of this theory explain how the predominant type of 150‐km radar echoes are generated in the ionosphere. First, peaks in the suprathermal electron velocity distribution drive a bump‐on‐tail like instability that excites the Upper‐Hybrid mode. This excited wave then couples nonlinearly to the ion‐acoustic mode, generating the ∼10 dB enhancement observed by radars. This theory also explains why higher frequency radars like ALTAIR do not observe these echoes.


Introduction
Kinetic plasma theory is built on a foundation of linearizing the Boltzmann equation and finding solutions algebraically.This technique has been highly successful in describing kinetic phenomena such as Landau damping (Nicholson, 1983) and instability growth rates (Longley et al., 2020).However, linear theory assumes the plasma is near equilibrium, with only small perturbations occurring.This assumption makes it impossible for different wave modes to interact with each other, limiting its application in studies of turbulence and instabilities.This paper develops a nonlinear kinetic theory to explain how two different wave modes can interact.
The motivation for this study is to explain the ionospheric phenomenon called "150-km echoes".These echoes are meter-scale irregularities observed consistently during the day by the Jicamarca radar in Peru, as well as at other equatorial radars.150-km echoes have a typical necklace shape (inset of Figure 1), where a 10-20 dB enhancement of the ion-acoustic mode is observed to follow electron density contours between 130 and 180 km in altitude during the day.Chau et al. (2023) provides a review of the observational history of 150-km echoes and the open questions surrounding them.The biggest open question is the simplest: How are 150-km echoes generated?Some recent progress has been made toward understanding the generation mechanism.Oppenheim and Dimant (2016), Longley et al. (2020Longley et al. ( , 2021)), and Green et al. (2023) developed theory and simulations to show how peaks in the photoelectron distribution can drive a bump-on-tail-like instability at the same locations that 150-km echoes are observed, but those studies did not provide a complete explanation of their generation mechanism.
Figure 1 illustrates the currently incomplete knowledge of how 150-km echoes may be generated.The growth rates of the photoelectron driven Upper-Hybrid instability (Longley et al., 2020) match several features of the observed echoes, most notably the gaps at specific density contours.Figure 2 demonstrates that the Upper-Hybrid mode excited by the photoelectron driven instability occurs at high frequencies (∼4 MHz) and short wavelengths (λ ≈ 20 cm).However, 150-km echoes are observed at the much lower ion-acoustic mode frequencies (∼kHz), with wavelengths of 3-m or larger.Therefore, it must be explained how the energy from the instability converts to energy in the observed ion-acoustic mode.This model does that.Longley et al. (2021).At an angle almost perpendicular to the magnetic field, the Upper-Hybrid mode is the most prominent feature at high frequencies (∼4 MHz).The Upper-Hybrid mode is driven unstable by peaks in the photoelectron distribution.However, 150 km echoes are observed at lower frequencies (∼kHz) and larger wavelengths.The theory developed in this paper describes how the instability at high frequencies couples nonlinearly to the low frequency waves (blue arrow).
Figure 1.The growth rate of the photoelectron driven Upper-Hybrid instability is calculated over a full day (Longley et al., 2020), and explains many features of 150-km echoes.The blue curves follow density contours where ω UH = nΩ ce , where n is the labeled integer.A typical observation of 150-km echoes and their characteristic necklace shape is shown in the inset (adapted from Chau & Kudeki, 2013).
The goal of this paper is to describe how the high frequency instability developed in Longley et al. (2020Longley et al. ( , 2021) ) can couple nonlinearly to the low frequency ion-acoustic mode.Section 2 describes the theoretical derivation of this nonlinear mode coupling.Section 3 shows how it is applied to the problem of 150-km echoes, giving a complete explanation of how the predominant "naturally enhanced incoherent scatter" type of echoes are generated.Section 4 describes the broader applications of this theoretical approach.

Linear Versus Nonlinear Theory
To understand the need for a nonlinear solution of the Boltzmann equation, we first start with a typical linear solution.The Boltzmann equation is a general description of kinetic plasma behavior: Physically, the Boltzmann equation is a total time derivative of the distribution function F s [t, x → , v → ] on the lefthand side, with a collision operator S on the right-hand side.In this paper the convention is used such that The process of linearization is to assume each variable quantity can be decomposed into a 0th order term, and a small first order perturbation: We assume that E 0 = 0 for this solution.In the ionosphere, the Earth's magnetic field is strong enough that the electrostatic approximation is made, taking B → 1 = 0. Putting Equations 2-4 into Equation 1 produces All the 0th order terms can be collected together: The zeroth order terms are simple and setting f 0s to an unchanging Maxwellian leads to the first set of brackets evaluating to 0. This leaves Typically, the next step in solving for density perturbations is to linearize Equation 7 by keeping only the first order terms (the bracketed terms).This linearization is justified by choosing the zeroth order distribution such that the first order distribution is a small perturbation, and since E 1 ∝ f 1s from Gauss' law, the term E → 1 • ∇ v f 1s is the product of two small quantities, and therefore is an even smaller quantity.Thus, only the terms in the brackets are retained, and the resulting equation is linear as it only contains single factors of f 1s or E 1 .
The basis of the mode coupling developed here is to retain the nonlinear term, E → 1 • ∇ v f 1s , in Equation 7. In the case of instability, the perturbed electric field can have significant amplitude and therefore the nonlinear term is no longer a small quantity.

The Fourier-Laplace Transform
We continue to follow the typical linear approach and Fourier-Laplace transform Equation 7, including the nonlinear term.In this paper we use the non-unitary convention for the Fourier transform, defined as As discussed in texts such as Nicholson (1983), the time variable must be Laplace transformed instead of Fourier transformed in order to accommodate Landau damping.The convention we use for the Laplace transform is In applying the Fourier-Laplace transform to Equation 7, the linear terms are straightforward and transform as In Equation 10 an initial value term results from the Laplace transform, as an integration by parts is used to evaluate the time derivative.
The Fourier-Laplace transforms in Equations 10-13 are straightforward because each term is linear, and therefore only one function of t or x is present.However, the nonlinear term in Equation 7 is a product of two functions that depend on both space and time.The Fourier transform of a product of functions is a convolution: Geophysical Research Letters 10.1029/2023GL107212 LONGLEY where ⋆ is the convolution operator.The convolution can be worked out to the single integral Notice that the Fourier transform of the product results in an integration over a dummy wavenumber κ, where one function is evaluated at κ and the other is evaluated at k κ.This integration over all wavenumbers (and later frequencies) allows different wave modes to interact nonlinearly.The Laplace transform also results in a convolution integral with the argument ω ν, where ν is the dummy integration frequency.
At this point a collision operator must be specified.Physically, Coulomb collisions are an important process, however an analytic solution for the Fokker-Planck collision operator is not supported in this framework (see Kudeki & Milla, 2011;Longley et al., 2019;Milla & Kudeki, 2011).Instead, we use the BGK operator as it has the advantage of being linear, and it adequately describes electron-neutral collisions: where ν s is the collision rate, and the zeroth and first order densities relate to the distribution functions as The full Fourier-Laplace transform of Equation 7 is then (17)

Recursively Defined Equations: Breaking Self-Consistency
Thomson scatter radars such as Jicamarca work by transmitting radio pulses into the ionosphere, then measuring the weak backscattered signal.Physically, the radio waves Thomson scatter off electrons, but a strong signal is only detected if a Bragg scattering condition is met.In plasma, this structuring comes from waves or instability fluctuations.Froula et al. (2011) derives the scattering power of a collisional plasma as where ν q is the collision frequency for the species being ensemble averaged.The quantity S(ω, k → ) is analogous to a scattering cross section in the radar equation: Equation 19 relates P rec , the power received by the radar, to constants (r 0 is the classical electron radius) and system parameters (transmitted power P tr and antenna area A).The remaining term S(ω, k → ) therefore contains all Geophysical Research Letters 10.1029/2023GL107212 the information of the plasma.Since the scattering power S depends directly on the electron density fluctuations n 1e , we can also use the final results to study density fluctuations and waves more generally.
To obtain the scattering spectra S(ω, k → ) in Equation 18, the perturbed electron density is solved for and then ensemble averaged.This is done by taking a system of equations consisting of the Boltzmann equation for electrons (Equation 17), the Boltzmann equation for ions (linearized version of Equation 17), and Gauss' law.The Fourier-Laplace transform of Gauss' law is The charge density therefore couples the electron and ion Boltzmann equations since If Gauss' law is directly substituted into Equation 17 and then integrated over velocity, the perturbed densities are The collisional integral U s , susceptibilities χ s ≡ H s /(1 + U S ) and the distribution-like terms M s ≡ 2ν s 〈|D s |〉 are derived in Froula et al. (2011).The nonlinear term defines the shorthand function: The system of Equations 21 and 22 presents a significant and insurmountable problem: the solution is a recursively defined equation.Take for example, the test case where only electrons are present (n 1i = 0) , Equation 21reduces to the form On the left-hand side every function is evaluated at (ω, k → ) , which is the desired frequency of the ion-acoustic mode with the wavenumber set by the radar's transmitted frequency.However, the right hand side involves the exact quantity that is trying to be solved for at (ω, k → ) , as well as There are no standard methods for solving recursively defined equations analytically, even for the simplest case of f(x) + f(x c) = a.A numeric solution is possible through a recursive algorithm, but such a solution is very sensitive to the initial value used to start the solution.
Physically, the resulting recursive equation makes sense.To know the density fluctuations at (ω,k), one must know the fluctuations at (ω ν,k κ) and (ν,κ).But by coupling the modes together, the density fluctuations at (ω, k) will change the fluctuations at (ω ν,k κ), which then changes the fluctuations at (ω,k) and so on.A general solution of Equations 21 and 22 is likely impossible.However, for the case of 150-km echoes we are interested in how the high-frequency Upper-Hybrid mode affects the low-frequency ion-acoustic mode.This is fortuitous, as only electron motions are relevant for the Upper-Hybrid mode.We make the argument that the low-frequency motions from the ion-acoustic mode do not affect the high-frequency Upper-Hybrid (and Langmuir) mode at all.This argument is based on the vastly different masses of electrons and ions and is used often in plasma physics (e.g., only electrons contribute to Debye shielding).
To gain tractability, the recursively defined system of Equations 21 and 22 must be separated by time scale.The nonlinear term is assumed to only be important for high frequencies-specifically the Upper-Hybrid frequency ω hf .
The perturbed density and distribution in the nonlinear term are then taken to be a separate function than , and these separate functions are denoted as The system of equations then becomes where G hf s is Equation 23 evaluated with f hf 1s .Now in Equations 25 and 26, the left-hand side is entirely functions evaluated at (ω, k → ) , and the right-hand side is the nonlinear term with the high-frequency (hf) functions being effectively constants when solving for For the problem of 150-km echoes, photoelectron peaks will drive the Upper-Hybrid mode unstable.The first attempt at solving Equations 25 and 26 calculated the high-frequency terms (n hf 1s and f hf 1s ) using the linear theory for Upper-Hybrid mode generation developed in Longley et al. (2021).This solution failed because the Longley et al. (2021) theory is linear, and therefore does not adequately describe the wave growth during an instability.The problem came from the substitution of the perturbed electric field E 1 using Gauss' law.
Reversing the substitution of Gauss' law only for the nonlinear term, the system of equations is To solve this system of equations we will assume a value for the perturbed electric field.This is similar to how quasilinear theory requires an external evaluation of the wave amplitude (Nicholson, 1983).The choice of electric field value is discussed further in Section 3.
Finally, we recognize that the ion distribution in G i is evaluated at the high-frequency mode.As argued before, the high ion mass will prohibit the ions from moving at high frequencies, and therefore f hf 1i will be zero at ν ≈ ω hf , effectively dropping the nonlinear term from the ion equation.Solving the system of Equations 27 and 28 with G i = 0 finally yields Geophysical Research Letters (29)

Solving for the Scattering Spectra
To obtain the scattering spectra, Equation 29 is squared and ensemble averaged: where CC denotes a complex conjugate, and the angle brackets denote the ensemble average defined in Froula et al. (2011): P(v) is the probability of finding the system in state v, and there is a velocity coordinate in the integral for each particle in the system (i.e., ∫dv 1 dv 2 …dv N ).
In Equation 30, the first two terms define the scattering spectra from linear theory (Froula et al., 2011;Kudeki & Milla, 2011).The functions M s are defined in the Supporting Information, and physically represent the scattering from a non-interacting gas of particles.The factor χ e /ϵ is very large at wave modes where ϵ → 0, and shows how the scattering spectra (and density fluctuations) are strongest when waves are present to create the structuring for Bragg scattering.
The remaining nonlinear term in Equation 30shows that the mode coupling is an additive term to the density fluctuations and can be evaluated independently from the linear terms.The challenge of the nonlinear term is its integration over all frequency and wavenumber space.First, we will simplify the number of integrals by assuming a spherical symmetry to the system, so the wavenumber integrations simplify as The frequency integration is approximated by assuming the integrand is shaped as a Lorentzian over frequency.This is partially justified in Longley et al. (2021), which showed that a Lorentzian describes the frequency dependence of density fluctuations of the Upper-Hybrid mode.Therefore, Effectively, the growth/damping rate γ of the Upper-Hybrid mode determines the width of the Lorentzian.Setting ν = ω hf as the center frequency of the Upper-Hybrid mode also sets κ Putting Equations 32 and 33 into Equation 30, with κ→k hf : Geophysical Research Letters The result involves the multiplication of two single variable integrals.While this double integral follows from squaring the nonlinear term, it lacks an obvious physical meaning.The derivation in this paper makes several approximations that are mostly justified, with the goal of obtaining an analytic solution.As written, Equation 34 is intractable unless the unjustified approximation of k hf,1 = k hf,2 is made.In this step, the integration over k hf is also discretized, and therefore taking k hf,1 = k hf,2 is equivalent to taking the diagonal terms in the product of summations: The electric field will be externally imposed, and therefore does not depend on the ensemble average and is factored out.Substituting back for G e yields where we define the shifted frequencies and wavenumbers as with (ω, k → ) being low-frequency values.

Summary of Solution
Equation 36 is the solution for the electron density fluctuations at low frequencies.It includes contributions from all possible driving wavenumbers, k hf .Since the integration over k hf was discretized, we can either evaluate the contributions from a single driving source or add the contributions from multiple sources to obtain a realistic result of the Upper-Hybrid mode being excited across a wide range of wavenumbers.
To obtain the scattering power, Equation 18 is combined with Equation 36: The ensemble average term in Equation 39 is defined as The remaining task is to evaluate Equation 40.Note that in Equation 40 there are three different frequency/ wavenumber combinations: the low frequency mode (ω, k → ) that came from solving for the density fluctuations; the shifted wavenumber k ), and the high frequency mode quencies and driven by photoelectrons, the linear theory from Longley et al. ( 2021) is appropriate for calculating f hf 1s as the high frequency response of electrons to perturbations.This term from Longley et al. ( 2021) is then where the functions in the leading factor are The results contain velocity integrations over multiple poles that must be dealt with.For example, the first term in the expanded ensemble average is Typically, linear solutions will only result in one or two simple poles in the denominator, and therefore the Plemelj theorem can be applied as an analytic solution for the v ∥ integral.However, the Plemelj theorem results in principal value integrals that are not easy for either non-Maxwellian distributions or multiple poles.For this reason, the integrations are evaluated numerically over v ∥ and v ⊥ (the integration over ϕ is already performed, resulting in 2π due to cylindrical symmetry).The poles lead to difficulties in the v ∥ integration, though working in the collisional regime means the pole is displaced away from the real axis.A numerical algorithm was constructed to evaluate the v ∥ integrals along the real line by choosing a velocity mesh with high resolution centered at each pole.
At this point, the physical derivation is complete.The mathematical solution for Equations 40 and 41 is lengthy but straight forward, so the results are listed in the Supporting Information.Readers interested in a step-by-step calculation of these terms can consult the materials in the "Open Research" statement of this manuscript, which also includes the code used to calculate spectra from Equation 39.

Application to 150-km Echoes
The nonlinear mode coupling developed in Section 2 applies to any set of waves that are substantially separated in frequency, meaning ω ≪ ω hf .This mode coupling framework explains how the photoelectron driven instability in Longley et al. (2020Longley et al. ( , 2021) ) can generate 150-km echoes.The generation mechanism is as follows: 1. Peaks are created in the photoelectron distribution.A collisional resonance with N 2 creates a dip at ∼2.5 eV that manifests as a peak at ∼5 eV.Additionally, EUV emission lines from the Sun create peaks in the range of 10-30 eV.See Oppenheim and Dimant (2016) and Chau et al. (2023).2. The photoelectron peaks drive a bump-on-tail like instability.This instability excites the Langmuir/Upper-Hybrid mode at high frequencies (ω hf ≈ ω pe , typically 5-10 MHz) and short wavelengths (∼20 cm).The linear theory of this instability is derived in Longley et al. (2020Longley et al. ( , 2021)).3.At high enough wave amplitudes, the unstable Upper-Hybrid mode couples nonlinearly to the low frequency ion-acoustic mode (ω ≈ ω ia , typically ∼kHz) through the mechanism derived in Section 2.
Observationally, 150-km echoes are a 10-20 dB enhancement of the ion-acoustic mode.They are observed where photoelectron production is the highest (130-180 km), and only in directions nearly perpendicular to Earth's magnetic field.To calculate the scattering spectra in Equation 39, we specify the zeroth order distribution as the summation of a thermal Maxwellian population, and a photoelectron population with a power-law tail and bumpon-tail features at V i = 5, 15, 25, and 45 eV.This is the same photoelectron distribution used in Longley et al. (2021) and is advantageous as it allows derivatives to be calculated analytically.
The theory in Section 2 requires the electric field amplitude of the high-frequency mode E 1 [ω * ,k * → ] to be specified externally.Determining the saturated amplitude of an instability is an open problem in plasma physics.Therefore, this remains a free parameter in the calculation.In Derghazarian et al. (2023) a similar problem of cross-wavelength coupling between lower-hybrid waves was studied.In that study, a range of wave amplitudes from 20 μV/m to 20 mV/m is used based on in situ measurements.Since the scattering power in Section 2 varies as S ∝ |E 1 | 2 , we do not need to investigate a range of wave amplitudes as each order of magnitude increase results in a 20 dB power increase.We make the choice of using a 1 μV/m electric field, which is an order of magnitude estimate well below the wave amplitudes studied in Derghazarian et al. (2023), and it produces realistic results.
Figure 3 shows the calculated ion-acoustic mode spectra for a single driving wavenumber of k hf = 24, adding the contributions of both the upshifted (positive k hf , wave traveling to radar) and downshifted (negative k hf ) Upper-Hybrid mode.The enhancement of the nonlinear spectra (S NL ) is calculated in dB as where the linear spectra are calculated from Kudeki and Milla (2011).In several panels of Figure 3, a noticeable left-right asymmetry is present in the nonlinear spectra.Mathematically, this results from the dependence on the resonant velocities y lf n = (ω nΩ ce iν e )/ k ∥ .Physically, this could be a preferential coupling between wave modes traveling in one direction versus another, but a full explanation of this asymmetry is currently not available.Nonetheless, Figure 7 in Chau et al. (2023) shows an asymmetry is present in the data.
The results in Figure 3 show that the level of the enhancement depends strongly on the radar wavenumber and the plasma density.The sensitivity to plasma density, through the ratio ω UH / Ω ce = ̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅ ̅ ω 2 pe + Ω 2 ce √ /Ω ce , is expected as Lehmacher et al. (2020) showed experimentally that forbidden gaps in 150-km echoes occur when ω UH /Ω ce is an integer.Analytic work in Longley et al. (2020) showed this same condition leads to the suppression of the instability by thermal Landau and cyclotron damping.
Observationally, 150-km echoes have only been observed by radars with transmit frequencies in the range of 30-50 MHz (Kudeki & Fawcett, 1993;Patra, Chaitanya, Rao, & Kamaraj, 2020).Notably, the ALTAIR radar (150 MHz) has never observed 150-km echoes despite its favorable location at the equator (Chau et al., 2023).In Figure 3, the enhancement is weakest at 150 MHz, but still present for the ω UH /Ω ce = 3.84 condition.Equation 39shows the imposed wave electric field is a function of radar wavenumber, E 1 [ω * ,k * → ] , but the results in Figure 3 use the same E 1 for each k* ≡ k k hf value.At higher radar wavenumbers, E 1 [ω * ,k * → ] is evaluated further away from the dispersion relation at (ω hf ,k hf ), and therefore the driving electric field for a 150 MHz radar should be smaller than the driving electric field for 30-50 MHz radars.This means that ALTAIR has not observed 150-km echoes due to its frequency being too high, which both inhibits the wave coupling and leads to lower values evaluated away from the dispersion relation.

Summary
The generation mechanism of 150-km echoes is now a solved problem.The nonlinear theory developed in this paper bridges the gap between decades of observations and the linear instability work in Longley et al. (2020Longley et al. ( , 2021)).As discussed in the review paper Chau et al. (2023), there are numerous outstanding questions with 150km echoes which this paper cannot explain.However, understanding the generation mechanism is the first step to understanding more complex details in the echoes.Furthermore, this work only explains the more common but weaker type of echo termed "naturally enhanced incoherent scatter" (NEIS) and does not apply to the stronger but less common "field aligned irregularities" (FAI), though the two types are often observed simultaneously (Chau & Kudeki, 2013;Patra, Chaitanya, Rao, & Reddy, 2020).Clearly, more work is needed to fully understand 150km echoes and their unique combination of plasma dynamics, neutral dynamics, and photochemistry.
The nonlinear mode coupling described in this paper works for any set of electrostatic waves with a separation of frequencies.The use of a driving wave electric field means this theory can be applied to HF heating problems with minimal modifications.In Derghazarian et al. (2021Derghazarian et al. ( , 2023) ) high altitude (∼2,000 km) echoes are observed by Jicamarca with strong enhancements in the lower-hybrid mode.An initial search was done to see if the lowerhybrid mode is excited by the photoelectron instability, but no such enhancements were found.
The use of an externally determined wave electric field was necessary but leaves an open question of what value should be used.The theory in this paper could be fit to the observed power in 150-km echoes, giving an experimental estimate of what the wave amplitude should be.Additionally, the method in this paper of breaking the self-consistency of the nonlinear term can be applied to the photoelectron instability theory, possibly determining the wave amplitude from theory.Both approaches would lead to new insights into how instabilities reach a saturated state.

Figure 2 .
Figure2.The kinetic dispersion relation using the theory inLongley et al. (2021).At an angle almost perpendicular to the magnetic field, the Upper-Hybrid mode is the most prominent feature at high frequencies (∼4 MHz).The Upper-Hybrid mode is driven unstable by peaks in the photoelectron distribution.However, 150 km echoes are observed at lower frequencies (∼kHz) and larger wavelengths.The theory developed in this paper describes how the instability at high frequencies couples nonlinearly to the low frequency waves (blue arrow).

Figure 3 .
Figure 3. Calculated ion-acoustic mode spectra, for different radars and plasma densities.Panels (a) and (b) show the results for a 30 MHz radar (5-m Bragg wavelength, k = 1.26), panels (c) and (d) show the results for a 50 MHz radar (3-m wavelength, k = 2.1), and panels (e) and (f) are for a 150 MHz radar (1-m wavelength, k = 6.3).The density is chosen such that ω UH = 3.84Ω ce in panels (a), (c), and (e), and ω UH = 4.0Ω ce in panels (b), (d), and (f).In each panel, the spectra are normalized to the maximum value of the linear theory (solid blue curve).