# Harmonic Radio Emission in Randomly Inhomogeneous Plasma

Anna Tkachenko, Vladimir Krasnoselskikh, Andrii Voshchepynets

► To cite this version:

Anna Tkachenko, Vladimir Krasnoselskikh, Andrii Voshchepynets. Harmonic Radio Emission in Randomly Inhomogeneous Plasma. The Astrophysical Journal, 2021, 908 (2), pp.126. 10.3847/1538-4357/abd2bd. insu-03169369

HAL Id: insu-03169369

https://insu.hal.science/insu-03169369v1

Submitted on 27 Sep 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Harmonic Radio Emission in Randomly Inhomogeneous Plasma

Anna Tkachenko¹ Ⓓ, Vladimir Krasnoselskikh¹,² Ⓓ, and Andrii Voshchepynets³,⁴ Ⓓ

¹ LPC2E, CNRS and University of Orléans, 3A Avenue de la Recherche Scientifique, Orléans, France; anna.tkachenko@cnrs-orleans.fr

² Space Sciences Laboratory, University of California at Berkeley, 7 Gauss Way, Berkeley, CA 94720, USA

³ The Swedish Institute of Space Physics (IRF), Rymdcampus1, Kiruna, Sweden

⁴ IRAP, 9 Avenue du Colonel Roche, Toulouse, France

Received 2020 September 15; revised 2020 November 2; accepted 2020 November 28; published 2021 February 18

## Abstract

In the present paper, we describe a theoretical model of the generation of harmonic emissions of type III solar radio bursts. The goal of our study is to fully take into account the most efficient physical processes involved in the generation of harmonic electromagnetic emission via nonlinear coupling of Langmuir waves in randomly inhomogeneous plasma of solar wind ($l + l' \rightarrow t$). We revisit the conventional mechanism of coalescence of primarily generated and back-scattered Langmuir waves in quasihomogeneous plasma. Additionally, we propose and investigate another mechanism that generates harmonic emission only in a strongly inhomogeneous plasma: the nonlinear coupling of incident and reflected Langmuir waves inside localized regions with enhanced plasma density (clumps), in the close vicinity of the reflection point. Both mechanisms imply the presence of strong density fluctuations in plasma. We use the results of a probabilistic model of beam-plasma interaction and evaluate the efficiency of energy transfer from Langmuir waves to harmonic emission. We infer that harmonic emissions from a quasihomogeneous plasma are significantly more intense than found in previous studies. The efficiency of Langmuir wave conversion into electromagnetic harmonic emission is expected to be higher at large heliospheric distances for the mechanism operating in quasihomogeneous plasma and at small heliocentric distances for the one operating in inhomogeneous plasma. The evaluation of emission intensity in quasihomogeneous plasma may also be applied for type II solar radio bursts. The radiation pattern in both cases is quadrupolar, and we show that emission from density clumps may efficiently contribute to the visibility of harmonic radio emission.

*Unified Astronomy Thesaurus concepts: Solar wind (1534); Interplanetary turbulence (830); Solar coronal waves (1995)*

## 1. Introduction

Radio emissions in the inner heliosphere, associated with extreme space weather events, such as solar flares, are generated via the plasma emission mechanism first suggested by Ginzburg & Zheleznyakov (1958). Generally, it can be described as consisting of two steps: a beam of electrons, accelerated at reconnection sites of solar flares or excited by coronal mass ejection (CMEs) driven shock waves, propagates from the Sun along opened magnetic field lines. On its way, it interacts with the ambient plasma, generating Langmuir waves (with frequencies close to the local plasma frequency $\omega_{pc}$) via bump-on-tail instability. In turn, these waves may transfer part of their energy into electromagnetic (EM) radio emission at the fundamental frequency (close to the local plasma frequency) and its second harmonic (about twice the plasma frequency). This mechanism is widely recognized as responsible for the generation of type II and type III solar radio bursts. Type IIIs are normally associated with solar flares and energetic electrons beams with typical velocities of 0.1–0.5c, whereas type IIs are attributed to CMEs and less energetic electron beams.

Most type III radio bursts exhibit harmonic structures. Up to decametric wavelengths, the fundamental and second harmonic (hereafter referred to simply as *harmonic*) components can often be distinguished when occurring simultaneously. This wavelength range is supposed to be associated with the coronal plasma. Within the interplanetary medium, at larger wavelengths, it is almost impossible to separate one component from the other, except for rare cases when electron beams are observed in situ by spacecraft (Kellogg 1980; Dulk et al. 1998). There is a big observational base of radio bursts of type III,

made by spacecraft (e.g., Wind, STEREO, Parker Solar Probe (PSP)), as well as by ground-based radio telescopes (e.g., LOFAR, Nançay), which covers various frequency ranges. Type IIIs, due to their mechanism of generation, are a powerful tool for diagnostics of coronal and solar wind plasma, as well as for tracking energetic electron beams (Mann et al. 2018). However, there remain open questions associated with type III radio bursts, and many of them are related to the role played by background density fluctuations of the solar wind and corona (Robinson & Cairns 1998; Reid & Ratcliffe 2014; Chen et al. 2018).

Random density fluctuations within the solar wind and solar corona are a well-known feature of heliospheric plasma (Neugebauer 1975; Goldstein et al. 1995; Shaikh & Zank 2010; Chen et al. 2012). The density spectrum typically exhibits a two-knee power law within the frequency domain of $10^{-3}-10^2$ Hz, with a breaking frequency around 0.6 Hz. Part of the spectrum below the break follows quite well the Kolmogorov power law, while the spectral index above the break depends on the conditions in the solar wind. Reported values of the spectral index above the break are typically in the range [-0.91, -0.38] (Celnikier et al. 1987). However there are reports (e.g., by Kellogg & Horbury 2005) that show that the spectrum of density fluctuations can be well described by a single power law with a spectral index of -1.37. Recent observations as well as in situ measurements by the PSP spacecraft confirm that the level of density fluctuations in the solar wind can go up to 7% of the average background density at ~36 $R_⊙$, and analysis based on comparison of Monte Carlo simulations with PSP observations of the decay times of type
III bursts predicts a growth of the level of density fluctuations toward the Sun, reaching up to 20% in the high corona (Krupar et al. 2020). These turbulent structures strongly affect the propagation and observed properties of radio emissions in coronal and solar wind plasmas (Kontar et al. 2019). At the same time, the presence of density fluctuations has a significant impact on the generation of Langmuir waves (Reid & Kontar 2010; Krafft et al. 2013; Bian et al. 2014; Voshchepynets et al. 2015), which are later converted into EM emission. As shown by Voshchepynets & Krasnoselskikh (2015), part of the density spectrum in the range $10^{-2}$–10 Hz is of particular importance for the bump-on-tail instability under typical conditions of the solar wind plasma.

According to the conventional plasma emission mechanism, EM emission at the fundamental frequency is generated due to the Rayleigh scattering of Langmuir waves by plasma thermal ions, whereas harmonic emission is the result of the Raman scattering of Langmuir waves (Ginzburg & Zheleznyakov 1958). However, some of the observed properties of type III radio emissions are not fully explained by the plasma emission mechanism. This has led to numerous revisions of the initial theory (e.g., see the review by Reid & Ratcliffe 2014). Several mechanisms have been proposed, aiming to cover the properties of fundamental emissions—among them the nonlinear wave–wave interaction of Langmuir, ion sound, and EM waves: $l \pm s \to t$ (e.g., see Gurnett & Frank 1978; Melrose 1987). The most important question concerning the efficiency of such a process is the presence of ion sound waves with wavevectors $k_s \approx \pm k_l$ required to satisfy the resonant condition $k_l \pm k_s \to k_l$, where $k_l = \omega_{p_e}/v_b \gg k_t = \omega_{p_e}/c$ and $v_b$ is the characteristic velocity of the beam. It means that either the broad spectrum of ion sound waves should contain the required waves or they should be generated by the decay instability; otherwise there is a high probability that some other process is responsible for the generation of EM emission at the fundamental frequency (Melrose 1987). An ion sound wave generated by the decay instability $l \to l' + s$ in general cannot be directly involved in the generation of EM emission since the secondary Langmuir wave has a wavevector almost opposite to the wavevector of the primary wave; the wavevector of the sound wave is typically too large $k_s \approx -2k_l$. Another mechanism that may explain EM emissions at the fundamental frequency is the linear mode conversion of Langmuir waves directly into EM waves in the presence of an increasing density gradient. The importance of encounters of Langmuir waves with regions with higher density for the generation of EM waves has been pointed out by a number of authors (e.g., see Hinkel-Lipsker et al. 1992 and references therein) and extensively studied by Mjølhus (1983, 1990), Kim et al. (2007, 2008, 2009, 2013), and Schleyer et al. (2013, 2014) for the case of a magnetized plasma. It was shown by Krasnoselskikh et al. (2019) that even a simple reflection of Langmuir waves on a density inhomogeneity can result in an efficient generation of fundamental EM emission. It raises an important question about the possible role of similar localized regions where the reflection process can take place for the generation of harmonic EM emission.

For harmonic emission with frequency of about $2\omega_{p_e}$, the wavevector of the EM wave is $k_t \approx \sqrt{3}\omega_{p_e}/c$, while the primary Langmuir wave's wavevector is $k_l \approx \omega_{p_e}/v_b$. Thus the wavevector of the EM wave is much smaller than the wavevector of

the Langmuir wave: $k_t \ll k_l$. It leads to the conclusion that the two coalescing Langmuir waves that produce an EM wave at $2\omega_{p_e}$ should be almost antiparallel to fulfill momentum conservation $k_l + k_{l'} = k_t$. In the original work by Ginzburg & Zheleznyakov (1958), harmonic emission was attributed to the induced scattering of beam-driven Langmuir waves on thermal ions and the subsequent coalescence of the forward-moving and scattered waves. However, a few decades later Melrose (1980a, 1980b) argued that ion scattering is not efficient enough to be consistent with the brightness temperatures of type IIIs observed in the corona. This has led to the investigation of the role of ion sound waves in producing back-scattered Langmuir waves via the electrostatic decay process $l \to l' + s$ (e.g., see Cairns 1987a and references therein). Based on this idea, Willes et al. (1996) derived the analytical solutions to describe the $l+l' \to t$ process for a broad class of Langmuir wave spectra, including the case of almost antiparallel waves (a head-on approximation). The ion sound waves playing a role in such processes should be generated in highly non-isothermal plasma with the electron temperature $T_e$ much higher than the ion temperature $T_i$ ($T_e \gg T_i$) (Chen et al. 1984). If this condition is not satisfied, the ion sound waves are quickly damped by resonant interactions with the ions by means of Landau damping. Observations indicate that in solar wind typically $T_e \sim T_i$ (e.g., Lin et al. 1986), which leads to the conclusion that electrostatic decay might be insufficient to account for the observed properties of harmonic type III emissions.

Among the other mechanisms suggested to explain harmonic emissions in plasma, there is radiation by localized bunches of Langmuir waves (Galeev & Krasnoselskikh 1976; Brejzman & Pekker 1978; Papadopoulos & Freund 1978; Goldman et al. 1980; Ergun et al. 2008; Malaspina et al. 2012), which considers radiation by nonlinear currents at twice the plasma frequency (antenna-type radiation). Some of these works imply the existence of a strong turbulence. However, in the presence of strong density fluctuations, this approach becomes invalid. Ergun et al. (2008) have suggested another mechanism, based on the idea that a significant fraction of the Langmuir waves are localized as eigenmodes in solar wind density cavities. To enable this mechanism, density irregularities should have the form of density holes to capture Langmuir waves inside them. This is not always the case, since it requires the presence of localized density depressions. More common topological features are density clumps with density gradients of varying steepness. The theory of harmonic emission in inhomogeneous plasma has already been partially considered, for example by Erokhin et al. (1974) in the one-dimensional case. Our goal is to revisit the theory of harmonic emission by taking into account density fluctuations within the solar wind and corona. In this case, the back-scattered Langmuir waves required for the generation of EM emission are produced by the reflection of forward-moving Langmuir waves from density fluctuations. We consider two different cases of harmonic EM emission generation: (1) emission that is produced in the vicinity of the reflection point inside the clump and (2) emission that is produced far from the clump, where the plasma can be considered quasihomogeneous.

Several studies have also investigated the generation of third and higher harmonics, applicable to type II and III radio bursts (Brejzman & Pekker 1978; Cairns 1987b).
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Schematic illustration of the generation of harmonic EM emission in a randomly inhomogeneous plasma. (a) An electron beam propagates from the Sun and generates a spectrum of Langmuir waves ($l$) with frequency $\omega_l$. The background electron plasma density decreases with distance from the Sun (blue line). (b) Langmuir waves are reflected ($l'$) after an encounter with density clumps. (c) Coalescence of two oppositely directed waves $k_1$ and $k_2$ from the spectra of Langmuir waves in quasihomogeneous plasma: the z-axis is directed along the electron beam direction, and the x-axis is an arbitrary perpendicular direction. (d) Coalescence of an incident Langmuir wave with its reflected part inside the conversion region: the z-axis is directed along the density gradient inside the clump, and the x-axis is an arbitrary perpendicular direction. In panels (c) and (d) $\psi$ is the acute angle between the wavevector of a forward-moving (incident) Langmuir wave and the z-axis, and $\theta$ is the acute angle between the wavevector of the harmonic EM wave and the z-axis. Generally, the z-axes in these two different cases do not coincide. However, since the wavevectors of beam-generated Langmuir waves are usually highly aligned with the beam direction and since we consider coalescence in both regions under head-on approximation, we may assume that these axes are roughly equivalent.</figcaption></figure>

The theoretical model of beam-plasma interaction in a plasma with random density fluctuations developed by Voshchepynets et al. (2015) and Voshchepynets & Krasnoselskikh (2015) shows that the reflection process is indeed very important when the level of density fluctuations becomes large enough to overcome the effect of linear dispersion (Kellogg et al. 1999), namely when

$$ \frac{\langle \Delta n \rangle}{n_0} > 3k_l^2 \lambda_D^2 = 3 \frac{k_B T_e}{E_b}, \qquad (1) $$

where $\langle \Delta n \rangle$ is the average level of density fluctuations, $n_0$ is the average background plasma density, $\lambda_D$ is the Debye length, $k_B$ is the Boltzmann constant, and $E_b$ is the characteristic kinetic energy of the beam. It is also shown that for the aforementioned conditions the portion of wave energy carried by reflected waves approaches 50% of the total wave energy, i.e., the energies of the primary and reflected waves are almost equal. Thus, the process of reflection of Langmuir waves on density fluctuations likely plays a crucial role for the generation of EM emissions in quasihomogeneous plasma, even if it is treated as a conventional coupling process that is not directly affected by density fluctuations. On the other hand, in the close vicinity of the reflection points, the electrostatic

fields are known to be enhanced, which may lead to the operation of the antenna-type mechanism in these localized areas.

In the present paper we study the process of generation of harmonic emission of type III solar radio bursts in a randomly inhomogeneous plasma. An electron beam resonantly generates a spectrum of Langmuir waves with a frequency $\omega_l$ very close to the local electron plasma frequency $\omega_{pe}$ and with wavevectors highly aligned with the beam direction (Malaspina & Ergun 2008; Krasnoselskikh et al. 2011). As the average density of the background plasma decreases, random density fluctuations provide local density enhancements (clumps). Langmuir waves may encounter these clumps, and if the electron plasma frequency inside these structures reaches $\omega_l$, the waves will be reflected in the opposite direction, forming a spectrum of backward-moving (reflected) Langmuir waves (see Figure 1(b)). The forward- and backward-moving Langmuir waves may then interact and produce harmonic EM emissions. We will formally distinguish two different regions of such interaction: (1) the quasihomogeneous plasma with the average local electron density $n_0$ located away from localized density perturbations (see Figure 1(c)) and (2) the locally inhomogeneous plasma inside the density clumps, confined between the
start of the positive density gradient and the reflection point, i.e., confined within the conversion region (see Figure 1(d)). For the first case, we evaluate the process of nonlinear coupling of Langmuir waves assuming a mirror-type reflection and a Gaussian spectrum of forward-moving and reflected Langmuir waves. The reflection process is taken into account by means of the coefficient $P_{\text{ref}}$, characterizing the part of energy carried by reflected waves. The coupling process itself is described similarly to the one in homogeneous plasma, assuming that it is not affected by density fluctuations (Section 2). For the second case, we consider the coalescence of a single forward-moving (incident) Langmuir wave with its reflected part inside the conversion region, in the close vicinity of the reflection point. We evaluate the electric fields of these Langmuir waves assuming a linear density gradient. This allows us to obtain perturbations of the density and velocity of the electrons caused by the presence of Langmuir waves, and consequently we may evaluate the excited nonlinear currents at frequencies around $2\omega_{\text{pe}}$ (Section 3). Next, we estimate the energy density of the EM emission produced by the aforementioned nonlinear currents for such a single event of Langmuir wave reflection within a single density clump. In order to obtain the value of the energy density of EM emissions that corresponds to the full spectrum of Langmuir waves and to multiple reflections from density clumps with different amplitudes of density fluctuations and different characteristic scales of density gradients, we average our result over the relevant parameters (Section 4). In both cases, the energy density of EM harmonic emissions is expressed in terms of the energy density of Langmuir waves. To obtain a quantitative evaluation and deduce scaling laws for the dependencies of EM wave intensity, we use the results of a probabilistic model of beam-plasma interaction to describe the generation of Langmuir waves in a randomly inhomogeneous plasma (Section 5). Finally, we estimate the efficiency of conversion of Langmuir waves into harmonic EM emission for both considered regions of interaction (Section 6). The assumptions used for the evaluation of harmonic emission from a quasihomogeneous plasma also allow us to apply the obtained results to type II solar radio bursts.

## 2. Nonlinear Coupling of Langmuir Waves and Generation of EM Emission in a Quasihomogeneous Plasma

Here we describe the process of generation of harmonic EM emission by nonlinear coupling of Langmuir waves in a quasihomogeneous plasma. We consider the coalescence of two oppositely directed Langmuir waves, one of which, $k_1$, is aligned with the beam (it belongs to a spectrum of primary generated waves) and the other, $k_2$, is supposed to be reflected by density fluctuations (it belongs to the spectrum of reflected waves; see Figure 1(c)). Even though density fluctuations play an important role here, providing reflected waves, we examine harmonic wave generation as happening in a quasihomogeneous plasma, i.e., far from localized density perturbations. The process of EM wave generation by the coupling of Langmuir waves is described by the following set of equations (Tsytovich 2012):

$$ \frac{dN_t(k_t)}{dt} = \int \int \frac{d^3k_1 d^3k_2}{(2\pi)^6} w_{ll}^t(k_1, k_2, k_t) [N_l(k_1)N_l(k_2) - N_l(k_1)N_l(k_t) - N_l(k_2)N_l(k_t)], \quad (2) $$

$$ w_{ll}^t(k_1, k_2, k_t) = \frac{(2\pi)^6 e^2 (k_t^2 - k_2^2)^2 [k_1 \times k_2]^2}{32\pi m_e^2 \omega_{\text{pe}} k_t^2 k_2^2} \times \delta(\omega_t - \omega_1 - \omega_2) \delta(k_t - k_1 - k_2), \quad (3) $$

where $N_l(k_1)$, $N_l(k_2)$, and $N_l(k_t)$ are the numbers of wave quanta of the Langmuir and EM waves; $k_1$, $\omega_1$, $k_2$, and $\omega_2$ are the wavevectors and frequencies of the Langmuir waves; $k_t$ and $\omega_t$ are the wavevector and frequency of the EM wave; and $e$ and $m_e$ are the electron charge and mass. It is reasonable to suppose that the number of quanta of the Langmuir waves is much larger than the number of quanta of the EM waves,

$$ N_l(k_1), N_l(k_2) \gg N_l(k_t); \quad (4) $$

then the second and third terms under the integral in Equation (2) may be neglected. The number of quanta is related to the wave energy density via the following formulas:

$$ W_l(k_l) = \omega_l N_l, \quad W_t(k_t) = \omega_t N_t. \quad (5) $$

We recall that Langmuir waves have a frequency $\omega_t$ very close to the local electron plasma frequency $\omega_{\text{pe}}$ of the region where they are excited. The spectrum of Langmuir waves is formed by two processes: direct excitation of Langmuir waves due to the bump-on-tail instability (thus it may be approximated as a Gaussian, centered at the resonant wavevector $k_b$ with width $\Delta k_b$) and reflection of Langmuir waves, assumed to be of the mirror type (i.e., also a Gaussian, centered at $-k_b$ having the same width in the wavevector space). Here we note that a similar approximation was used by Willes et al. (1996), but the results obtained here are very different, as will be discussed later. Thus, the spectrum of Langmuir waves has the following form (for details see Appendix A):

$$ N_l(k_{\parallel}, k_{\perp}) = \frac{1}{\pi^{3/2} \Delta k_b^3} N \left[ (1 - P_{\text{ref}}) \exp\left(-\frac{(k_{\parallel} - k_b)^2 + k_{\perp}^2}{\Delta k_b^2}\right) + P_{\text{ref}} \exp\left(-\frac{(k_{\parallel} + k_b)^2 + k_{\perp}^2}{\Delta k_b^2}\right) \right]. \quad (6) $$

Here we choose the parallel direction as the direction of beam propagation, $N$ is the total number of Langmuir wave quanta, and $P_{\text{ref}}$ is the reflection coefficient that defines the redistribution of wave energy between primary (forward-moving) and reflected waves. In our calculations, we choose the $z$-axis to be directed along the direction of beam propagation, and the $x$-axis is along the second component of the EM wavevector. After the direct but slightly cumbersome calculations presented in Appendix A, one can obtain the following equation for the energy density of harmonic EM emission in wavevector space generated by the aforementioned Langmuir waves:

$$ \frac{dW_t^{\text{homog}}(k_t)}{dt} = \frac{P_{\text{ref}}(1 - P_{\text{ref}})\omega_{\text{pe}}}{192\pi^5} \left(\frac{k_b}{\Delta k_b}\right)^4 \times \frac{k_t^2}{k_b^5} \frac{w_t}{n_0 k_B T_c} W_t \sin^2 \psi \cos^2 \psi, \quad (7) $$

where the superscript “homog” indicates that we perform this estimation for a quasihomogeneous plasma, $W_t$ is the total energy density of Langmuir waves, and the angle $\psi$ is the angle between the vector $k_1$ and the $z$-axis. Let us return to the discussion of the fact that this result is very different from the
one obtained by Willes et al. (1996, their Equation (29)). Willes et al. (1996) performed their calculation supposing that the following inequality is satisfied:

$$ \frac{2k_r k_l \sin \psi \sin \theta}{\Delta k^2} \gg 1, \qquad (8) $$

and their approximate estimation is strongly dependent on this assumption. We carry out exact calculations without this assumption, and the result obtained shows that the major input comes from the region in wavevector space where this inequality is not satisfied. It is worth noting here that there is a quite simple interpretation of this disagreement. The above assumption corresponds to neglecting currents parallel to the direction of propagation of Langmuir waves. On the other hand, the multiplier $\sin^2 \psi \cos^2 \psi$ indicates that the major directions of the emission comprise angles $\pi/4$ and $3\pi/4$ with the electron beam direction, unambiguously showing the quadrupolar character of the emission and the major input of these parallel electric currents. Integrating over $d^3\mathbf{k}_l$ one can evaluate

$$ \frac{dW_l^{\text{homog}}}{dt} = \frac{P_{\text{ref}}(1 - P_{\text{ref}})\omega_{\text{pe}}}{7200\pi^4} \left(\frac{k_b}{\Delta k_b}\right)^4 \frac{\mathbf{k}_l^5}{k_b^5} \frac{W_l}{n_0 k_B T_e} W_l. \quad (9) $$

Observations show that the dynamics of the burst on its initial stage is quite similar to exponential growth from some noise level until the maximum intensity is reached. Thus the Langmuir wave dynamics may be presented in the form

$$ W_l = W_{\text{noise}} \exp(\gamma t). \qquad (10) $$

The characteristic growth factor until instability saturation is supposed to be of the order of $(\gamma t_s) = \Lambda = \ln(n\lambda_D^3)$; thus

$$ \frac{W_l}{W_{\text{noise}}} = \exp\Lambda, \qquad (11) $$

and the saturation occurs at a time $t_s = \Lambda/\gamma$. Consequently the equation for $W_l^{\text{homog}}$ may be solved as

$$ W_l^{\text{homog}}(t) = \frac{P_{\text{ref}}(1 - P_{\text{ref}})\omega_{\text{pe}}}{7200\pi^4\gamma} \left(\frac{k_b}{\Delta k_b}\right)^4 \frac{\mathbf{k}_l^5}{k_b^5} \frac{W_l^2}{n_0 k_B T_e}. \quad (12) $$

Here $\gamma$ may be evaluated as $\gamma_{\text{lin}}/\Lambda$, where $\gamma_{\text{lin}}$ is the linear increment of the bump-on-tail instability in an inhomogeneous plasma with random fluctuations. Parameters $W_l$ and $\gamma_{\text{lin}}$ will be estimated according to a probabilistic model of beam-plasma interaction in a randomly inhomogeneous plasma in Section 5.

### 3. Description of the Fields and Currents in the Vicinity of Reflection Points

Let us now consider the generation of EM waves that may come from localized regions where reflection of Langmuir waves occurs (see Figure 1(d)). In these regions, the field amplitudes and the corresponding currents are strongly affected by inhomogeneity, and here we present a model that intends to describe them precisely. We consider a Langmuir wave of frequency $\omega_l$ that encounters a density clump, formed due to random density fluctuations, inside which the density linearly increases toward the center. The assumption of a linear density gradient inside the clump is made for the sake of simplicity. As the wave enters the clump, the component of its wavevector in

the direction of the density gradient decreases, and in the point where the local plasma frequency reaches the frequency of the Langmuir wave, this wave will undergo a mirror-type reflection. After reflection, the two waves, the incident $l$ and the reflected $l'$, will perturb the electron trajectories, creating variations of the velocity and particle density that produce currents with frequencies close to $2\omega_{\text{pe}}$ in the vicinity of the reflection point. These localized currents represent a source that can generate an EM wave $t$ with a frequency around $2\omega_{\text{pe}}$. Here we shall consider this process in more detail. We choose a density gradient inside the clump that is directed along the $z$-axis, while we assume that along the other two directions, $x$ and $y$, density variations are negligible, allowing us to reduce the number of dimensions and consider our problem only in the $(x, z)$-plane. The generation of harmonic EM waves in this model is only possible for certain values of the angle of incidence $\psi$, since momentum conservation implies constraints on its value, determined by the ratio of the beam velocity to the speed of light: $|\psi_{\text{max}}| \approx \sqrt{3} v_b / 2c$ (see Appendix C).

We begin our calculations with the system of equations for plasma oscillations (Zakharov 1972):

$$ \Delta\Phi_l^{\text{inhom}}(\mathbf{r}, t) = 4\pi e \delta n_l, \qquad (13) $$

$$ \frac{\partial}{\partial t} \delta n_l + \nabla((n_0 + \delta n)\mathbf{v}_e) = 0, \qquad (14) $$

$$ \frac{\partial \mathbf{v}_e}{\partial t} = \frac{e}{m_e} \nabla \Phi_l^{\text{inhom}}(\mathbf{r}, t) - 3 v_T^2 \nabla \frac{\delta n_l}{n_0}, \qquad (15) $$

where

$$ n_e = n_0 + \delta n + \delta n_l. $$

Here $\delta n$ is the given low-frequency plasma inhomogeneity, $\delta n_l$ is the high-frequency density variation caused by Langmuir oscillations, and $\Phi_l^{\text{inhom}}(\mathbf{r}, t)$ is the corresponding high-frequency part of the electrostatic potential. Here and further below, the superscript “inhom” refers to an inhomogeneous plasma inside a density clump.

We assume that in a quasihomogeneous plasma outside the clump, the high-frequency part of the electrostatic potential has the form of a plane wave: $\Phi_l^{\text{homog}}(\mathbf{r}, t) = \Phi_0^{\text{homog}} \exp(-i(\omega_{\text{pe}}t + \eta) + ik_{l,x}x + ik_{l,z}z)$, where $\Phi_0^{\text{homog}}$ is the amplitude and $\eta$ is the phase difference between the wave in the homogeneous plasma and that in the density clump. As we assume that all the parameters vary only along the $z$-axis, the solution for $\Phi_l^{\text{inhom}}(\mathbf{r}, t)$ can be written in the form $\Phi_l^{\text{inhom}}(\mathbf{r}, t) = \Phi_0^{\text{inhom}} \phi(z) \exp(-i\omega_{\text{pe}}t + ik_{l,x}x)$, where $\phi(z)$ is an unknown function of $z$.

Without loss of generality, one can assume that the density gradient within the clump is linear; thus the density inhomogeneity profile has the following form (density starts to increase at $z=0$):

$$ \frac{\delta n}{n_0} = \theta(z) \frac{z}{L}, $$

where $\theta(z)$ is the Heaviside step function and $L$ is a term that can be interpreted as the characteristic scale of the density gradient. By introducing a new dimensionless variable

$$ z = \left( \frac{k_l^2 L^2}{3k_l^2 \lambda_D^2} \right)^{1/3} \left( \frac{z}{L} - 3k_l^2 \lambda_D^2 \cos^2 \psi \right), \quad (16) $$
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Absolute value of total electric field of incident and reflected Langmuir waves, normalized to the amplitude of the electric field in homogeneous plasma $E_0$, plotted vs. a schematic illustration of the corresponding electron density profile. The parameters are $T_e = 100$ eV, $v_b/c = 0.1$, and $\psi = 5^\circ$. (a) $k_L L = 5 \times 10^2$. (b) $k_L L = 5 \times 10^3$.</figcaption></figure>

one can find from Equations (13)–(15) a well-known Airy equation for the electrostatic potential:

$$ \frac{d^2}{dz^2} \phi(\bar{z}) - \bar{z} \phi(z) = 0. \qquad (17) $$

It is convenient to present the solution by making use of Hankel functions $H_l^{(n)}$ ($n = 1, 2$). They allow one to easily separate the incident and reflected waves. Within the conversion region, the electrostatic potential is written in the following form:

$$ \phi(\bar{z}_0 < \bar{z} < 0) = \frac{1}{2} \sqrt{\frac{(-\bar{z})}{3}} \left( e^{-i\pi/6} H_{1/3}^{(2)} \left( \frac{2}{3}(-\bar{z})^{3/2} \right) + e^{i\pi/6} H_{1/3}^{(1)} \left( \frac{2}{3}(-\bar{z})^{3/2} \right) \right) = \phi_i + \phi_r. \qquad (18) $$

Here $\bar{z}_0 = -(3k_l^2 \lambda_D^2)^{2/3} (k_L L)^{2/3} \cos^2 \psi$ denotes the beginning of the density gradient, the reflection of the Langmuir wave occurs at $\bar{z} = 0$, and the first term corresponds to the incident (i) wave and the second term to the reflected (r) wave. We evaluate the amplitude of the electrostatic potential $\Phi_l^{\text{inhom}}$ (for details see Appendix B) as a function of the electrostatic potential in a homogeneous plasma:

$$ \Phi_0^{\text{inhom}} = 2\sqrt{\pi}(3k_l^2 \lambda_D^2 k_L L)^{1/6} \cos^{1/2} \psi \Phi_0^{\text{homog}}. \qquad (19) $$

The corresponding electric field components can be found as $E_l^{\text{inhom}} = -\nabla \Phi_l^{\text{inhom}}(\mathbf{r}, t)$ (Figure 2), which are expressed in

terms of $E_0$—the amplitude of the electric field in the homogeneous plasma, which can be found as $E_0 = k_l \Phi_0^{\text{homog}}$.

Having obtained the expressions for two separate components of the electric field, attributed to the incident and reflected Langmuir waves, one can derive the electron density and velocity perturbations excited in the plasma by each wave, using simple linear relations (see Appendix B). Thus, the nonlinear current resulting from the superposition of these perturbations should be written in the following form:

$$ J_{2\omega_{\text{pe}}}(\mathbf{r}) = -e(\delta n_{l_i}(\mathbf{r})\delta v_r(\mathbf{r}) + \delta n_{l_r}(\mathbf{r})\delta v_i(\mathbf{r})). \qquad (20) $$

This nonlinear current represents a localized source of generation of EM emission at about $2\omega_{\text{pe}}$ (see Figure 3).

## 4. Emission from Localized Density Perturbations

As a harmonic EM emission generated in the corona or interplanetary medium is mainly observed at distances much larger that its source size and wavelength, we will consider the EM field of this emission as such at a large distance from the source, which implies a decomposition (Landau & Lifshitz 2013):

$$ \begin{aligned} J_{2\omega_{\text{pe}}}(\mathbf{r}, t - \frac{|\mathbf{R} - \mathbf{r}|}{c}) &\approx J_{2\omega_{\text{pe}}}(\mathbf{r}, t - \frac{R}{c} + \frac{rn}{c}) \\ &= J_{2\omega_{\text{pe}}}(\bar{z})e^{-2i\omega_{\text{pe}}t + 2i\omega_{\text{pe}}R/c - 2i\omega_{\text{pe}}rn/c + 2ik_{l_x}x}. \end{aligned} \qquad (21) $$
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Absolute value of harmonic current density, normalized to the square of the amplitude of the electric field in homogeneous plasma $E_0$, inside the conversion region ($\tilde{z}_0 < \tilde{z} < 0$). The parameters are $T_e = 100$ eV, $v_b/c = 0.1$, and $\psi = 5^\circ$. (a) $k_l L = 5 \times 10^2$. (b) $k_l L = 5 \times 10^3$.</figcaption></figure>

Here $\omega_{\text{pe}} \text{rn/c}$ may be interpreted as the ratio of the emission source size to the wavelength. The Liénard–Wiechert potential of the current is

$$A_{2\omega_{\text{pe}}} (\mathbf{r}) e^{-2i\omega_{\text{pe}}t} = \frac{\sqrt{\epsilon}}{c} \int \frac{J_{2\omega_{\text{pe}}} (\mathbf{r}, t - \frac{|\mathbf{R} - \mathbf{r}|}{c})}{|\mathbf{R} - \mathbf{r}|} d^3\mathbf{r}. \quad (22)$$

It is well known that a dipolar emission is linearly proportional to the amplitude of oscillations of the center of mass of charged particles. For a system that consists only of electrons (we neglect ion motions for Langmuir waves) the center of mass may not undergo any displacement; thus dipolar emission is absent (Landau & Lifshitz 2013). The absence of a dipolar component in the mechanism of generation of EM emission by localized wave packets was first pointed out by Galeev & Krasnoselskikh (1976). So, we imply that the quadrupole component of the emission is dominant. For the very same reason, magnetic dipolar radiation is also absent. Thus we can use the decomposition

$$\frac{1}{|\mathbf{R} - \mathbf{r}|} \approx 4\pi \sum_{l=0}^{\infty} \sum_{m=-l}^{l} \frac{1}{2l+1} \frac{r^l}{R^{l+1}} Y_l^{m*}(\theta_r, \phi_r) Y_l^m(\theta_R, \phi_R), \quad (23)$$

and keep only the terms corresponding to $l=2$ that account for the quadrupolar emission. The magnetic field of the EM wave is

$$H_{t2\omega_{\text{pe}}} = 2i[k_l A_{2\omega_{\text{pe}}}], \quad (24)$$

which can be rewritten as

$$|H_{t2\omega_{\text{pe}}}| = H_y = 2i(-k_{t_x}A_{z2\omega_{\text{pe}}} + k_{t_z}A_{x2\omega_{\text{pe}}}). \quad (25)$$

We note that the term proportional to $A_{x2\omega_{\text{pe}}}$ will be the most important as the current $J_{x2\omega_{\text{pe}}}$ prevails over $J_{z2\omega_{\text{pe}}}$ for larger angles (see Figure 3), and for very small angles the product $k_{t_x}A_{z2\omega_{\text{pe}}} \sim \sin \psi A_{x2\omega_{\text{pe}}}$ vanishes. The radiant energy density of

the emission is

$$W_t^{\text{inhom}} = \frac{|H_{t2\omega_{\text{pe}}}|^2}{8\pi}. \quad (26)$$

An approximate analytical expression for $W_t^{\text{inhom}}$ from a single localized density clump is (for details see Appendix C)

$$W_t^{\text{inhom}} = 4 \cdot 10^2 \pi \epsilon \frac{v_T^6}{v_b^6} \left( \frac{k_l^2 \lambda_D^2}{\delta n / n_0} \right)^2 \frac{\omega_{\text{pe}}^2 L^2}{c^2} \sin^2 \psi \frac{W_l}{n_0 k_B T_e} W_l, \quad (27)$$

where $v_T$ is the thermal velocity of electrons in the plasma and $W_l = E_0^2/8\pi$ is the energy density of Langmuir waves. This expression is derived under the approximation $3k_l^2 \lambda_D^2 \ll \delta n / n_0$ and $\sin \psi \ll 1$.

On their way through the inhomogeneous solar wind, Langmuir waves can encounter density clumps of different sizes and magnitudes. As shown earlier, both of these parameters can strongly affect the harmonic emission from inside the density clump. For instance, Figure 2 demonstrates that with the growth of $L$ (and the decrease of the amplitude of density fluctuations $\delta n / n_0 \sim L^{-1}$) the conversion region increases, allowing the perpendicular current $J_{x2\omega_{\text{pe}}}$ to grow significantly larger than $J_{z2\omega_{\text{pe}}}$ for the major part of the angles $\psi$. To take this into account, one can estimate statistically the averaged value of $W_t^{\text{inhom}}$. The size of the source region of type III bursts is typically much larger than the characteristic scale of density fluctuations within the solar wind (Reid & Ratcliffe 2014), and thus, the number of encounters is large enough to justify averaging:

$$\langle W_t^{\text{inhom}} \rangle_{\psi, \delta n, L} = P_{\text{ref}} \int P(\psi) P(\delta n) P(L) W_t^{\text{inhom}}(\psi, \delta n, L) \times d\psi \, d\delta n \, dL, \quad (28)$$

where $P(\psi)$, $P(\delta n)$, and $P(L)$ are the probability distribution functions of the angles of incidence, the amplitudes of the density fluctuations within the clumps, and the scales of gradients inside the clumps, respectively. We suppose that the
angles $\psi$ are distributed uniformly between 0 and $\psi_{max}$, the amplitudes $\delta n$ follow a normal distribution with zero mean and standard deviation $\langle \Delta n \rangle$, and the scales follow a distribution adopted from Krasnoselskikh et al. (2019; see also Appendix D):

$$P_L(L) = \frac{1}{\sqrt{2\pi}} \frac{L_{sc}}{L^2} \exp\left[-\frac{L_{sc}^2}{2L^2}\right], \quad (29)$$

where $L_{sc}$ is the characteristic scale of density gradients for normally distributed density fluctuations. It can be approximated as a function of the level of density fluctuations:

$$L_{sc} \approx 1.4 \cdot \left( \frac{\langle \Delta n \rangle}{n_0} \right)^{-1}, \quad [\text{km}]. \quad (30)$$

As shown by Krasnoselskikh et al. (2019), random density fluctuations, described this way in terms of $P(\delta n)$ and $P(L)$, reproduce the interval of the density spectrum from $10^{-2}$ Hz to 530 Hz measured within the solar wind. The reason for choosing this part of the spectrum is that the density fluctuations that may affect the beam-plasma interaction have characteristic scales that on the one hand should be much larger than the wavelength of the Langmuir waves and, on the other hand, must be significantly smaller than the relaxation length of the beam-plasma interaction. The averaged value of the energy density of the harmonic EM emission is (for details see Appendix D)

$$\left\langle W_t^{\text{inhom}} \right\rangle_{\psi, \delta n, L} = 25\sqrt{6} \epsilon P_{\text{ref}} \frac{v_T^3}{c^3} \frac{v_T^3}{v_b^3} \left( \frac{k_l^2 \lambda_D^2}{\langle \Delta n \rangle / n_0} \right)^2 \times \frac{\omega_{P_e}^2 L_{sc}^2}{c^2} \frac{W_l}{n_0 k_B T_e} W_l. \quad (31)$$

To evaluate the efficiency of the generation mechanism, one should establish the relations between the parameters of the electron beam and the characteristics of the spectra and temporal evolution of the Langmuir waves generated via the beam-plasma interaction. We present hereafter some of the results of the probabilistic model of beam-plasma interaction in a plasma with random density fluctuations.

## 5. Electron Beam–Plasma Interaction in Randomly Inhomogeneous Plasma

Before proceeding with EM harmonic emission, we should evaluate the wave energy density of Langmuir waves. The problem of beam-plasma interaction may be analyzed by means of quasilinear theory (QLT), which takes into account the process of generation of Langmuir waves due to the bump-on-tail instability of energetic electrons, and the following modification of the electron distribution function that eventually leads to the formation of a plateau. An important condition for a QLT description of this process consists in exact resonance between the wave and the particle: in the one-dimensional case they interact only when the particle velocity is exactly equal to the phase velocity of the wave (Vedenov et al. 1962; Drummond & Pines 1964). There have been several reports exploring two-dimensional (e.g., by Ziebell et al. 2008, 2011) and three-dimensional (e.g., by Harding et al. 2020)

quasilinear wave-particle interactions, that account for the angular diffusion of the wavevectors of Langmuir waves (Nishikawa & Ryutov 1976; Krasnoselskikh et al. 2007). However, the one-dimensional approach remains justified here due to the weak angular dispersion of the beam velocities and associated Langmuir waves reported at around 1 au (Ergun et al. 1998; Malaspina & Ergun 2008; Krasnoselskikh et al. 2011). For the case of a homogeneous plasma, the QLT predicts the formation of a plateau in the electron velocity distribution function in the range of velocities from the beam velocity to the thermal velocity of the plasma. The process of “plateauing” is accompanied by a transformation of the free kinetic energy of the electrons into the potential energy of Langmuir waves. Sturrock (1964) applied the QLT description of beam-plasma interaction to solar radio bursts. The analysis of beam-plasma interaction under conditions relevant to the solar corona and solar wind resulted in the so-called “Sturrock paradox”: the relaxation of the beam should have stopped after a very short distance, about 100 km. However, Langmuir waves and associated beams were observed up to the Earth’s orbit. Later, satellite measurements showed that such beams are observed in the solar wind even at distances of about 5 au.

Recent studies (Kellogg & Horbury 2005; Krucker et al. 2009; Ratcliffe et al. 2012; Voshchepynets & Krasnoselskikh 2013) have demonstrated that there is an important characteristic of the solar wind that should be taken into account when analyzing the aforementioned processes. The solar wind is quite strongly inhomogeneous, filled with random density fluctuations $\delta n/n_0$ that may be quite intense, about several percent of the background plasma density at 1 au. Under such conditions the phase velocity $V_{ph}$ of Langmuir waves varies, since the probability distribution of the phase velocity is determined by the probability distribution of the density fluctuations due to the relation

$$\omega_{P_e}^2(n) = \omega_{P_e}^2(n_0) \left(1 + \frac{\delta n}{n_0}\right) = \omega^2(1 - 3v_T^2/V_{ph}^2), \quad (32)$$

and consequently, the wave along its path resonantly interacts with electrons of different velocities. There are several complementary models that describe wave-particle interaction in randomly inhomogeneous plasma. One is a Hamiltonian numerical model where the background plasma is described by Zakharov equations (Krafft et al. 2013; Krafft & Volokitin 2014; Volokitin & Krafft 2016, 2018, 2020), and the beam and its interaction with waves are modeled by a particle-in-cell code. In this model the system is periodic and is chosen to be long enough to incorporate several modes of density fluctuations. The second model is based on the use of the probability distribution of the wave’s phase velocity in plasma with random density fluctuations (Voshchepynets et al. 2015). The resonant wave-particle interaction takes into account the wave interaction with particles having different velocities. The probability distribution becomes a statistical weighting function that results in a natural widening of the resonance conditions.

It was shown by Voshchepynets & Krasnoselskikh (2015) that there exist two regimes of beam-plasma interaction, depending on the ratio of two important parameters of the problem: the dispersion $k_l^2 \lambda_D^2$ and the density fluctuation level $\langle \Delta n \rangle / n_0$. When the level of density fluctuations is small with respect to dispersion effects, relaxation occurs very similarly to the homogeneous case. In the opposite situation, the dynamics
<figure><img src="image_6.png" /><figcaption>**Figure 4.** (a) Temporal evolution of the energy density of Langmuir waves $W_l^{\text{temp}}$ in units of the initial energy density of the electron beam $W_{b0}$, obtained from a probabilistic model (Voshchepynets et al. 2015). Numerical solutions are found under the following conditions: electron density in the beam $n_b = 10^{-5}n_0$, $v_b = 16v_T$, and six different levels of average density fluctuations $(\Delta n)/n_0$. (b) Maximum values of Langmuir wave energy density (see panel (a)) vs. average level of density fluctuations, for three different ratios of initial beam velocity $v_b$ to electron thermal velocity $v_T$. The simulation results (solid lines) are fitted by a power law (dashed lines). (c) Time of Langmuir waves' energy growth during beam relaxation (see panel (a)) in units of growth rate in a homogeneous plasma $\gamma_{\text{homog}}$ vs. average level of density fluctuations, for three different ratios of initial beam velocity $v_b$ to electron thermal velocity $v_T$.</figcaption></figure>

of the instability is quite different. First of all, waves grow much slower (Figures 4(a), (c)) since the increment significantly decreases. Initially the waves' energy grows, reaches a maximum, and begins to decrease, as shown in Figure 4(a). The most surprising result is the transfer of a significant part of the wave energy to electrons with energies higher than the energy of the beam.

Both descriptions—probabilistic models and models based on Zakharov's equations—give very similar results (Voshchepynets et al. 2017). It was shown by Voshchepynets & Krasnoselskikh (2015) that the relaxation process consists of two stages. During the first stage, the major relaxation process occurs and its characteristic time may be determined as in conventional QLT. However, at the end of the first stage the system does not reach a stable state, only a marginally stable state, when the increment of wave growth is $10^{-4}-10^{-5}$ smaller than the initial increment. During the second stage, the system exists in this quasi-stable state and still generates waves that are significantly above noise level but have much smaller amplitudes than during the first stage of relaxation. This allows one to explain simultaneous observations of strong Langmuir waves and positive slopes on electron velocity distribution functions at large distances from the Sun (Lin et al. 1981).

Krasnoselskikh et al. (2019) recently showed that density fluctuations may also change the mechanism of generation of radio emission at a fundamental frequency close to the local plasma frequency. Basically, the mechanism consists in a direct conversion of electrostatic Langmuir waves to EM waves when the Langmuir waves are reflected from a density clump. The estimated efficiency of such a transformation may become as high as $10^{-4}$ for density fluctuations of the order of several percent.

Harmonic emission is produced via a nonlinear process, and consequently, its efficiency depends on the Langmuir wave amplitude, which can be evaluated with the help of a probabilistic model. Figure 4 shows the evolution of the energy density of Langmuir waves generated in a randomly inhomogeneous plasma by an electron beam. Results are provided for typical physical parameters in the source region of type III solar radio bursts: beam electron density $n_b = 10^{-5}n_0$, beam velocity $v_b = 16v_T$, and six levels of average density fluctuations $(\Delta n)/n_0 = 0.01, 0.02, 0.03, 0.04, 0.07$, and $0.1$. Langmuir waves' energy density is shown as a ratio of the initial energy density of the beam $W_{b0} = n_b m_e v_b^2 / 2$. It is worth mentioning that the probabilistic model describes the temporal evolution of Langmuir waves, while in the present study the problem should be treated as a spatial one, with a boundary
condition corresponding to a continuous ejection of electron beams. It is well known that the solutions of these two problems are rather similar, except that the quasi-equilibrium saturation state corresponds to a redistribution of energy fluxes rather than of the energy itself. This implies a higher level of electrostatic wave energy in the spatial problem with respect to the temporal problem. The wave energy flux moves with the waves' group velocity, and the wave energy density may be found from the solution of the temporal problem:

$$W_l^{\text{SBP}} = \frac{V_{\text{ph}}}{V_{\text{gr}}} W_l^{\text{temp}}, \qquad (33)$$

where $W_l^{\text{SBP}}$ is the wave energy in the framework of the spatial boundary problem, $W_l^{\text{temp}}$ is the wave energy corresponding to quasi-saturation in the framework of the temporal problem, and $V_{\text{gr}}$ is the group velocity of Langmuir waves. Since for the waves generated by the beam one has $V_{\text{ph}}/V_{\text{gr}} = v_b^2/3v_T^2$, for beams having velocities of 10–15$v_T$, $V_{\text{ph}}/V_{\text{gr}}$ may vary from 33 to 75. It also leads to an intensification of the waves in a relatively small region of space.

In order to estimate $E_0$, we hereafter use the maximum energy density of Langmuir waves $W_{lmax} \equiv \max(W_l)$, reached during beam relaxation in a plasma with density fluctuations (see Figure 4(a)) and obtained in the framework of the spatial boundary problem:

$$\begin{aligned} W_{lmax} &= \max(W_{lmax}^{\text{SBP}}) = \max\left(\frac{v_b^2}{3v_T^2} W_{lmax}^{\text{temp}}\right) \\ &= \frac{v_b^2}{3v_T^2} \chi(\langle\Delta n\rangle/n_0) W_{b0}, \end{aligned} \qquad (34)$$

where $\chi$ is the characteristic coefficient that shows the ratio of the maximal energy density reached during the relaxation process to the initial energy density of the beam $W_{b0}$ (Figure 4(b)) at a given level of density fluctuations.

## 6. Efficiency of Conversion of Langmuir Waves into Harmonic EM Emission

Using the results of the previous section, one can obtain $W_{lmax} = (v_b^2/v_T^2)\chi(\langle\Delta n\rangle/n_0)n_b m_e v_b^2/6$, and the linear increment of the growth of Langmuir waves may be estimated as $\gamma_{\text{lin}} = \xi(\langle\Delta n\rangle/n_0)\omega_{p_e}(n_b/n_0)(v_b^2/\Delta v_b^2)$, where $\chi(\langle\Delta n\rangle/n_0)$ is a coefficient characterizing the ratio of the maximum energy density of Langmuir waves in an inhomogeneous plasma to the initial electron beam energy and $\xi(\langle\Delta n\rangle/n_0)$ is the ratio of the increment of the instability in the inhomogeneous case to the increment in the homogeneous case. In order to evaluate it in computer simulations of beam–plasma interaction, we make a direct evaluation of the time of instability development, raising time $t_r$ for different beam velocities and levels of density fluctuations as shown in Figure 4. The values of $\chi(\langle\Delta n\rangle/n_0)/\xi(\langle\Delta n\rangle/n_0)$ typically vary from 10 to 70. The efficiency of conversion of Langmuir waves into harmonic

emission in a quasihomogeneous plasma is

$$\begin{aligned} K_{2\omega_{p_e}}^{\text{homog}} = \frac{W_{lmax}^{\text{homog}}}{W_{lmax}} &= \frac{\chi(\langle\Delta n\rangle/n_0)}{\xi(\langle\Delta n\rangle/n_0)} \Lambda \frac{m_e \Delta v_b^2}{k_B T_e} \\ &\times \frac{P_{\text{ref}}(1 - P_{\text{ref}})}{2400\pi^4} \left(\frac{v_b}{v_T}\right)^2 \left(\frac{k_b}{\Delta k_b}\right)^4 \frac{k_l^5}{k_b^5}. \end{aligned} \qquad (35)$$

Taking into account that $k_t = \sqrt{3}\omega_{p_e}/2c$, $k_b = \omega_{p_e}/v_b$, and $\Delta k_b/k_b \simeq \Delta v_b/v_b$, we obtain

$$\begin{aligned} K_{2\omega_{p_e}}^{\text{homog}} = \frac{\sqrt{3}}{38400\pi^4} \Lambda P_{\text{ref}}(1 - P_{\text{ref}}) \frac{\chi(\langle\Delta n\rangle/n_0)}{\xi(\langle\Delta n\rangle/n_0)} \\ \times \left(\frac{v_b}{c}\right)^9 \left(\frac{v_b}{\Delta v_b}\right)^2 \left(\frac{m_e c^2}{k_B T_e}\right)^2. \end{aligned} \qquad (36)$$

As a next step, we estimate the efficiency of conversion of a Langmuir wave into harmonic EM emission in the vicinity of the reflection points:

$$\begin{aligned} K_{2\omega_{p_e}}^{\text{inhom}} = \frac{\langle W_t^{\text{inhom}} \rangle_{\psi, \delta n, L}}{W_{lmax}} &= \frac{25\sqrt{6}}{3} \epsilon \chi(\langle\Delta n\rangle/n_0) \\ &\times P_{\text{ref}} \frac{n_b}{n_0} \frac{v_b}{c} \left(\frac{v_T}{c}\right)^2 \left(\frac{k_l^2 \lambda_D^2}{\langle\Delta n\rangle/n_0}\right)^2 \frac{\omega_{p_e}^2 L_{sc}^2}{c^2}. \end{aligned} \qquad (37)$$

The dependence of the efficiency coefficients $K_{2\omega_{p_e}}^{\text{homog}}$ and $K_{2\omega_{p_e}}^{\text{inhom}}$ on the plasma temperature and electron plasma frequency respectively versus the electron beam velocity is presented in Figure 5.

## 7. Radiation Directivity and Intensity in the Solar Corona and Wind

### 7.1. Radiation Directivity

The directivity of the harmonic radiation of type III radio bursts in a homogeneous plasma has been extensively discussed since the suggestion of the plasma emission mechanism by Ginzburg & Zheleznyakov (1958). According to this mechanism, the harmonic emission is always quadrupolar, and its angular range of visibility depends on the value of the beam velocity, typically larger for smaller values of beam velocity (Zheleznyakov & Zaitsev 1970). Disregarding revisions of the initial theory, the general directivity characteristics predicted for such emission remain unchanged.

In the present paper, we revisit the emission mechanism in a quasihomogeneous plasma and also consider a physical process of emission from localized regions (clumps) where reflection of Langmuir waves occurs. Similar to the plasma emission mechanism, this harmonic radiation is quadrupolar. But unlike the plasma emission mechanism, it produces EM emission in the parallel and perpendicular directions with respect to the electron beam direction (see Figure 6). In this context, we can deduce that the direction of the density gradient inside the density clump and the direction of the electron beam should be roughly similar based on the following conditions: (1) the angular range of directivity of Langmuir waves generated by the electron beam is quite narrow, up to 20° with respect to the beam direction (Krasnoselskikh et al. 2011), and (2) harmonic
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Efficiency of conversion of beam-generated Langmuir waves to harmonic EM emission in homogeneous (left panel) and inhomogeneous (right panel) plasma. The radial distance from the Sun in solar radii was inferred from temperature and ion density scaling for slow solar wind by Meyer-Vernet & Issautier (1998). Left panel: $K_{2\omega_{p_e}}^{\text{homog}}$ vs. electron beam velocity $v_b$ and electron temperature $T_e$. The parameters are $\chi/\xi = 1$, $\Lambda = 10$, $P_{\text{ref}} = 0.5$, and $\Delta v_b/v_b = 1/3$. Right panel: $K_{2\omega_{p_e}}^{\text{inhom}}$ vs. electron beam velocity $v_b$ and electron plasma frequency $f_{p_e}$. The parameters are $\epsilon = 0.75$, $P_{\text{ref}} = 0.5$, $n_b = 2.5 \times 10^{-5} n_0$, $T_e = 100$ eV, and $(\Delta n)/n_0 = 0.1$ (corresponding to $\chi \approx 0.1$, $L_{sc} \approx 1.4 \times 10^6$ cm).</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 6.** A schematic radiation pattern of harmonic EM emission produced via the plasma emission mechanism in homogeneous plasma (multicolored) and via Langmuir wave coalescence inside density clumps (orange). Here the ratio between homogeneous and inhomogeneous parts is chosen to be 5:1.</figcaption></figure>

emission in our model can be generated only for an incident Langmuir wave that is highly aligned with the density gradient direction, since the angle of incidence $\psi$, which allows the successful production of harmonics, is strictly determined by the beam velocity, $\psi \approx \sqrt{3} v_b/2c$, and is small ($\lesssim 10^\circ$).

We can schematically compare the two aforementioned radiation patterns by positioning them against the direction of the electron beam (Figure 6). Here we do not precisely represent the magnitudes of the relevant intensity or angular range of radiation but rather make a formal comparison of the

major emission directions given by these two mechanisms combined. Disregarding the fact that propagation effects, such as refraction, absorption, and reflection (of backward-moving emission), cause a widening of the angular range where harmonic emission from homogeneous plasma is visible, they are insufficient for explaining a widespread visibility of harmonic emission of type III radio bursts (Thejappa et al. 2007). As such, the mechanism of generation of radio emission inside density clumps adds parallel- and perpendicular-directed radiation to the conventional harmonic radiation pattern and contributes to the general visibility of harmonic emission of type IIIs.

## 7.2. Radiation Intensity in the Corona and Solar Wind

Observations of type III radio bursts indicate that fundamental–harmonic pairs represent the majority of radio bursts in the high-frequency range. The fundamental usually begins below 100 MHz, while the harmonic can begin at as high as ~500 MHz (Dulk & Suzuki 1980). At the same time, the rarity of fundamental emission in the 100–500 MHz range remains unexplained, as absorption due to inverse bremsstrahlung becomes significant only above around 500 MHz for fundamental emission and above 1 GHz for harmonic emission (Reid & Ratcliffe 2014). This might indicate that in the corona and its vicinity the efficiency of conversion of Langmuir waves into EM radio emission is higher for certain mechanisms of harmonic radiation generation. Below, we formally discuss the intensity of the radiation, recalling its simple relation to the wave energy density $I \sim W_t$.

We have revisited a well-known result for harmonic emission from a quasihomogeneous plasma and obtained an analytical result that is different from the one obtained by Willes et al. (1996). We have shown that the assumptions used by Willes et al. (1996) are not always justified and they lead to a sufficient underestimation of the EM wave amplitude. We have performed here direct calculations for the general case. According to our results, the intensity of such emissions is much higher than previously predicted. Such emissions are
more efficiently produced for higher ratios of $v_b$ to $c$ and lower electron temperatures $T_e$ (see the left panel of Figure 5). As electron temperature decreases with heliocentric distance, we can infer a dependence of $T_e$ on radial distance from the Sun, making use of one of the solar wind models (e.g., Meyer-Vernet & Issautier 1998), and make a comparison between the parameter domains of the domination of harmonic emissions from a homogeneous plasma and from density clumps.

Harmonic emission from density clumps is most intense at smaller heliocentric distances. Its efficiency of generation is higher for lower ratios of $v_b$ to $c$ and for smaller levels of density fluctuations (see the right panel of Figure 5). Here we note that we have applied a constant level of electron temperature $T_e$ and density fluctuations throughout the whole plasma frequency interval. Recent studies show that, according to in situ measurements of PSP, the level of density fluctuations at around $36 R_\odot$ is about 0.06–0.07, and it is predicted to grow up to ~0.2 at a distance of a few solar radii (Krupar et al. 2020). On the right panel of Figure 5 we assume a level of density fluctuations $\langle\Delta n\rangle/n_0 = 0.1$.

As we compare the efficiency of conversion of Langmuir waves into harmonic emission and thus indirectly the relevant intensities, we see that emission from density clumps can become as important as emission from quasihomogeneous plasma at around a few solar radii, closer to the low corona.

### 7.3. Summary

It is widely accepted that electron density fluctuations in the solar wind affect the propagation of radio emission. On the other hand, there are very few studies of the impact of these inhomogeneities on the process of generation of such radio emissions. In the present paper, we consider the generation of harmonic radio emission via the $l + l' \rightarrow t$ process under two different circumstances: in a quasihomogeneous plasma (through the coalescence of two nearly oppositely propagating Langmuir waves) and inside structures formed by density fluctuations with increasing density gradients (via the coalescence of a Langmuir wave with its reflected part in the vicinity of the reflection point). For the first generation process we make the following assumptions: (1) coalescence takes place in a homogeneous plasma, (2) the spectrum of forward-moving and reflected Langmuir waves is Gaussian, (3) the population of reflected waves is the result of the reflection of a part of forward-moving waves from density irregularities, and the reflection process is taken into account by means of a

coefficient $P_{\text{ref}}$, and (4) the two coalescing Langmuir waves meet head-on. For the second generation process, we use the system of equations proposed by Zakharov (1972). The assumptions we make when deriving the solution are as follows: (1) $3k_l^2 \lambda_D^2 \ll \langle\Delta n\rangle/n_0$, i.e., linear dispersion is less significant compared to the effect of density fluctuations, (2) $k_l L_{sc} \gg 1$, i.e., the characteristic scales of density gradients inside density clumps are significantly larger than the wavelength of the Langmuir wave, (3) $\psi \ll 1$, i.e., the incident Langmuir wave should be closely aligned with the direction of the density gradient and, consequently, almost antiparallel to the reflected wave (HOA), (4) clumps are approximately spherical with the electron density increasing linearly toward their center, (5) the radio emission is formed within the conversion region, and (6) the quadrupole component of the harmonic emission is dominant. In both cases, the plasma is unmagnetized and we consider only the conversion of Langmuir waves into harmonic emission. This way we obtain analytical expressions for the energy density of harmonic radio emission in both cases. And finally, we estimate the efficiency of conversion of beam-generated Langmuir waves into harmonic EM emission from both regions of emission.

## 8. Conclusions

(1) A direct calculation of the generation of harmonic EM emission via the process of coupling of primary beam-generated Langmuir waves with reflected waves in a quasihomogeneous plasma yields a higher radiation intensity than found previously (e.g., by Willes et al. 1996).

(2) The new model of generation of harmonic emission inside density clumps close to the region of reflection of Langmuir waves demonstrates the efficiency of conversion of Langmuir waves into EM waves, which is under certain conditions comparable with the aforementioned quasihomogeneous plasma emission and even prevails at smaller heliocentric distances.

(3) EM radiation from density clumps may be important for the visibility of harmonic emissions.

V.K. acknowledges financial support from CNES through grants “Search Coil for Solar Orbiter” and “Parker Solar Probe” and financial support from NASA through grant 80NSSC20K0697. The authors are grateful to Didier Mourenas for useful discussions.
## Appendix A
### Nonlinear Coupling of Langmuir Waves in a Quasihomogeneous Plasma

The process of EM wave generation by the coupling of Langmuir waves is described by the following set of equations (Tsytovich 2012):

$$ \frac{dN_t(\mathbf{k}_t)}{dt} = \iint \frac{d^3\mathbf{k}_1 d^3\mathbf{k}_2}{(2\pi)^6} w_{ll}^t(\mathbf{k}_1, \mathbf{k}_2, \mathbf{k}_t) [N_l(\mathbf{k}_1)N_l(\mathbf{k}_2) - N_l(\mathbf{k}_1)N_l(\mathbf{k}_t) - N_l(\mathbf{k}_2)N_l(\mathbf{k}_t)], \quad (A1) $$

$$ w_{ll}^t(\mathbf{k}_1, \mathbf{k}_2, \mathbf{k}_t) = \frac{(2\pi)^6 e^2 (k_1^2 - k_2^2)^2 [\mathbf{k}_1 \times \mathbf{k}_2]^2}{32\pi m_e^2 \omega_{P_e} k_t^2} \delta(\omega_t - \omega_1 - \omega_2) \delta(\mathbf{k}_t - \mathbf{k}_1 - \mathbf{k}_2). \quad (A2) $$

For our problem the number of quanta of Langmuir waves is supposed to be much larger than the number of quanta of EM waves; thus Equation (A1) may be simplified as follows:

$$ \frac{dN_t(\mathbf{k}_t)}{dt} = \iint \frac{d^3\mathbf{k}_1 d^3\mathbf{k}_2}{(2\pi)^6} w_{ll}^t(\mathbf{k}_1, \mathbf{k}_2, \mathbf{k}_t) N_l(\mathbf{k}_1) N_l(\mathbf{k}_2). \quad (A3) $$

Aiming to evaluate the number of quanta of EM waves we choose a distribution that consists of two symmetric parts, primarily generated Langmuir waves and reflected ones that have similar distribution functions but different intensities. Like Willes et al. (1996), we set both distributions to be Gaussian and centered at $\mathbf{k}_b$ ($\mathbf{k}_b = \omega_{P_e}/v_b$) for forward-moving waves and at $-\mathbf{k}_b$ for reflected waves. This suggests that the reflections are of the mirror type. Under such conditions the total distribution can be written as follows:

$$ N_l(\mathbf{k}_\parallel, \mathbf{k}_\perp) = (1 - P_{\text{ref}})N \exp\left(-\frac{(\mathbf{k}_\parallel - \mathbf{k}_b)^2 + \mathbf{k}_\perp^2}{\Delta k^2}\right) + P_{\text{ref}}N \exp\left(-\frac{(\mathbf{k}_\parallel + \mathbf{k}_b)^2 + \mathbf{k}_\perp^2}{\Delta k^2}\right), \quad (A4) $$

where N is normalized according to

$$ N = \frac{1}{\pi^{3/2} \Delta k^3} \frac{W_l}{\omega_{P_e}}. \quad (A5) $$

This way Equation (A3) takes the form

$$ \frac{dN_t(\mathbf{k}_t)}{dt} = \frac{e^2 P_{\text{ref}} (1 - P_{\text{ref}}) N^2}{32 \pi m_e^2 \omega_{P_e} (\pi^3 \Delta k^6)} \times \iint d^3\mathbf{k}_1 d^3\mathbf{k}_2 \frac{(k_1^2 - k_2^2)^2 [\mathbf{k}_1 \times \mathbf{k}_2]^2}{k_t^2 k_2^2} \\ \times \delta(\omega_t - \omega_1 - \omega_2) \delta(\mathbf{k}_t - \mathbf{k}_1 - \mathbf{k}_2) \exp\left[-\frac{(k_1 \parallel - k_b)^2 + k_t^2}{\Delta k^2}\right] \exp\left[-\frac{(k_2 \parallel + k_b)^2 + k_2^2}{\Delta k^2}\right]. \quad (A6) $$

We use a reference frame where the z-axis is in a parallel direction along the direction of propagation of a beam that is generating primary Langmuir waves (or along the magnetic field, which is not taken into account here but is still present in the solar wind). The x-axis is a perpendicular direction that is chosen to be along the second component of the $\mathbf{k}_t$ vector. The three vectors $\mathbf{k}_t, \mathbf{k}_1$, and $\mathbf{k}_2$ make a triangle (see Figure 7) according to momentum conservation so that

$$ \mathbf{k}_2 = \mathbf{k}_t - \mathbf{k}_1. \quad (A7) $$

The natural assumption here is $|\mathbf{k}_t| \ll |\mathbf{k}_{1,2}|$; thus

$$ \mathbf{k}_2^2 = \mathbf{k}_1^2 - 2k_t \parallel k_1 - 2k_1 \perp k_t \perp \cos \varphi \approx \mathbf{k}_1^2 - 2k_t k_1 (\cos \psi \cos \theta + \sin \psi \sin \theta \cos \varphi), \quad (A8) $$

where $\psi$ and $\theta$ are the angles between the z-axis and vectors $\mathbf{k}_1$ and $\mathbf{k}_t$, respectively, $\varphi$ is the angle between the projections of $\mathbf{k}_1$ and $\mathbf{k}_t$ to the plane perpendicular to the z-axis (xy-plane). The vector $\mathbf{k}_t$ without loss of generality is chosen to be in the xz-plane. Consequently the multiplier from Equation (A6) may be rewritten as

$$ \frac{(k_1^2 - k_2^2)^2 [\mathbf{k}_1 \times \mathbf{k}_2]^2}{k_t^2} = 4k_t^2 \times R(\psi, \theta, \varphi), \quad (A9) $$

where the trigonometric expression $R(\psi, \theta, \varphi)$ stands for

$$ R(\psi, \theta, \varphi) = (\sin^2 \psi \sin^2 \theta \cos^2 \varphi + 2 \sin \psi \cos \psi \\ \times \sin \theta \cos \theta \cos \varphi + \cos^2 \psi \cos^2 \theta) \times (-\sin^2 \theta \sin^2 \psi \cos^2 \varphi - 2 \sin \theta \\ \times \cos \theta \sin \psi \cos \psi \cos \varphi + (\cos^2 \psi \sin^2 \theta + \sin^2 \psi)). \quad (A10) $$

After the integration of Equation (A6) over $\mathbf{k}_2$ we obtain

$$ \frac{dN_t(\mathbf{k}_t)}{dt} = \frac{e^2 P_{\text{ref}} (1 - P_{\text{ref}}) N^2 \mathbf{k}_t^2}{8 \pi m_e^2 \omega_{P_e} (\pi^3 \Delta k^6)} \times G, \quad (A11) $$
<figure><img src="image_9.png" /><figcaption>Figure 7. Schematic illustration of wavevectors $k_1$ and $k_t$ against the direction of the electron beam (z-axis).</figcaption></figure>

where we use the notation

$$G = \int_{-\pi}^{\pi} d\varphi \int_{0}^{\pi} \sin \psi d\psi \int k^2 dk \delta(\omega_t - \omega_1 - \omega_2) R(\psi, \theta, \varphi) \\ \times \exp\left(-\frac{(k_1 \parallel - k_b)^2 + k_1^2}{\Delta k^2}\right) \exp\left(-\frac{(k_r \parallel - k_1 + k_b)^2 + (k_1 \perp + k_r^2 - 2k_1 k_r \cos \varphi)^2}{\Delta k^2}\right). \quad (A12)$$

The first integration over $k_1$ should take into account the delta function over frequencies. In the first-order approximation, neglecting terms that are linear on $k_t$, the argument of the delta function may be rewritten as follows:

$$\omega_t - \omega_1 - \omega_2 = \omega_t - 2\omega_{pe} - 3\omega_{pe}\lambda_D^2 k_1^2 = -3\omega_{pe}\lambda_D^2(k_1^2 - k_r^2), \quad (A13)$$

where

$$k_r^2 = \frac{\omega_t - 2\omega_{pe}}{3\omega_{pe}\lambda_D^2}. \quad (A14)$$

Thus the delta function has the form

$$\delta(\omega_t - \omega_1 - \omega_2) = \frac{1}{6\omega_{pe}\lambda_D^2 k_r} \delta(k_1 - k_r), \quad (A15)$$

and consequently

$$G = \frac{k_r}{6\omega_{pe}\lambda_D^2} \int_{-\pi}^{\pi} d\varphi \int_{0}^{\pi} \sin \psi d\psi R(\psi, \theta, \varphi) \exp\left(-\frac{2(k_r^2 - 2k_r k_b \cos \psi + k_b^2)}{\Delta k^2}\right). \quad (A16)$$

The expression $R(\psi, \theta, \varphi)$ may be rewritten as follows:

$$\begin{aligned} R(\theta, \psi, \varphi) = & -\frac{1}{8} \cos 4\varphi \sin^4 \theta \sin^4 \psi - \cos 3\varphi \sin^3 \theta \cos \theta \sin^3 \psi \cos \psi \\ & + \frac{1}{2} \cos 2\varphi [-\sin^4 \theta \sin^4 \psi + \sin^4 \theta \sin^2 \psi \cos^2 \psi + \sin^4 \psi \sin^2 \theta \\ & - 5 \sin^2 \theta \cos^2 \theta \sin^2 \psi \cos^2 \psi] + \cos \varphi [-3 \sin^3 \theta \cos \theta \sin^3 \psi \cos \psi \\ & - 2 \sin \theta \cos^3 \theta \sin \psi \cos^3 \psi + 2 \sin^3 \theta \cos \theta \sin \psi \cos^3 \psi + 2 \sin \theta \cos \theta \sin^3 \psi \cos \psi] + \\ & - \frac{3}{8} \sin^4 \theta \sin^4 \psi + \frac{1}{2} [\sin^4 \theta \sin^2 \psi \cos^2 \psi + \sin^4 \psi \sin^2 \theta \\ & - 5 \sin^2 \theta \cos^2 \theta \sin^2 \psi \cos^2 \psi] + (\cos^2 \psi \sin^2 \theta + \sin^2 \psi) \cos^2 \psi \cos^2 \theta. \end{aligned} \quad (A17)$$
The next step consists in integration over $d\varphi$, taking into account the following relation:

$$ \int_{0}^{2\pi} d\varphi \cos n\varphi \exp(Z \cos \varphi) = I_n(Z), \qquad (A18) $$

where $I_n(Z)$ is a modified Bessel function of the nth order. We should point out a very important difference of our calculation with the one by Willes et al. (1996). The assumption made by Willes et al. (1996) consists in the inequality

$$ \frac{2k_r k_t \sin \psi \sin \theta}{\Delta k^2} \gg 1, \qquad (A19) $$

and as will be seen later it is not satisfied for this calculation. Indeed, the inequality

$$ \frac{2k_r k_t}{\Delta k^2} \gg 1 \qquad (A20) $$

may be satisfied since $\Delta k \ll k_r$ (the spectrum of generated waves may be considered to be rather narrow), but the multiplier $\sin \psi \sin \theta$ may be quite small as we use a head-on approximation that naturally comes from our problem statement. Integration over $d\varphi$ results in

$$ G = \frac{k_r}{6\omega_{\text{pe}}\lambda_D^2} \int_0^\pi \sin \psi d\psi \exp\left(-\frac{2(k_r^2 - 2k_r k_b \cos \psi + k_b^2)}{\Delta k^2}\right) \\ \times \exp\left(-\frac{2k_r \cos \theta(k_b - k_r \cos \psi)}{\Delta k^2}\right) \times \left\{ -\frac{1}{8} \sin^4 \theta \sin^4 \psi I_4(\Xi) \right. \\ - \sin^3 \theta \cos \theta \sin^3 \psi \cos \psi I_3(\Xi) + \frac{1}{2}[-\sin^4 \theta \sin^4 \psi + \sin^4 \theta \sin^2 \psi \cos^2 \psi + \sin^4 \theta \sin^2 \theta \\ - 5 \sin^2 \theta \cos^2 \theta \sin^2 \psi \cos^2 \psi] I_2(\Xi) + [-3 \sin^3 \theta \cos \theta \sin^3 \psi \cos \psi - 2 \sin \theta \cos^3 \theta \sin \psi \cos^3 \psi \\ + 2 \sin^3 \theta \cos \theta \sin \psi \cos^3 \psi + 2 \sin \theta \cos \theta \sin^3 \psi \cos \psi] I_1(\Xi) \\ + [-\frac{3}{8} \sin^4 \theta \sin^4 \psi + \frac{1}{2}[\sin^4 \theta \sin^2 \psi \cos^2 \psi + \sin^4 \theta \sin^2 \theta \\ - 5 \sin^2 \theta \cos^2 \theta \sin^2 \psi \cos^2 \psi] \\ + (\cos^2 \psi \sin^2 \theta + \sin^2 \psi)\cos^2 \psi \cos^2 \theta] I_0(\Xi)\}, \qquad (A21) $$

where $\Xi = \frac{2k_r k_t \sin \psi \sin \theta}{\Delta k^2}$. The last step in our calculation consists in the calculation of the following integral:

$$ \int_{-1}^{1} d(\cos \psi) \exp\left(\frac{4k_r k_b \cos \psi}{\Delta k^2}\right) g(\cos \psi), \qquad (A22) $$

where the factor $4k_r k_b / \Delta k^2$ is very large. The standard asymptotic estimation of such an integral when the function has the maximum on the given interval and when the parameter $Y \gg 1$ is, according to Wasow (2018),

$$ \int_{-1}^{1} dt \exp(Yt) g(t) dt = \frac{1}{Y} \exp(Y) g(1) + O(Y^{-2}). \qquad (A23) $$

Accordingly Equation (A21) reduces to

$$ G = \frac{\Delta k^2}{6\omega_{\text{pe}}\lambda_D^2 k_b} \sin^2 \theta \cos^2 \theta \exp\left[-\frac{2(k_r - k_b)^2}{\Delta k^2} - \frac{2k_r \cos \theta(k_b - k_r)}{\Delta k^2}\right]. \qquad (A24) $$

Taking into account that the maximum of this expression corresponds to $k_r = k_b$, we estimate the result as

$$ \frac{dN_t(k_t)}{dt} = \frac{P_{\text{ref}}(1 - P_{\text{ref}})\omega_{\text{pe}}^2}{768\pi^5} \left(\frac{k_b}{\Delta k}\right)^4 \frac{k_t^2}{k_b^5} \frac{N_t^2}{n_0 m_e v_t^2} \sin^2 \theta \cos^2 \theta, \qquad (A25) $$

or

$$ \frac{dW_t^{\text{homog}}(k_t)}{dt} = \frac{P_{\text{ref}}(1 - P_{\text{ref}})\omega_{\text{pe}}}{768\pi^5} \left(\frac{k_b}{\Delta k}\right)^4 \frac{k_t^2}{k_b^5} \frac{W_t}{n_0 k_B T_e} W_t \sin^2 \theta \cos^2 \theta. \qquad (A26) $$

In order to evaluate the characteristic energy of the EM wave, we integrate over $d^3k_t$

$$ \frac{dW_t^{\text{homog}}}{dt} = \frac{P_{\text{ref}}(1 - P_{\text{ref}})\omega_{\text{pe}}}{7200\pi^4} \left(\frac{k_b}{\Delta k}\right)^4 \frac{k_t^5}{k_b^5} \frac{W_t}{n_0 k_B T_e} W_t. \qquad (A27) $$

As the temporal evolution of the Langmuir wave energy density may be presented in the form

$$ W_t(t) = W_{\text{noise}} \exp(\gamma t), \qquad (A28) $$
the total energy density of EM harmonic emission will be

$$W_l^{\text{homog}} = \frac{P_{\text{ref}}(1 - P_{\text{ref}})\omega_{\text{pe}}}{7200\pi^4\gamma} \left(\frac{k_b}{\Delta k}\right)^4 \frac{\mathbf{k}_l^5}{k_b^5} \frac{W_l}{n_0 k_B T_e} W_l. \quad (A29)$$

A conventional estimate of the time of growth, namely $\frac{1}{\gamma}$, will be

$$\frac{1}{\gamma} = \frac{\Lambda}{\gamma_{\text{lin}}}, \quad (A30)$$

where $\gamma_{\text{lin}}$ is the linear increment of the instability of Langmuir waves.

## Appendix B
Equations for Electrostatic Potential, Electric Field, and Current Density

We start with the equation system for plasma oscillations with $k_l \lambda_D \ll 1$ derived by Zakharov (1972):

$$\Delta\Phi_l^{\text{inhom}}(\mathbf{r}, t) = 4\pi e \delta n_l, \quad (B1)$$

$$\frac{\partial}{\partial t} \delta n_l + \nabla((n_0 + \delta n)\mathbf{v}_e) = 0, \quad (B2)$$

$$\frac{\partial \mathbf{v}_e}{\partial t} = \frac{e}{m_e} \nabla \Phi_l^{\text{inhom}}(\mathbf{r}, t) - 3v_T^2 \nabla \frac{\delta n_l}{n_0}, \quad (B3)$$

where electron density has the following form:

$$n_e = n_0 + \delta n + \delta n_l.$$

After applying a method of small perturbations to our system and performing simple transformations, we obtain

$$\frac{\partial^2}{\partial t^2} \frac{\Delta \Phi_l^{\text{inhom}}(\mathbf{r}, t)}{4\pi e} + n_0 \cdot \left( \frac{e}{m} \Delta \Phi_l^{\text{inhom}}(\mathbf{r}, t) - \frac{3v_T^2}{4\pi e n_0} \Delta(\Delta \Phi_l^{\text{inhom}}(\mathbf{r}, t)) \right) + \nabla \left( \delta n \left( \frac{e}{m_e} \nabla \Phi_l^{\text{inhom}}(\mathbf{r}, t) - \frac{3v_T^2}{4\pi e n_0} \nabla(\Delta \Phi_l^{\text{inhom}}(\mathbf{r}, t)) \right) \right) = 0. \quad (B4)$$

Bearing in mind that each quantity in Equation (B4) except $\Phi_l^{\text{inhom}}(\mathbf{r}, t)$ is assumed to be independent of the spatial coordinates or to be a very slowly varying function of those coordinates relative to $\Phi_l^{\text{inhom}}(\mathbf{r}, t)$, we rewrite Equation (B4) as

$$\Delta \left( -\frac{\omega_l^2}{4\pi e} \Phi_l^{\text{inhom}}(\mathbf{r}, t) + \frac{e n_0}{m} \Phi_l^{\text{inhom}}(\mathbf{r}, t) - \frac{3v_T^2}{4\pi e} \Delta \Phi_l^{\text{inhom}}(\mathbf{r}, t) + \frac{e \delta n}{m} \Phi_l^{\text{inhom}}(\mathbf{r}, t) - \frac{3v_T^2 \delta n}{4\pi e n_0} \Delta \Phi_l^{\text{inhom}}(\mathbf{r}, t) \right) = 0$$

or, after a few transformations,

$$\left( \omega_l^2 - \omega_{\text{pe}}^2 \left( 1 + \frac{\delta n}{n_0} \right) + 3v_T^2 \Delta \right) \Phi_l^{\text{inhom}}(\mathbf{r}, t) = 0. \quad (B5)$$

The density gradient profile within the clump is assumed to be linear, directed along the z-axis, and may be expressed by means of the Heaviside step function $\theta(z)$:

$$\frac{\delta n}{n_0} = \theta(z) \frac{z}{L}.$$

Now, since there is no change of parameters along the x-axis and since we choose a solution with the form $\Phi_l^{\text{inhom}}(\mathbf{r}, t) = \Phi_0^{\text{inhom}}\phi(z)\exp(-i\omega_{\text{pe}}t + ik_{l,x}x)$,

$$\Delta \rightarrow \frac{d^2}{dz^2} - k_{l,x}^2,$$

and taking into account that in homogeneous plasma

$$\omega_l^2 = \omega_{\text{pe}}^2 (1 + 3k_{l,x}^2 \lambda_D^2 + 3k_{l,0}^2 \lambda_D^2),$$

where $k_{l,0} = k_l \cos \psi$, we obtain

$$\left( 3k_l^2 \lambda_D^2 \cos^2 \psi - \theta(z) \frac{z}{L} + 3\lambda_D^2 \frac{d^2}{dz^2} \right) \phi(z) = 0. \quad (B6)$$
Then, considering only the region of positive values of z hereafter, which denotes the region with increasing density, we have

$$ \left( \frac{d^2}{dz^2} + k_l^2 \cos^2 \psi - \frac{1}{3 \lambda_D^2 L} \right) \phi(z) = 0. \qquad (B7) $$

We introduce a new dimensionless variable,

$$ \bar{z} = \left( \frac{k_l^2 L^2}{3 k_l^2 \lambda_D^2} \right)^{1/3} \left( \frac{z}{L} - 3 k_l^2 \lambda_D^2 \cos^2 \psi \right), \qquad (B8) $$

or, if we introduce a new parameter $\alpha = 3k_l^2 \lambda_D^2 k_l L$, we may rewrite the previous expression:

$$ \bar{z} = \alpha^{-1/3} k_l z - \alpha^{2/3} \cos^2 \psi. \qquad (B9) $$

Equation (B7) reduces to the Airy equation:

$$ \frac{d^2}{d\bar{z}^2} \phi(\bar{z}) - \bar{z} \phi(\bar{z}) = 0. \qquad (B10) $$

It is convenient to write down the solution in the form of Hankel functions $H_l^{(n)}$ ($n=1, 2$) in order to easily separate the incident and reflected waves. Within the density clump we have two regions: the conversion region, $\bar{z}_0 < \bar{z} < 0$, where $\bar{z}_0 = -\alpha^{2/3} \cos^2 \psi$ ($z=0$) corresponds to the point where density starts to increase, and the region behind the reflection point, $\bar{z} > 0$. Inside the conversion region the solution is

$$ \phi(\bar{z}_0 < \bar{z} < 0) = \frac{1}{2} \sqrt{\frac{(-\bar{z})}{3}} \left( e^{-i\pi/6} H_{1/3}^{(2)}\left(\frac{2}{3}(-\bar{z})^{3/2}\right) + e^{i\pi/6} H_{1/3}^{(1)}\left(\frac{2}{3}(-\bar{z})^{3/2}\right) \right) = \phi_i + \phi_r, \qquad (B11) $$

In our case the Hankel function of the second kind corresponds to the incident (i) wave and that of the first kind to the reflected (r) wave. After the conversion region, the wave simply damps according to the solution

$$ \phi(\bar{z} > 0) = \frac{1}{\pi} \sqrt{\frac{\bar{z}}{3}} K_{1/3}\left(\frac{2}{3}\bar{z}^{3/2}\right). \qquad (B12) $$

In order to evaluate the amplitude $\Phi_0^{\text{inhom}}$ of the electrostatic potential, we use a Wentzel-Kramers-Brillouin (WKB) approximation to solve Equation (B10). First we assume that the solution of Equation (B10) has the form

$$ \phi(\bar{z})_{\text{WKB}} = \phi_0(\bar{z}) e^{-ik_l \Psi(\bar{z})}, \qquad (B13) $$

where the amplitude $\phi_0(\bar{z})$ and phase $\Psi(\bar{z})$ vary slowly with $\bar{z}$. We substitute this solution into Equation (B10) and obtain ($\phi$ denotes $d/d\bar{z}$)

$$ \phi_0'' - 2ik_l \Psi' \phi_0' - ik_l \Psi'' \phi_0 - k_l^2 (\Psi')^2 \phi_0 - \bar{z} \phi_0 = 0. \qquad (B14) $$

Afterward, we divide the whole equation by $k_l^2 = 4\pi^2/\lambda_l^2$ and note that according to our assumption, $\phi_0$ and $\Psi$ change noticeably only on scales $l \gg \lambda_l$. For this reason we can make an estimation of $\phi_0'' \sim \phi_0/l^2$, $\phi_0' \sim \phi_0/l$, and $\Psi'' \sim \Psi'/l$. Equation (B14) takes the form

$$ \frac{\lambda_l^2}{4\pi^2} \frac{\phi_0}{l^2} - 2i \frac{\lambda_l}{2\pi} \frac{\phi_0}{l} \Psi' - i \frac{\lambda_l}{2\pi} \frac{\Psi'}{l} \phi_0 - (\Psi')^2 \phi_0 - \frac{\lambda_l^2}{4\pi^2} \bar{z} \phi_0 = 0. \qquad (B15) $$

We may find an approximate solution by assigning terms of different orders of $\lambda_l/l$ equal to zero:

$$ ((\Psi')^2 + \frac{\lambda_l^2}{4\pi^2} \bar{z}) \phi_0 = 0, \qquad (B16) $$

$$ \phi_0' + \frac{\Psi''}{2\Psi'} \phi_0 = 0. \qquad (B17) $$

From Equation (B16) we obtain ($\phi_0 \neq 0$)

$$ \Psi' = \pm \sqrt{-\frac{\lambda_l^2}{4\pi^2} \bar{z}}, \qquad (B18) $$

$$ \Psi = \pm \sqrt{\frac{\lambda_l^2}{4\pi^2}} \int_{\bar{z}_0}^{\bar{z}} \sqrt{-\bar{z}} d\bar{z} = \mp \frac{\lambda_l}{2\pi} \left( \frac{2}{3}(-\bar{z}_0)^{3/2} - \frac{2}{3}(-\bar{z})^{3/2} \right). \qquad (B19) $$

Afterward we proceed to the solution of Equation (B17):

$$ \phi_0 = \frac{C}{\sqrt{\Psi'}}, \qquad (B20) $$
where C is an integration constant. Coming back to solution (B13) and substituting (B18) into (B20) we obtain

$$ \phi(\bar{z})_{\text{WKB}} = \frac{C\sqrt{k_l}}{(-\bar{z})^{1/4}} e^{\pm i\left(\frac{2}{3}(-\bar{z}_0)^{3/2} - \frac{2}{3}(-\bar{z})^{3/2}\right)}. \qquad (B21) $$

Sewing is performed at the point $\bar{z} = \bar{z}_0$ ($z=0$); thus we rewrite the previous expression:

$$ \phi(\bar{z} = \bar{z}_0)_{\text{WKB}} = \frac{C\sqrt{k_l}}{(-\bar{z}_0)^{1/4}}. \qquad (B22) $$

The complete solution for the electrostatic potential in the WKB approximation is

$$ \Phi_l^{\text{inhom}}(\bar{z})_{\text{WKB}} = \Phi_0^{\text{inhom}} \frac{C\sqrt{k_l}}{(-\bar{z}_0)^{1/4}} e^{\pm i\left(\frac{2}{3}(-\bar{z}_0)^{3/2} - \frac{2}{3}(-\bar{z})^{3/2}\right) - i\omega_{\text{pe}}t + i\zeta + ik_l \sin \psi x}, \qquad (B23) $$

where $\zeta$ is the phase difference between the exact solution and the solution under the WKB approximation for $\Phi_l^{\text{inhom}}(\bar{z})$. Now we may proceed to sewing the incident waves from homogeneous plasma with our approximate WKB solution (we omit $e^{-i\omega_{\text{pe}}t+ik_l \sin \psi x}$ terms, common for both waves):

$$ \Phi_{\text{incid}}^{\text{homog}}(\bar{z}_0) = \Phi_0^{\text{homog}} e^{i\eta}, \qquad (B24) $$

$$ \Phi_{\text{incid}}^{\text{inhom}}(\bar{z}_0)_{\text{WKB}} = \Phi_0^{\text{inhom}} \frac{C\sqrt{k_l}}{(-\bar{z}_0)^{1/4}} e^{i\zeta}, \qquad (B25) $$

and putting the equality

$$ \Phi_0^{\text{homog}} e^{i\eta} = \Phi_0^{\text{inhom}} \frac{C\sqrt{k_l}}{(-\bar{z}_0)^{1/4}} e^{i\zeta}, \qquad (B26) $$

we obtain $\zeta = \eta$ and

$$ C = \frac{\Phi_0^{\text{homog}}}{\Phi_0^{\text{inhom}}} \frac{(-\bar{z}_0)^{1/4}}{\sqrt{k_l}}. \qquad (B27) $$

Thus a solution for the incident wave in the WKB approximation has the form

$$ \Phi_{\text{incid}}^{\text{inhom}}(\bar{z})_{\text{WKB}} = \Phi_0^{\text{homog}} \frac{(-\bar{z}_0)^{1/4}}{(-\bar{z})^{1/4}} e^{-i\left(\frac{2}{3}(-\bar{z}_0)^{3/2} - \frac{2}{3}(-\bar{z})^{3/2}\right) - i\omega t + ik_l \sin \psi x + i\eta}. \qquad (B28) $$

Now we need to figure out the phase $\eta$. In order to do this, we use an asymptotic expansion of the exact solution for the incident wave (see Equation (B11)) for large values of the argument:

$$ \Phi_{\text{incid}}^{\text{inhom}}(\bar{z}) = \Phi_0^{\text{inhom}} \frac{(-\bar{z})^{-1/4}}{2\sqrt{\pi}} e^{-i\left(\frac{2}{3}(-\bar{z})^{3/2} - \frac{1}{4}\pi\right)}. \qquad (B29) $$

We set the expressions for $\Phi_{\text{incid}}^{\text{inhom}}(\bar{z})$ and $\Phi_{\text{incid}}^{\text{inhom}}(\bar{z})_{\text{WKB}}$ equal at the point $\bar{z} = \bar{z}_0 = -\alpha^{2/3} \cos^2 \psi$ and obtain

$$ \Phi_0^{\text{homog}} \cos^{1/2} \psi e^{i\frac{2}{3}\alpha \cos^3 \psi + i\eta} = \Phi_0^{\text{inhom}} \frac{1}{2\sqrt{\pi}\alpha^{1/6}} e^{i\frac{\pi}{4}}. \qquad (B30) $$

From this equation system we obtain

$$ \eta = \frac{\pi}{4} - \frac{2}{3}\alpha \cos^3 \psi, \qquad \Phi_0^{\text{inhom}} = 2\sqrt{\pi}\alpha^{1/6} \cos^{1/2} \psi \Phi_0^{\text{homog}}. \qquad (B31) $$

And the final expression for the electrostatic potential inside a density clump is

$$ \Phi_l^{\text{inhom}}(\bar{z}_0 < \bar{z} < 0) = \Phi_0^{\text{homog}} \sqrt{\pi} \alpha^{1/6} \cos^{1/2} \psi \sqrt{\frac{(-\bar{z})}{3}} \left( e^{-i\pi/6} H_{1/3}^{(2)}\left(\frac{2}{3}(-\bar{z})^{3/2}\right) + e^{i\pi/6} H_{1/3}^{(1)}\left(\frac{2}{3}(-\bar{z})^{3/2}\right) \right) e^{-i\omega t + ik_l \sin \psi x}. \qquad (B32) $$

This solution is obtained under the assumption $(-\bar{z}_0) \gg 1$, as the WKB approximation can be only applied in the wave zone, far from the reflection point. The criterion $(-\bar{z}_0) \gg 1$ may be rewritten as $\alpha^{2/3} \cos^2 \psi \gg 1$, or $\alpha \gg 1$. The specific limitations on the angle of incidence $\psi$ will be discussed in Appendix C. The corresponding electric field components can be calculated from the equation $E^{\text{inhom}} = -\nabla\Phi_l^{\text{inhom}}(\mathbf{r}, t)$ and are

$$ E_{l_x}^{\text{inhom}} = -\frac{d}{dx} \Phi_l^{\text{inhom}}(\mathbf{r}, t) = -ik_{l_x} \Phi_l^{\text{inhom}}(\mathbf{r}, t) = E_{l_{x_i}}^{\text{inhom}} + E_{l_{x_r}}^{\text{inhom}}, \\ = -i \sin \psi \sqrt{\pi} \alpha^{1/6} \cos^{1/2} \psi \sqrt{\frac{(-\bar{z})}{3}} \left( e^{-i\pi/6} H_{1/3}^{(2)}\left(\frac{2}{3}(-\bar{z})^{3/2}\right) + e^{i\pi/6} H_{1/3}^{(1)}\left(\frac{2}{3}(-\bar{z})^{3/2}\right) \right) E_0 e^{-i\omega_{\text{pe}}t + ik_l \sin \psi x} = E_{l_{x_i}}^{\text{inhom}} + E_{l_{x_r}}^{\text{inhom}}, \qquad (B33) $$

$$ E_{l_z}^{\text{inhom}} = -\frac{d}{dz} \Phi_l^{\text{inhom}}(\mathbf{r}, t) = -\alpha^{-1/3} k_l \frac{d}{dz} \Phi_l^{\text{inhom}}(\mathbf{r}, t) \\ = -\sqrt{\pi} \alpha^{-1/6} \cos^{1/2} \psi \sqrt{\frac{(-\bar{z})}{3}} \left( e^{i\pi/6} H_{2/3}^{(2)}\left(\frac{2}{3}|\bar{z}|^{3/2}\right) + e^{-i\pi/6} H_{2/3}^{(1)}\left(\frac{2}{3}|\bar{z}|^{3/2}\right) \right) E_0 e^{-i\omega_{\text{pe}}t + ik_{l_z} \sin \psi x} = E_{l_{z_i}} + E_{l_{z_r}}, \qquad (B34) $$
where we take into account $\Phi_0^{\text{homog}} = E_0^{\text{homog}} k_l^{-1}$. $E_0^{\text{homog}} \equiv E_0$ is the amplitude of the electric field in a homogeneous plasma.

The current that is excited by electron density and velocity perturbations caused by the incident and reflected Langmuir waves should be written in the following form:

$$J_{2\omega_{\text{pe}}}(\mathbf{r}) = -e(\delta n_{l_i}(\mathbf{r})\delta\mathbf{v}_r(\mathbf{r}) + \delta n_{l_r}(\mathbf{r})\delta\mathbf{v}_i(\mathbf{r})), \qquad (B35)$$

where $\delta n_l$ and $\delta\mathbf{v}$ can be expressed from the simplest linear relations:

$$-i\omega_{\text{pe}}\delta\mathbf{v}_{i,r} = -eE_{i,r}/m_e, \qquad (B36)$$

$$-i\omega_{\text{pe}}\delta n_{l,r} + \text{div}(n_0\delta\mathbf{v}) = 0 = -i\omega_{\text{pe}}\delta n_{l,r} + ik_{l,r}n_0\delta\mathbf{v}_{x,r} + n_0\frac{d\delta\mathbf{v}_{z,r}}{dz}, \qquad (B37)$$

$$\delta\mathbf{v}_{i,r} = -ie\frac{E_{i,r}}{\omega_{\text{pe}}m_e}, \qquad (B38)$$

$$\delta n_{l,r} = -\frac{i en_0 k_{l,r}}{\omega_{\text{pe}}^2 m_e} E_{l,x,r} - \frac{en_0}{\omega_{\text{pe}}^2 m_e} \frac{dE_{l,z,r}}{dz}. \qquad (B39)$$

The analytical expression for $J_x$ is

$$\begin{aligned} J_x^{\text{inhom}} &= i \frac{e^3 n_0}{\omega_l^3 m_e^2} \left( (-ik_l \sin \psi E_{x_i} - \frac{\partial E_{z_i}}{\partial z}) E_{x_r} + (-ik_l \sin \psi E_{x_r} - \frac{\partial E_{z_r}}{\partial z}) E_{z_i} \right) \\ &= -\frac{1}{6} \frac{e}{m_e c} \frac{c}{v_b} \cos \psi \sin \psi \alpha^{1/3} (-\bar{z}) (\sin^2 \psi + \alpha^{-2/3} (-\bar{z})) E_0^2 H_{1/3}^{(2)}(\frac{2}{3}(-\bar{z})^{3/2}) H_{1/3}^{(1)}(\frac{2}{3}(-\bar{z})^{3/2}) e^{-2i\omega_{\text{pe}}t + 2ik_{l,r}x}, \end{aligned} \qquad (B40)$$

where we use the fact that Langmuir waves are generated at the local plasma frequency ($\omega_l \approx \omega_{\text{pe}}$) by the electron beam under the resonance condition $k_l \approx \omega_{\text{pe}}/v_b$.

The analytical expression for $J_z$ is

$$\begin{aligned} J_z^{\text{inhom}} &= i \frac{e^3 n_0}{\omega_{\text{pe}}^2 m_e^2} \left( (-ik_l \sin \psi E_{x_i} - \frac{\partial E_{z_i}}{\partial z}) E_{z_r} + (-ik_l \sin \psi E_{x_r} - \frac{\partial E_{z_r}}{\partial z}) E_{z_i} \right) \\ &= i \frac{1}{12} \frac{e}{m_e c} \frac{c}{v_b} \cos \psi (-\bar{z})^{3/2} (\sin^2 \psi + \alpha^{-2/3} (-\bar{z})) E_0^2 \\ &\quad \times \left( H_{-1/3}^{(2)}(\frac{2}{3}(-\bar{z})^{3/2}) H_{2/3}^{(1)}(\frac{2}{3}(-\bar{z})^{3/2}) + H_{-1/3}^{(1)}(\frac{2}{3}(-\bar{z})^{3/2}) H_{2/3}^{(2)}(\frac{2}{3}(-\bar{z})^{3/2}) \right) e^{-2i\omega_{\text{pe}}t + 2ik_{l,r}x}. \end{aligned} \qquad (B41)$$

## Appendix C
Emission from a Single Clump

As soon as we have currents in the form

$$J(\mathbf{r}, t') = J(z)e^{-2i\omega_{\text{pe}}t' + 2ik_{l,r}x}, \qquad (C1)$$

we can calculate the Liénard–Wiechert potential of the harmonic field (Landau & Lifshitz 2013; Jackson 2007),

$$A_{2\omega_{\text{pe}}}(\mathbf{r})e^{-2i\omega_{\text{pe}}t} = \frac{\sqrt{c}}{c} \int \frac{J_{2\omega_{\text{pe}}}(\mathbf{r}, t - \frac{|\mathbf{R} - \mathbf{r}|}{c})}{|\mathbf{R} - \mathbf{r}|} d^3\mathbf{r}, \qquad (C2)$$

and the corresponding current is

$$J_{2\omega_{\text{pe}}}(\mathbf{r}, t - \frac{|\mathbf{R} - \mathbf{r}|}{c}) \approx J_{2\omega_{\text{pe}}}(\mathbf{r}, t - \frac{\mathbf{r}}{c} + \frac{\mathbf{m}\mathbf{n}}{c}) = J_{2\omega_{\text{pe}}}(z)e^{-2i\omega_{\text{pe}}t + 2i\omega_{\text{pe}}R/c + 2i\omega_{\text{pe}}\mathbf{r}\mathbf{n}/c + 2ik_{l,r}x}, \qquad (C3)$$

where $2\omega_{\text{pe}}\mathbf{r}\mathbf{n}/c$ stands for the ratio of the source size to the wavelength of the EM wave, with a factor of $4\pi/\sqrt{3}$. We rewrite $\mathbf{r}\mathbf{n} = x \sin \theta + z \cos \theta$, where $\theta$ is the angle between the wavevector of the EM emission and the z-axis. There are specific constraints for the value of $\psi$ for which the generation of harmonic emission is possible. As along the x-axis, physical parameters do not change according to our assumption, the x-components of the wavevectors of incident and reflected Langmuir waves should be equal $k_{l,r} = k_{l,r'}$, and the momentum conservation of the three-wave interaction should be applied along this axis:

$$k_{l,r} = k_{l,r'} + k_{l,r'}^i, \quad k_{l,r} = \frac{\sqrt{3} \omega_{\text{pe}}}{c} \sin \theta, \quad k_{l,r} = k_{l,r}^i = \frac{\omega_{\text{pe}}}{v_b} \sin \psi. \qquad (C4)$$

This way we obtain the relation between the angles $\psi$ and $\theta$ and a limitation for the angle $\psi$:

$$\sin \psi = \frac{\sqrt{3}}{2} \frac{v_b}{c} \sin \theta, \quad |\psi_{\text{max}}| \approx \frac{\sqrt{3}}{2} \frac{v_b}{c}. \qquad (C5)$$
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Schematic illustration of an integration volume in the vicinity of the reflection point where the Langmuir waves' coalescence into harmonic EM emission takes place. The density clump is approximately represented as a sphere with a characteristic radius $R_s = L\delta n/n_0$, while the integration volume may be roughly described as a cylinder with radius $x_{source}$ and height $z_{source}$, which are the transverse and longitudinal characteristic sizes of the conversion region, i.e., the region where nonlinear currents $J_{2\omega_{pe}}(r, t)$ are excited.</figcaption></figure>

Thus we may rewrite the term proportional to the source size by taking into account $k_t \sin \theta = 2k_t \sin \psi$:

$$2i\omega_{pe} \mathbf{rn}/c = 2i \frac{\omega_{pe}}{c} (x \sin \theta + z \cos \theta) = \frac{4i v_b}{\sqrt{3} c} k_t x \sin \psi + 2i \frac{\omega_{pe}}{c} z \cos \theta, \quad (C6)$$

and the current will have the form

$$J_{2\omega_{pe}}(z) e^{-2i\omega_{pe}t + 2i\omega_{pe}R/c - 2i\omega_{pe}rn/c + 2ik_l x} \approx J_{2\omega_{pe}}(z) e^{-2i\omega_{pe}t + 2i\omega_{pe}R/c - 2i\frac{\omega_{pe}}{c}z \cos\theta + 2ik_l x}. \quad (C7)$$

We use a decomposition for $\frac{1}{|R-r|}$ in spherical coordinates in order to easily separate the quadrupolar term $l=2$:

$$\frac{1}{|R - r|} \approx 4\pi \sum_{l=0}^{\infty} \sum_{m=-l}^{l} \frac{1}{2l + 1} \frac{r^l}{R^{l+1}} Y_l^{m*}(\theta_r, \phi_r) Y_l^m(\theta_R, \phi_R), \quad (C8)$$

and apply it to Equation (C2):

$$A_{2\omega_{pe}}(\mathbf{r}) e^{-2i\omega_{pe}t} = \frac{4\pi\sqrt{\epsilon}}{5c} \frac{1}{R^3} e^{-2i\omega_{pe}t + 2i\omega_{pe}R/c} \sum_{m=-2}^{2} Y_2^m(\theta_R, \phi_R) \int J_{2\omega_{pe}}(z) e^{-2i\frac{\omega_{pe}}{c}z \cos\theta + 2ik_l x} Y_2^{m*}(\theta_r, \phi_r) r^2 d^3r. \quad (C9)$$

It is convenient for us to perform integration (DLMF 2019; Prudnikov et al. 1986) in a cylindrical coordinate system roughly over a cylindrical volume (see Figure 8):

$$A_{2\omega_{pe}}(\mathbf{r}) = \frac{4\pi\sqrt{\epsilon}}{5c} \frac{e^{2i\frac{\omega_{pe}}{c}R}}{R^3} \sum_{m=-2}^{2} Y_2^m(\theta_R, \phi_R) \\ \times \int_0^{z_{\text{source}}} J_{2\omega_{pe}}(z) e^{-2i\frac{\omega_{pe}}{c}z \cos\theta} dz \int_0^{x_{\text{source}}} \rho(\rho^2 + z^2) d\rho \int_0^{2\pi} e^{2ik_x \rho} \cos\phi Y_2^{m*}(\tan^{-1}(\rho/z), \phi) d\phi, \quad (C10)$$

where $z_{source}$ is the distance from the start of the density gradient to the reflection point, i.e., the width of the conversion region, and it should be equal to $3k_l^2 \lambda_D^2 L \cos^2 \psi$ (see Section 3 and Appendix B). As we suppose that the characteristic radius of a spherical density clump is $R_s = L\delta n/n_0$, we may find that under the condition $z_{source} \ll R_s$, which is equivalent to $3k_l^2 \lambda_D^2 \ll \delta n/n_0$, the transverse size of the source region is $x_{source} = \sqrt{6k_l^2 \lambda_D^2/(\delta n/n_0)} L\delta n/n_0$. In order to calculate the magnetic field component of the harmonic emission, we use the following expression:

$$H_{f2\omega_{pe}} = 2i[k_t A_{2\omega_{pe}}], \quad (C11)$$

or

$$|H_{f2\omega_{pe}}| = H_y = 2i(-k_{t_x} A_{z2\omega_{pe}} + k_{t_z} A_{x2\omega_{pe}}). \quad (C12)$$

After calculating all the integrals and keeping only the most significant terms we obtain the $A_x$ component, which has the greatest contribution:

$$A_{x2\omega_{pe}} = -i \frac{\sqrt{\pi \epsilon}}{8} \frac{e^{2iR\omega_{pe}/c - 2i\alpha \cos\theta \cos^2\psi v_b/c + 2i\alpha \cos\theta v_b/c}}{R^3} \frac{e^{c^2 - \frac{1}{v_b^2} \cos\theta}}{m_e c^2 v_b^2 \cos\theta} \frac{1}{k_l^2} \left( \sqrt{\frac{6k_l^2 \lambda_D^2}{n_0}} \frac{\delta n}{n_0} k_l L \right)^3 E_0^2 \\ \times \left[ 2 Y_2^0(\theta_R, \phi_R) J_1\left(2\sqrt{\frac{6k_l^2 \lambda_D^2}{n_0}} k_l L \sin\psi\right) + \sqrt{6} \left(Y_2^{-2}(\theta_R, \phi_R) + Y_2^2(\theta_R, \phi_R)\right) J_3\left(2\sqrt{\frac{6k_l^2 \lambda_D^2}{n_0}} k_l L \sin\psi\right) \right]. \quad (C13)$$
Afterward we calculate the magnetic field component:

$$|H_{t2\omega_{pe}}| \approx 2ik_{t2}A_{x2\omega_{pe}}. \qquad (C14)$$

And finally we can calculate the radiation energy density

$$W_t^{\text{inhom}} = \frac{|H_{t2\omega_{pe}}|^2}{8\pi} \qquad (C15)$$

and integrate it by angles $\theta_R$ and $\phi_R$ using the orthonormality of the spherical harmonics

$$\int_0^\pi \int_0^{2\pi} Y_l^m Y_{l'}^{m'} d\Omega = \delta_{l,l'} \delta_{m,m'}. \qquad (C16)$$

The result has the form

$$W_t^{\text{inhom}} = \frac{1}{8\pi} \frac{3\pi\epsilon}{80} \frac{e^2}{R^6} \frac{c^2}{m_c^2 c^4} \frac{1}{v_b^2 k_l^8} \left( \sqrt{\frac{6k_l^2 \lambda_D^2}{n_0}} \frac{\delta n}{n_0} k_l L \right)^6 \\ \times \left( 4 J_1 \left( 2 \sqrt{\frac{6k_l^2 \lambda_D^2}{n_0}} \frac{\delta n}{n_0} k_l L \sin \psi \right) + 12 J_3 \left( 2 \sqrt{\frac{6k_l^2 \lambda_D^2}{n_0}} \frac{\delta n}{n_0} k_l L \sin \psi \right) \right) E_0^4. \qquad (C17)$$

We want to simplify it slightly by using limiting forms of the Bessel functions for a small value of their argument (DLMF 2019; this approximation is always correct for rather small values of the angle $\psi$):

$$4J_1^2\left(2\sqrt{\frac{6k_l^2 \lambda_D^2}{n_0}} \frac{\delta n}{n_0} k_l L \sin \psi\right) + 12J_3^2\left(2\sqrt{\frac{6k_l^2 \lambda_D^2}{n_0}} \frac{\delta n}{n_0} k_l L \sin \psi\right) \approx 24 \frac{\delta n}{n_0} k_l^2 L^2 k_l^2 \lambda_D^2 \sin^2 \psi. \qquad (C18)$$

We put the observation point $R$ at the border of a density clump to estimate the emission that is detected when it leaves the source region. Since we imply $z_{\text{source}} \ll R_s$, this approximation that is used for the decompositions in Equations (C3) and (C8) is valid and we may set $R \approx R_s = L\delta n/n_0$. Taking into account the approximate expression for Bessel functions and $k_l^2 \lambda_D^2 \approx v_T^2/v_b^2$ for beam-generated Langmuir waves, after a few simple transformations we obtain an expression for the energy density of harmonic emission for a single density clump:

$$W_t^{\text{inhom}} = 4 \cdot 10^2 \pi \epsilon \frac{v_T^6}{v_b^6} \left( \frac{k_l^2 \lambda_D^2}{\delta n / n_0} \right)^2 \frac{\omega_{pe}^2 L^2}{c^2} \sin^2 \psi \frac{W_l}{n_0 k_B T_e} W_l. \qquad (C19)$$

## Appendix D
### Statistically Averaged Emission

In order to account for different angles of incidence, amplitudes of random density fluctuations, and gradient scales, we perform statistical averaging over the $\psi$, $\delta n$, and $L$ probability density functions:

$$\langle W_t^{\text{inhom}} \rangle_{\psi, \delta n, L} = P_{\text{ref}} \int P(\psi) P(\delta n) P(L) W_t^{\text{inhom}}(\psi, \delta n, L) d\psi d\delta n dL. \qquad (D1)$$

We assume a uniform distribution over the angles of incidence, $P(\psi) = \pi^{-1}$. In order to perform the integration we recall that the value of $\psi$ is limited by $\psi_{\text{max}} \approx \sqrt{3} v_b / 2c \ll 1$:

$$\langle W_t^{\text{inhom}} \rangle_{\psi} = \dots \frac{1}{\pi} \int_0^{\psi_{\text{max}}} \sin^2 \psi d\psi \approx \dots \frac{1}{\pi} \int_0^{\psi_{\text{max}}} \psi^2 d\psi = \dots \frac{1}{3\pi} \psi_{\text{max}}^3 \approx \frac{\sqrt{3} v_b^3}{8\pi c^3}. \qquad (D2)$$

We assume that fluctuations follow a normal distribution with zero mean value and standard deviation $\langle \Delta n \rangle$:

$$P_{\delta n}(\delta n) = \frac{1}{\sqrt{2\pi} \langle \Delta n \rangle} \exp \left[ -\frac{\delta n^2}{2 \langle \Delta n \rangle^2} \right]. \qquad (D3)$$

Then, averaging over the amplitudes of density fluctuations can be performed as follows:

$$\langle W_t^{\text{inhom}} \rangle_{\delta n} = \dots \frac{n_0^2}{\sqrt{2\pi} \langle \Delta n \rangle} \int_0^{\infty} \delta n^{-2} \exp \left[ -\frac{\delta n^2}{2 \langle \Delta n \rangle^2} \right] d\delta n = \dots - \sqrt{2} \left( \frac{\langle \Delta n \rangle}{n_0} \right)^{-2}. \qquad (D4)$$

Let us assume that all fluctuations have the same size $P_{\text{ref}}$. In this case, the probability of finding a fluctuation with density variation $\delta n$ should be equal to the probability of finding a fluctuation with gradient $L$:

$$P_{\delta n}(\delta n) d\delta n = P_L(L) dL. \qquad (D5)$$
Since $\delta n/n_0 = P/L$, we may write

$$ \frac{\partial \delta n}{\partial L} = - \frac{n_0 P}{L^2}, \qquad (D6) $$

and

$$ P_L(L)dL = P_{\delta n}(\delta n(L)) \left| \frac{\partial \delta n}{\partial L} \right| dL = \frac{1}{\sqrt{2\pi}} \frac{n_0}{\langle \Delta n \rangle L^2} \frac{P}{L^2} \exp \left[ - \frac{P^2 n_0^2}{2 \langle \Delta n \rangle^2 L^2} \right] dL, \qquad (D7) $$

or after substitution of $P/(\langle \Delta n \rangle/n_0) = L_{sc}$,

$$ P_L(L)dL = \frac{1}{\sqrt{2\pi}} \frac{L_{sc}}{L^2} \exp \left[ - \frac{L_{sc}^2}{2L^2} \right] dL. \qquad (D8) $$

Now we can perform averaging over L:

$$ \left\langle W_t^{\text{inhom}} \right\rangle_L = \dots \frac{L_{sc}}{\sqrt{2\pi}} \int_0^\infty \exp\left[-\frac{L_{sc}^2}{2L^2}\right] dL = \dots \frac{L_{sc}^2}{2}. \qquad (D9) $$

Putting them all together, we obtain

$$ \left\langle W_t^{\text{inhom}} \right\rangle_{\psi, \delta n, L} = 25\sqrt{6} \epsilon P_{\text{ref}} \frac{v_T^3 v_T^3}{c^3 v_B^3} \left( \frac{k_T^2 \lambda_D^2}{\langle \Delta n \rangle / n_0} \right)^2 \frac{\omega_{P_e}^2 L_{sc}^2}{c^2} \frac{W_f}{n_0 k_B T_e} W_L. \qquad (D10) $$

## ORCID iDs

Anna Tkachenko https://orcid.org/0000-0003-3705-9690

Vladimir Krasnoselskikh https://orcid.org/0000-0002-6809-6219

Andrii Voshchepynets https://orcid.org/0000-0001-8307-781X

## References

* Bian, N. H., Kontar, E. P., & Ratcliffe, H. 2014, JGRA, 119, 4239

* Brejzman, B., & Pekker, L. 1978, PhLA, 65, 121

* Cairns, I. H. 1987a, JPIPh, 38, 179

* Cairns, I. H. 1987b, JPIPh, 38, 199

* Celnikier, L., Muschietti, L., & Goldman, M. 1987, A&A, 181, 138

* Chen, C., Salem, C., Bonnell, J., Mozer, F., & Bale, S. 2012, PhRvL, 109, 035001

* Chen, F. F. 1984, Introduction to Plasma Physics and Controlled Fusion, Vol. 1 (Boston, MA: Springer)

* Chen, X., Kontar, E. P., Yu, S., et al. 2018, ApJ, 856, 73

* DLMF 2019, NIST Digital Library of Mathematical Functions, Release 1.0.25, http://dlmf.nist.gov/

* Drummond, W., & Pines, D. 1964, AnPhy, 28, 478

* Dulk, G., & Suzuki, S. 1980, A&A, 88, 203

* Dulk, G. A., Leblanc, Y., Robinson, P. A., Bougeret, J.-L., & Lin, R. P. 1998, JGR, 103, 17223

* Ergun, R., Larson, D., Lin, R., et al. 1998, ApJ, 503, 435

* Ergun, R., Malaspina, D., Cairns, I. H., et al. 2008, PhRvL, 101, 051101

* Erokhin, N., Moiseev, S., & Mukhin, V. 1974, NucFu, 14, 333

* Galeev, A., & Krasnoselskikh, V. 1976, ZhPmR, 24, 558

* Ginzburg, V., & Zheleznyakov, V. 1958, SvA, 2, 653

* Goldman, M. V., Reiter, G. F., & Nicholson, D. R. 1980, PhFl, 23, 388

* Goldstein, M. L., Roberts, D. A., & Matthaeus, W. 1995, ARA&A, 33, 283

* Gurnett, D. A., & Frank, L. A. 1978, JGR, 83, 58

* Harding, J. C., Cairns, I. H., & Melrose, D. B. 2020, PhPl, 27, 020702

* Hinkel-Lipsker, D., Fried, B., & Morales, G. 1992, PhFIB, 4, 559

* Jackson, J. D. 2007, Classical Electrodynamics (New York: Wiley)

* Kellogg, P. 1980, ApJ, 236, 696

* Kellogg, P., Goetz, K., Monson, S., & Bale, S. 1999, JGRA, 104, 17069

* Kellogg, P. J., & Horbury, T. 2005, AnGeo, 23, 3765

* Kim, E.-H., Cairns, I. H., & Johnson, J. R. 2013, PhPl, 20, 122103

* Kim, E.-H., Cairns, I. H., & Robinson, P. A. 2007, PhRvL, 99, 015003

* Kim, E.-H., Cairns, I. H., & Robinson, P. A. 2008, PhPl, 15, 102110

* Kim, E.-H., Johnson, J. R., Cairns, I. H., & Lee, D.-H. 2009, in AIP Conf. Proc. 1187, Radio Frequency Power in Plasmas: Proceedings of the 18th Topical Conference, ed. V. Bobkov & J.-M. Noterdaeme (Melville, NY: AIP), 13

* Kontar, E. P., Chen, X., Chrysaphi, N., et al. 2019, ApJ, 884, 122

* Krafft, C., & Volokitin, A. 2014, EPJD, 68, 370

* Krafft, C., Volokitin, A., & Krasnoselskikh, V. 2013, ApJ, 778, 111

* Krasnoselskikh, V. V., Dudok de Wit, T., & Bale, S. D. 2011, AnGeo, 29, 613

* Krasnoselskikh, V. V., Lobzin, V. V., Musatenko, K., et al. 2007, JGRA, 112, A10109

* Krasnoselskikh, V., Voshchepynets, A., & Maksimovic, M. 2019, ApJ, 879, 51

* Krucker, S., Oakley, P. H., & Lin, R. P. 2009, ApJ, 691, 806

* Krupar, V., Szabo, A., Maksimovic, M., et al. 2020, ApJS, 246, 57

* Landau, L. D., & Lifshitz, E. M. 2013, Course of Theoretical Physics (Amsterdam: Elsevier)

* Lin, R., Levedahl, W., Lotko, W., Gurnett, D., & Scarf, F. 1986, ApJ, 308, 954

* Lin, R., Potter, D., Gurnett, D., & Scarf, F. 1981, ApJ, 251, 364

* Malaspina, D. M., Cairns, I. H., & Ergun, R. E. 2012, ApJ, 755, 45

* Malaspina, D. M., & Ergun, R. E. 2008, JGRA, 113, A12108

* Mann, G., Breitling, F., Vocks, C., et al. 2018, A&A, 611, A57

* Melrose, D. B. 1980a, SSRv, 26, 3

* Melrose, D. B. 1980b, Plasma Astrophysics: Nonthermal Processes in Diffuse Magnetized Plasmas, Vol. 2 (New York: Gordon and Breach)

* Melrose, D. B. 1987, SoPh, 111, 89

* Meyer-Vernet, N., & Issautier, K. 1998, JGR, 103, 29705

* Mjølhus, E. 1983, JPIPh, 30, 179

* Mjølhus, E. 1990, RaSc, 25, 1321

* Neugebauer, M. 1975, JGR, 80, 998

* Nishikawa, K., & Ryutov, D. 1976, JPSJ, 41, 1757

* Papadopoulos, K., & Freund, H. 1978, GeoRL, 5, 881

* Prudnikov, A. P., Brychkov, Y. A., & Marichev, O. I. 1986, Integrals and Series: Special Functions, Vol. 2 (Boca Raton, FL: CRC Press)

* Ratcliffe, H., Bian, N. H., & Kontar, E. P. 2012, ApJ, 761, 176

* Reid, H. A., & Kontar, E. P. 2010, ApJ, 721, 864

* Reid, H. A. S., & Ratcliffe, H. 2014, RAA, 14, 773

* Robinson, P., & Cairns, I. 1998, SoPh, 181, 363

* Schleyer, F., Cairns, I. H., & Kim, E.-H. 2013, PhPl, 20, 032101

* Schleyer, F., Cairns, I. H., & Kim, E.-H. 2014, JGRA, 119, 3392

* Shaikh, D., & Zank, G. 2010, MNRAS, 402, 362

* Sturrock, P. 1964, NASSP, 50, 357

* Thejappa, G., MacDowall, R., & Kaiser, M. 2007, ApJ, 671, 894

* Tsytovich, V. 2012, Nonlinear Effects in Plasma (New York: Springer),

* Vedenov, A. A., Velikhov, E. P., & Sagdeev, R. Z. 1962, Quasi-linear Theory of Plasma Oscillations (Vienna: International Atomic Energy Agency)

* Volokitin, A. S., & Krafft, C. 2016, ApJ, 833, 166

* Volokitin, A. S., & Krafft, C. 2018, ApJ, 868, 104
Volokitin, A. S., & Krafft, C. 2020, ApJL, 893, L47

Voshchepynets, A., & Krasnoselskikh, V. 2013, AnGeo, 31, 1379

Voshchepynets, A., & Krasnoselskikh, V. 2015, JGRA, 120, 10139

Voshchepynets, A., Krasnoselskikh, V., Artemyev, A., & Volokitin, A. 2015, ApJ, 807, 38

Voshchepynets, A., Volokitin, A., Krasnoselskikh, V., & Krafft, C. 2017, JGRA, 122, 3915

Wasow, W. 2018, Asymptotic Expansions for Ordinary Differential Equations
(New York: Courier Dover)

Willes, A., Robinson, P., & Melrose, D. 1996, PhPl, 3, 149

Zakharov, V. E. 1972, JETP, 35, 908

Zheleznyakov, V., & Zaitsev, V. 1970, SvA, 14, 250

Ziebell, L., Gaelzer, R., Pavan, J., & Yoon, P. 2008, PPCF, 50, 085011

Ziebell, L., Yoon, P., Pavan, J., & Gaelzer, R. 2011, PPCF, 53, 085004