Particle-in-Cell Simulations of High-Power THz Generator Based on the Collision of Strongly Focused Relativistic Electron Beams in Plasma

Based on particle-in-cell simulations, we propose to generate sub-nanosecond pulses of narrowband terahertz radiation with tens of MW power using unique properties of kiloampere relativistic (2 MeV) electron beams produced by linear induction accelerators. Due to small emittance of such beams, they can be focused into millimeter and sub-millimeter spots comparable in sizes with the wavelength of THz radiation. If such a beam is injected into a plasma, it becomes unstable against the two-stream instability and excites plasma oscillations that can be converted to electromagnetic waves at the plasma frequency and its harmonics. It is shown that several radiation mechanisms with high efficiency of power conversion (∼1%) come into play when the radial size of the beam–plasma system becomes comparable with the wavelength of the emitted waves.


Introduction
In recent years, the efforts of many scientific groups have been aimed at creating radiation sources in the terahertz (THz) frequency range (0.3-10 THz). The reason for this is a large number of different applications of such radiation in science and technology (see reviews [1,2] and references therein). The leading role in generating narrowband THz radiation is presently played by large and expensive accelerating facilities such as free-electron lasers [3][4][5]. The generation of high-power narrowband THz radiation using more compact devices continues to be a challenge for several reasons. In particular, the conversion of laser radiation using nonlinear crystals [6][7][8] faces the problem of changing the properties and even destruction of a nonlinear medium when going to high powers. In the region of less than 1 THz, a lot of success has been achieved by devices based on vacuum electronics, for example, backward wave oscillators (BWO) [9]. For development in this area, breakthrough technological solutions are needed to create channels of a very small cross-section with a very high surface quality, without unnecessary impurities. In addition, it is also necessary to generate beams with very high currents and very small cross sections (cathodes, corresponding optics). The structural elements of such devices are often of the order of the wavelength, which can lead to overheating at high powers. Increasing the frequency in devices based on the use of magnetic fields often requires very high field values. In recent years, great progress in filling the THz gap has been made by gyrotrons [10][11][12]. However, an increase in the radiation frequency in such devices is associated with the need to use magnetic fields of tens of Tesla. In addition, transition to high frequencies in this case is accompanied by a significant decrease in maximum power and radiation efficiency.
The problem of low thresholds for the destruction of a nonlinear medium is naturally overcome in plasma. Plasma allows to maintain electromagnetic oscillations of very large amplitude, the frequency of which is tied to the plasma density. If we ensure their conversion into electromagnetic radiation of the THz range, then, on the one hand, we will avoid the need to use superstrong magnetic fields to generate THz radiation, and on the other hand, we will be able to change the radiation frequency by simply changing the plasma density. Such advantages of plasma have stimulated the development of plasmabased THz radiation sources, most of which are based on laser-plasma systems [13][14][15][16][17][18][19].
A promising way in the development of powerful narrowband THz radiation sources is the use of kiloampere electron beams. The power of such beams reaches several gigawatts; therefore, conversion into radiation of even a percent of this power makes it possible to produce electromagnetic (EM) emission at the level of hundreds of MW. In addition, such beams have a sufficiently long duration to provide a narrow radiation spectrum. Among the THz generation schemes developed in this direction, it is worth noting a powerful planar FEL driven by a large-size sheet electron beam forming by the high-current accelerator «ELMI» with beam energy 1 MeV, current 7 kA and duration 3 µs [20,21]. The possibility of constructing a free-electron maser based on high-current electron beams generated in linacs is also being considered [22,23].
Long-pulse electron beams can also be effectively used to generate high-power EM radiation. In experiments on the injection of a sub-relativistic electron beam into plasma at the GOL-3 facility (BINP SB RAS), an unexpectedly high level of electromagnetic emission has been observed when the beam diameter decreases to a size of the order of the radiation wavelength [24]. To explain the experimental results, the theory of the beam-driven plasma antenna [25,26] has been proposed. Using this model, it has become possible not only to shed light on the main mechanism of generation of EM waves in these beamplasma experiments, but also to propose new schemes for highly efficient generation of EM radiation at harmonics of the plasma frequency in a plasma with a pre-modulated density [27,28]. At present, experiments on the generation of radiation by relativistic beams in a plasma with regular density inhomogeneities [29] are being carried out at the GOL-PET facility (BINP SB RAS). In addition to processes in a single beam-plasma system, a significant increase in the level of electromagnetic emission is possible in the presence of colliding electron beams in plasma. Significant enhancement of the second harmonic emission is observed not only in the case of beams that excite counterpropagating plasma oscillations with different transverse profiles of electrostatic potential [30], but also with a special selection of the parameters of the beam-plasma system, at which three-wave interactions between the most unstable plasma oscillations and electromagnetic waves are realized at the linear stage of the beam-plasma instability [31].
In this paper, using particle-in-cell (PIC) simulations, we study the possibility to generate high-power THz radiation in plasma by the counterstreaming electron beams with a few kA currents and relativistic energies ∼2 MeV typical to the linear induction accelerators [32]. In Section 2, we describe our 2D3V axially symmetric PIC model with open boundary conditions allowing for the steady-state injection of electron beams into a plasma. In Section 3, simulation results on the injection of both a single beam and two oppositely propagating beams are presented. This section also discusses the methodology for processing simulation results and analyzes them. A summary of the results obtained is presented in Section 4.

Numerical Model
To simulate the problem, we use the PIC method in which distribution functions of real particles is replaced by the ones of model particles Here, r j , v j are the coordinate and velocity of the particle with the index j. Model particles have a spatial charge distribution characterized by the shape function S(r, r j ) ρ(r, t) = ∑ j q j S(r, r j (t)), (2) and interact with each other only through the computational grid thereby contributing to the charge and current density.
To solve the equations of particle motion, we use the scheme of Boris [33]. This is an explicit scheme in which the motion of a particle is split into separate motion in the electric field and in the magnetic field. Electric currents are calculated using the scheme of Esirkepov [34] adapted for the cylindrical geometry. In this scheme, at each step the charge density of the particle at the start and end points of movement is calculated by Formula (2) for a given form factor, and then used in the difference analogue of the continuity equation to obtain the current density. It automatically satisfies the difference continuity equation and, therefore, accurately fulfils the Gaussian difference law. The charge density of model particles are described by the fourth-order shape function: Here, d = |r − r j | and h is the grid step. We have chosen the fourth-order of the form factor in order to decrease numerical noise associated, among other things, with the numerical Cherenkov instability.
To solve Maxwell's equations, we use the finite-difference scheme (FDTD) proposed in [35], which provides the second order of accuracy in time and space. The problem is considered in the axially symmetric r-z geometry, therefore the particles and grid cells are represented by uniform rings.
The computational domain contains a continuous plasma column bounded in the radial direction by vacuum ( Figure 1). The dimensions of the area are chosen equal to L z = 100c/ω p and L r = 32c/ω p , the grid and time steps are h r = h z = 0.02c/ω p and τ = 0.01ω −1 p . The plasma width is L plasma = 20c/ω p . We use 49 model particles per cell for each particle species. In the z-direction, we set open boundary conditions described in [36]. Outside the area, we use the perfectly matched layers (PML) [37] in all directions to eliminate reflections of the produced EM waves.
We assume that the particles of the injected beams moving without the action of external forces would have a Gaussian distribution of density n(x, y) and momenta f (p x , p y ) in the focal region z = 0: The dispersions σ r , σ p of these distributions determine the focal spot sizes and focusing angles. The quantity n b defines the beam density. Along the z-axis, beam particles have the same momentum p z . Therefore, to simulate the injection of focusing beams, at each step, we first set the required distribution of particles at the focus in 3D3V space (for this purpose, it is sufficient to consider only one cell in the z-direction), after which we move the particles with a given momentum without external forces to the boundary at which injection will take place ( Figure 2): n(x, y, z) → n (x , y , z ) :   Here, L = L z /2 is the distance from the focus to the border. Then we make the transition from 3D3V (x , y , z , p x , p y , p z ) geometry to the axially symmetric system (r , z , p r , p φ , p z ) and inject the resulting particles into the simulation region, having previously changed the direction of the pulse: The momentum spread σ p = R L p z is acquired due to focusing from the radius R on the distance L. These dimensions are dictated by the intended experiment. The resulting radial distribution of beam density (Figure 3) is described by the formula: In our calculations, we consider two types of beams. The focal size of the first beam is σ (1) r ≡ σ 1 = 3.325c/ω p , while its peak density is determined by the value n (1) b = 0.02n 0 , where n 0 is the plasma density. We will call this beam as a "thin" one. The second beam ("wide") has σ (2) r ≡ σ 2 = 6.65c/ω p and n (2) b = 0.01n 0 . The energy of beam particles is chosen equal to 2 MeV which corresponds to the velocity v b = 0.979c (γ ≈ 4.9). In this case, the thin beam is characterized by the power P 1 = 3.665 GW and current I 1 = 2.9 kA. For the wide beam, we obtain P 2 = 7.33 GW and I 2 = 5.8 kA. Further we will operate with dimensionless quantities. So all lengths will be measured in units of c/ω p , times in ω −1 p , electromagnetic fields in m e cω p /e. For the convenience of conversion into dimensional quantities, Figure 4 shows the dependence of these quantities, as well as the linear frequency of radiation at the harmonics of the plasma frequency, on the initial plasma density n. In particular, for the density n 0 = 1.2512 × 10 15 cm −3 : • the unit of length is c/ω p = 150 µm and the focal spot-size is σ 2 = 0.1 cm; • the radiation frequency at the third harmonic of plasma frequency corresponds to • the unit of time ω −1 p equals to 0.5 ps, and the total duration of a simulation run is 1000ω −1 p = 0.5 ns; • the typical dimensionless electric field amplitude in radiated waves constitutes E rad = 0.005 with the maximum E rad = 0.01, which in dimensional units correspond to 0.17 MV/cm and 0.34 MV/cm.

Single Beam Injection
Let us first consider the process of focusing a single beam into a spot of the width σ 1 = 3.325 c/ω p located in the center of the system. Figure 5 shows how the beam density and the corresponding profile of longitudinal electric field in the unstable plasma wave excited by the beam evolve in time. When the beam is injected, the intense two-stream instability develops in the focal region. This leads to the formation of a modulated beam structure already by the time tω p = 150. In this case, the characteristic distance between the forming bunches corresponds to the wavelength of the most unstable mode (5.5 c/ω p ). The location of beam focus in plasma, however, does not remain stationary and demonstrates the complicated dynamics.
Initially, the most intense excitation of plasma oscillations occurs at the geometrical focus of the beam, where it reaches the highest density. Later, however, one can see that the beam is focused closer to the injector (tω p = 300). The reason for this refocusing is that a plasma compensates the charge of the beam which, in vacuum, tends to push particles from the axis via the radial electric field, but do not completely compensate the current which compresses the beam radially by its own azimuthal magnetic field. As a result, the beam turns out to be compressed to high densities before it reaches the focus that would be realized in a vacuum. Over time, the plasma wave becomes less regular and the bunch structure of the beam is less pronounced. Let us study spectral characteristics of beam-driven plasma oscillations in more detail. For this, we carry out a two-dimensional Fourier transform of the dependence of the longitudinal electric field E z (z, t) near the axis of the system on the longitudinal coordinate and time. Figure 6 shows the resulting spectrum. The dominant beam-driven mode has a longitudinal wavenumber |k | = |k 0 | ≈ 1.06 ω p /c and frequency ω 0 ≈ 0.96 ω p . A large amplitude of plasma oscillations leads to the development of nonlinear processes accompanied by excitation of harmonic oscillations (nω 0 , nk 0 ). The frequency of plasma waves throughout the simulation remains constant with good accuracy, as evidenced by the narrow line-width in the Fourier spectrum. Contrary, the spread in wavenumbers is more significant.
Electromagnetic emission in this case is found to be very weak. We only observe negligible electromagnetic radiation at harmonics 2ω 0 and 3ω 0 with the efficiency of beamto-radiation power conversion at the level of 10 −2 % and 3 × 10 −3 %, respectively. The situation is not changed drastically even for mobile ions. In the case of a thin beam, one can expect the switching on of the highly efficient mechanism of plasma antenna when the ion density get an opportunity to be modulated in the longitudinal direction by largeamplitude plasma waves. But longitudinal modulation in the case of a single beam turns out to be a slower process than radial displacement of plasma from the beam region by the ponderomotive force of plasma oscillations. Emission can be enhanced only if we introduce either the guiding magnetic field eliminating the radial plasma dynamics or a counterpropagating beam that makes possible to generate radiation by the collision of different-size waves. Indeed, if the thin beam is guided by the strong magnetic field B z = 1.1 (≈12.5 T for the plasma density n 0 = 1.25 × 10 15 cm −3 ), the maximum value of the total radiation power reaches 1.6% of the power injected by the beam (Figure 7). EM radiation appears when the ion density is modulated in the axial direction with the period appropriate for the antenna mechanism. It is also seen that the produced radiation is strongly dominated by the fundamental frequency ω 0 .

Processing Technique of Radiation Signals
Before proceeding to a discussion of simulation results obtained for the two-beam systems, let us consider the technique we use for processing the recorded radiation signals. Figure 8 (top) shows a part of the time dependence of a typical field signal E z recorded at one spatial point on the boundary of the computational domain. Due to the significant nonlinearity of the oscillations excited by the beams, the radiation contains not only the second, but also higher harmonics. In order to investigate their relative contributions separately, we use Fourier filtering. For each radiation signal, we carry out a direct Fourier transform and obtain its spectrum ( Figure 8 on the left). Then we carry out the inverse Fourier transform separately for parts of the spectrum in the range nω p ± 0.5 ( Figure 8 on the right). Figure A1 in the Appendix A shows the temporal dependences of the field amplitude E z at points on the boundary of the vacuum region for each radiation harmonic in each of the simulated two-beam systems. A similar procedure is carried out for the B φ field. Knowing the fields of individual harmonics E z and B φ on the boundary of the vacuum system, one can find their contribution to the normal component of the Poynting vector E z · B φ , from which we obtain the power of radiation leaving the system in the radial direction at a given harmonic (Figure 8 bottom). By integrating its value along the boundary of the vacuum region at each moment of time, we obtain the temporal dependence of the total power radiated in all directions. The ratio of the obtained value to the power of the injected beams will give us the efficiency of radiation generation.

Counter Injection of Two Symmetric Thin Beams
Let us study the process of radiation generation in a plasma with two colliding beams. Let us first consider the case of two thin identical beams, as well as the case of asymmetric collision of thin and wide beams in a plasma with immobile (infinitely heavy) ions. Then, we estimate the effect of mobile hydrogen ions with a real mass on the radiation process in both cases.
Let two colliding electron beams be focused to the center of the system into spots with the same sizes σ 1 = σ 2 = 3.325 c/ω p . Figure 9 shows the spectrum of oscillations excited by the beams. In addition to oscillations (nω 0 , ±nk 0 ) arising from each beam, oscillations at the combined frequencies and wavenumbers are observed in the spectrum. Resonance with vacuum EM waves can be achieved only for those plasma oscillations that travel along the plasma column with superluminal phase velocities. For this reason, neither the most unstable modes (ω 0 , ±k 0 ) directly excited by each of the beams at the Cherenkov resonance, nor their harmonics can be responsible for EM emission. Superluminal satellites in plasma can be excited only by nonlinear interaction of counterpropagating beam modes. For example, radiation at the second harmonic is generated due to the coalescence (ω 0 , k 0 ) + (ω 0 , −k 0 ) → (2ω 0 , 0), while the third and forth harmonics arise in the processes: (ω 0 , ∓k 0 ) + (2ω 0 , ±2k 0 ) → (3ω 0 , ±k 0 ), (18) The longitudinal phase velocities of the oscillations involved in these nonlinear interactions are given in Table 1. The process (17), however, degenerates if the dominant beam modes are excited by symmetric drivers and run along the beam axis. Indeed, the symmetry of the system in this case is so high that it leads to the vanishing of the electric current quadratic in the field. As shown in recent works [15,30], the asymmetry necessary for the generation of radiation in this case can be provided due to differing profiles of the electrostatic potential in colliding waves Φ 1 (r ⊥ ) = Φ 2 (r ⊥ ). The generating current then has the form and the radiation it produces is directed strictly across the axis of the system. Knowing the frequency and longitudinal wavenumber of oscillations, as well as the dispersion law of electromagnetic oscillations in vacuum ω = kc, one can find the angle at which radiation should be generated in the processes (18)- (22): Here, k ⊥ is the transverse wavenumber of emitted EM waves. For emission at the frequency 3ω 0 , Then the angle of radiation is Similarly, we obtain the angle of the fourth harmonic radiation: The frequencies of the plasma oscillations responsible for the radiation are determined unambiguously, while the wavenumbers have a certain scatter, which should lead to an uncertainty in the radiation angle. Further, based on the results shown in Figure 6, we will consider the longitudinal wavenumber k to be in the range [0.95; 1.2]. Thus, the third and forth harmonic radiations are predicted to be emitted at angles Let us focus on the radiation processes observed in PIC simulations in such a beamplasma system. Figure 10 shows simulation results in different moments of time. At the initial stage (t = 200ω −1 p ), the particles of both beams are focused into the center where the largest amplitude of plasma oscillations is achieved. It can be seen that radiation is generated in the center of the system at the third and fourth harmonics. The radiation angles in this case correspond to those theoretically predicted in (27). In the given moment of time, the local transverse profiles of plasma oscillations excited by each beam coincide with good accuracy over the entire length of the system, and radiation at the second harmonic is not observed. Later (t = 262ω −1 p ), the beams are refocused in the same manner as discussed in the single beam case and the regions of intense relaxation for each of the beams are shifted towards their own injectors. Two spaced regions with a high amplitude of plasma oscillations are formed. Each region becomes a source of radiation at the third and fourth harmonics. By the time moment t = 345ω −1 p , the relaxation regions of the beams are significantly removed from the center of the system. In each given region, the counterpropagating beam-driven waves have different transverse structures. This is a necessary condition for switching on the mechanism of radiation generation near the second harmonic of the plasma frequency (23). The appearance of such radiation is clearly observed in Figure 10. In the center of the system, a clear zero of the radiation field is visible, since this mechanism does not work here due to the symmetry of the plasma waves. In the main relaxation zone of each beam, the excitation of counterpropagating oscillations by the opposite beam is reduced. This leads to a relative decrease in the amplitude of nonlinear harmonic oscillations in plasma and to a decrease in the radiation power at frequencies 3ω 0 and, in particular, 4ω 0 . Let us analyze the frequency spectrum of the generated radiation. Figure 11 shows the Fourier spectrum of radiation detected at the boundary of the vacuum region and the spatial distribution of the energy emitted at each harmonic during the entire simulation time. It is seen that each harmonic is characterized by a narrow spectral line (∆ω/ω ∼ 1%), and the dominant role is played by the second harmonic radiation. The spatial profile of radiation intensity at the frequency 4ω 0 , integrated over the entire simulation time, corresponds to the early stage in Figure 10. The fourth harmonic emission is the result of high-order nonlinearity which requires high amplitude of plasma oscillations and is efficiently realized at the early stages of the process, when the relaxation regions of each beam are well overlapped. The spatial distribution of the third harmonic is more uniform, while the second harmonic radiation demonstrates the pronounced peaks from the newly formed foci of colliding beams.
Note that distribution functions of beams and plasma in numerical simulations are described using a finite number of model particles. To model a finite temperature, one should set a scatter in particle momenta using a random number generator. For these reasons, all numerical simulations, even with the same macro parameters, may differ in details from each other. Since the beam-plasma instability grows from these noise fields, dynamics of counterpropagating unstable waves is not fully symmetric. It explains the asymmetric spatial distribution for the second harmonic emission that is observed in Figure 11. The issue of repeatability for numerical simulations of beam-plasma interaction has been discussed in more details in [31].

Counter Injection of Different-Size Beams
Let us now consider the generation of radiation by colliding beams focused into different-size spots. The size of the left beam remains equal to σ 1 = 3.325 c/ω p , while the right beam is chosen twice wider σ 2 = 6.65 c/ω p . Results of simulations of such a system in different moments of time are presented in Figure 12. In this case, colliding beams excite plasma oscillations with different transverse structures from the very beginning; therefore, in contrast to the case of identical thin beams, the second harmonic emission here is observed immediately after the development of the two-stream instability. EM waves at higher harmonics are also visible at the initial stage, but they are mainly emitted along the direction of thin beam propagation. It is explained by the dominance at this stage of a plasma wave excited by the thin beam. Further development of two-stream instability leads to an increase in the amplitude of plasma waves. The radiation at high harmonics is amplified and begins to significantly exceed the fraction of radiation at the second harmonic. It is interesting that, by the later time moment t = 580ω −1 p , this radiation changes direction and begins to significantly dominate over other emission processes. Figure 13 shows the spectrum of the generated EM radiation and the spatial distribution of its energy accumulated during the entire simulation time. According to the theory, the radiation at the frequency 2ω 0 was expected to dominate due to the interaction of counterpropagating plasma oscillations with different transverse structures. However, such radiation was effectively generated only at the initial stage of interaction. At later times, the greatest contribution to the total radiation power is made by the third harmonic.

Ion Effects
Let us find out how the parameters of EM emission at the plasma frequency harmonics will change if hydrogen ions with a real mass (m i = 1836m e ) are used in the scheme with two thin beams. Results of PIC simulations in this case are presented in Figure 14. It can be seen from the ion density maps that, in addition to the displacement of the plasma along the radius, the formation of longitudinal modulation of the plasma density with the wavenumber 2k 0 is also observed. In contrast to the case of a single beam, the longitudinal modulation in the system of counterstreaming beams grows much faster due to the fact that superposition of colliding beam-driven oscillations form a standing plasma wave. The ponderomotive force in such a wave is modulated in space, therefore, by displacing electrons from the nodes of the wave and creating an uncompensated charge, it leads to modulation of the ion density. At the early stages, when the ion density perturbation can be considered as a small one, its presence should enhance the radiation at the second harmonic, since each beam-driven wave is able to radiate by the additional mechanism of the plasma antenna [26,27,38]. However, already by the time t = 220ω −1 p , the ion modulation becomes so deep that it significantly weakens the beam instability [39], leading to a significant decrease in the radiation intensity at all harmonics ( Figure 15).
Thus, the duration of the radiation during counter-beam injection is limited to the times at which the ion dynamics becomes significant. For plasma densities n 0 ≈ 1.2 × 10 15 cm −3 this duration is a fraction of a nanosecond. With increasing ion mass (in this work, we have considered the lightest ions-hydrogen nuclei), the duration of radiation can be increased by increasing the duration of beam injections.

Radiation Efficiency
Let us estimate quantitatively the efficiency of THz generation η = P rad /P b in the system considered. Figure 16 shows what fraction of the total power injected with the beams P b is converted to the radiation power P rad at each harmonic for all regimes discussed in the work. The power of the thin beam reaches the value P 1 = 3.67 GW and the power of the wide beam-P 2 = 7.33 GW. In the scheme of a single beam, efficient radiation can be generated only if the focused beam is guided by a strong magnetic field (Figure 7). In this case, the magnetic field prevents the radial dynamics of a plasma, which creates favorable conditions for the antenna mechanism at the stage of axial plasma density modulation. The maximum efficiency of beam-to-radiation power conversion reaches 1.6% and the spectrum of produced radiation is dominated by the fundamental frequency ω 0 . In the regime of counterstreaming beams, the highest efficiency is achieved for the collision of thin identical beams with the total injection power 2P 1 = 7.33 GW (Figure 16a). The peak efficiency of the second harmonic emission in this case reaches 1.4% (103 MW), 0.6% (44 MW) goes to the third harmonic and 0.3% (22 MW) goes to the forth one. In the case of beams with different focal spot sizes (Figure 16b), emission at the third harmonic with the conversion efficiency of up to η ≈ 0.5% (55 MW) dominates. Transition to plasmas with light ions such as hydrogen results in a significant decrease of emission at the third and forth harmonics and most of radiation is generated at the second harmonic due to the head-on collision of different-size waves (Figure 16c,d).

Discussion & Conclusions
Recent studies have shown that a thin beam-plasma system (with dimensions comparable to the radiation wavelength) can convert a significant fraction (up to 5-10%) of the injected beams power into radiation power at harmonics of the plasma frequency due to the mechanisms of the plasma antenna and head-on collision of plasma waves with different potential profiles. In this work, using 2D3V PIC simulations in axially symmetric geometry we have shown how efficiently these mechanisms work for specific parameters of LIA developed in the Budker Institute. The key property of this LIA is that a small enough beam emittance is achieved in such an accelerator at relatively small energies 2 MeV and relatively high currents 1.5-2 kA. Such a small emittance allows to focus an electron beam into a spot of 1 mm diameter which is comparable with the wavelength of EM waves at harmonics of the plasma frequency (lying in THz range for the typical plasma density 10 15 cm −3 ) and suited for switching on highly-efficient radiation mechanisms in plasma. High currents corresponding to the beam-to-plasma density ratio n b /n p = 10 −3 -10 −2 are needed for efficient development of the two-stream instability. The energy of beam particles should be weakly relativistic (1-2 MeV) in order to excite longitudinally propagating plasma waves capable of participating in the discussed emission processes. Moreover, a relatively low energy of beam particles allows to shorten the distance of beam-plasma interaction to the mm-scale making it possible to realize these radiation schemes in gas jets. An increase in energy will lead to an increase in the beam relaxation length (since the growth rate of the two-stream instability is inversely proportional to the relativistic factor) making it longer than the focusing region. Moreover, at high relativistic energies, the spectrum of unstable waves ceases to be dominated by longitudinally propagating modes and one should look for other efficient radiation mechanisms such as three-wave interaction discussed in [31].
PIC modeling of a single beam injection into a plasma without a magnetic field showed that the beam-plasma instability developing near the focus of the beam leads to effective beam bunching with the period of the most unstable wave (several c/ω p ), while this wave itself (being subluminal) cannot be converted into an electromagnetic one. The possibility of such a conversion opens up if the longitudinal density modulation appears in the plasma with a period comparable to the length of the dominant beam-driven wave. However, in the absence of a magnetic field, the radial displacement of plasma from the beam region under the action of the ponderomotive force of plasma oscillations turns out to be a faster process than axial modulation, which does not allow efficient use of the plasma antenna mechanism in this case. Nevertheless, the phenomenon of deep beam bunching near its focus can be used to generate THz radiation due to the Smith-Purcell mechanism [40,41] or coherent diffraction and transition radiation. A short plasma section (e.g., a supersonic gas jet) in which a periodic sequence of bunches is formed can be used for this purpose. Then these bunches passes through a hole in a metal screen or falls on a metal foil and generate radiation [42]. Switching on a strong magnetic field prevents the radial displacement of the plasma from the beam region, while the growing longitudinal modulation of the plasma density leads to a significant increase in the radiation power at the plasma frequency. The maximum efficiency of such radiation reaches 1.6%.
Another possibility of efficient THz generation is associated with the use of two colliding LIA beams. In addition to the nonlinear interaction of the most unstable modes (ω 0 , ±k 0 ), leading to the generation of radiation at the second harmonic, it becomes possible in such a system to generate superluminal oscillations at the third (3ω 0 , ±k 0 ) and the fourth (4ω 0 , ±2k 0 ) harmonics of the plasma frequency. The maximum efficiency of the second harmonic radiation that is achieved by switching on the mechanism of head-on collision of waves with locally different transverse structures, is 1.4%, which in absolute values corresponds to a power of 100 MW. The maximum amplitude of the radiation field in this case reaches E z ≈ 10 −2 , which corresponds to the value of 340 kV/cm for the plasma density n 0 ≈ 1.2 × 10 15 cm −3 . The most promising development of this scheme is the use of an external magnetic field. This will allow the beams to maintain a given transverse size at a large enough distance. In this case, the use of different-size beams will make it possible to generate radiation near the second harmonic of the plasma frequency at a greater distance and for a longer time, which will significantly increase the efficiency of the scheme.
The duration of radiation in two-beam schemes is limited by effects of ion dynamics. For the plasma density 1.25 × 10 15 cm −3 , this duration is fractions of a nanosecond. With an increase in the mass of ions (in this work, we considered only the lightest hydrogen ions and infinitely heavy ions), the duration of the radiation increases. In the case of a single beam injection into a magnetic field (Section 3.1), ion dynamics is a necessary condition for switching on the mechanism of the beam-plasma antenna. At much lower plasma densities, this mechanism can operate at much longer times (see experimental results from [24] and simulations of such experiments in [39]).
Compared to other THz sources based on the use of highly relativistic electron beams [43], plasma-based sources at lower beam energies are attractive for several reasons. The first advantage is the engineering simplicity. In addition to the particle accelerator, only the plasma section is needed, which in fact is a simple gas cell or supersonic gas jet. Second, it is possible to tune radiation frequency within a wide range by a simple changing of the plasma density. Third, in such a scheme it is possible to achieve a THz spectrum width of the order of one percent or less with a high electric field of the order of MV/cm.
The methods developed in this work for simulating the interaction of relativistic electron beams with plasma are relevant not only for studying promising sources of THz radiation. The need to investigate this kind of beam-plasma systems also arises in the development of modern X-ray radiography systems based on linear induction accelerators [44,45].

Institutional Review Board Statement: Not applicable.
Data Availability Statement: All data included in this study are available upon request by contact with the corresponding author. The code used is available at https://bitbucket.org/berendeev/ luthien_rz, accessed on 1 May 2021.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: Appendix A Figure A1. Dependence of the field E z on time for each harmonic of emission at points on the boundary of the computational domain. To construct these curves for each harmonic, the coordinate of the maximum field (z max , t max ) was found. Then the time dependence E z (z = z max , t) was plotted. (a) same beams with σ 1 = σ 2 = 3.325c/ω p ; (b) different beams with σ 1 = 3.325c/ω p and σ 2 = 6.65c/ω p ; (c,d)-same with (a,b), but in presence of hydrogen ions.