© The Authors 2025

1. Introduction

The interaction of electron beams with plasmas is one of the fundamental processes of plasma physics. Electron beams have a wide range of applications in laboratory plasma systems; for example, for plasma heating (Burdakov et al. 2007), emission generation (Arzhannikov et al. 2020), ionisation of neutral gas, and acceleration of electrons (Soldatkina et al. 2022). Recent experiments (Schönlein et al. 2016) on irradiation of metal wires demonstrate the possibility of creating electron beams inside such targets and heating them to temperature of dozens of electron volts. Such experiments should allow one to get insights into warm density matter evolution, which is relevant not only for research in the area of inertial fusion, but also for studies of dense stars. Being accelerated during various high-energy events in the solar atmosphere (Aschwanden 2002; Klein et al. 2005; Jebaraj et al. 2023; Khotyaintsev et al. 2019), electron fluxes cause a wide variety of electromagnetic (EM) emission processes (Ginzburg & Zheleznyakov 1958; Melrose 1980; Robinson & Cairns 2000; Reid & Ratcliffe 2014) and also possibly have an impact on artificial satellites (Hapgood 2016).

Great progress in the theoretical description of the beam-plasma interactions has been made in over a century of plasma research. However, an analytical description of these processes is strictly limited to different simplified cases, while many practically important regimes can be investigated only through numerical simulations, mainly using the particle-in-cell (PIC) approach (Dawson 1983). With rapid growth in computing power, it has become possible to simulate vast plasma systems using different numerical schemes with high resolution on large timescales. However, for real systems to be adequately described by means of simulations, it is also necessary to choose the correct model of simulation box. By now, the majority of PIC simulations of solar beam-plasma systems (e.g. Thurgood & Tsiklauri 2016; Henri et al. 2019; Zhang et al. 2022; Krafft & Savoini 2023; Li et al. 2024) are still performed in simplified geometry with periodic boundary conditions, which correspond to the infinite plasma case.

The model of infinite plasma is actually assumed in the majority of analytical estimates of different plasma processes; for example, in the classic theory of beam-plasma instability (Akhiezer & Fainberg 1949; Bohm & Gross 1949). That makes using a similar model in numerical simulations the best choice for the verification of many theoretical concepts. The model of infinite plasma can be realised for any spatial geometry via implementation of the periodic boundary conditions for fields and particles. Figure 1 demonstrates the scheme of the simulation box in 2D geometry for infinite plasma. Due to the periodic boundary conditions, the particles actually move on the surface of a torus, while the spectrum of EM fields become restricted by the spatial size of the simulation box. Typical initial conditions for beam-plasma systems in this model are as follows. Particles of the plasma species (ions and electrons) are distributed uniformly over the simulation box. In the same manner, electrons of the beam are placed in all space with the initial temperature distribution and directed velocity. Meanwhile, electrons of plasma have not only a temperature spread, but also some mean velocity in the opposite direction to the beam in order to construct a compensating current. After such an initiation, we can observe the evolution of the system in time. The beam drives various instabilities (Bret et al. 2008). Simulations usually correctly reproduce predictions of the theories on the linear stage. Then, one can observe the vast majority of non-linear processes. However, such initial conditions do not actually correspond to real beam-plasma systems. Whether in experimental or natural conditions, we always have a definite source for the beam particles, which propagate through plasma. For example, solar radio bursts can be generated by particle fluxes, which originate far from the radio emission source. On the one hand, this means that these beams will have a complicated distribution function due to their interaction with plasma before the emission site. On the other hand, the emission site is a plasma region of finite size with an inflow of beam particles from one side and an outflow from the other. These particles interact with the plasma, build up different waves, and leave the area of interest. The duration of typical radio bursts is much longer than the plasma period in the corresponding area. So, from the point of view of now-available beam-plasma simulations with a length of up to 10 000 inverse plasma frequencies $\left({\omega }_{p0}^{-1}\right)$, such a system will look like the constant injection of an electron beam in a finite size plasma column with open boundaries in the direction of beam propagation (Fig. 2). In this model, we correctly reproduce the real situation of a constant inflow of new particles into the system, while in the model of infinite plasma the electron beam has only a limited amount of initial energy and throughout the entire non-linear stage realistic energy income is absent. Another major difference between these models is the possibility that EM waves leave the system in the case of open boundary conditions, while in infinite plasma these waves remain in the simulation box and continue to interact with particles and other waves.

Fig. 1. Schematic representation of a simulation box for infinite plasma. The cyan arrow indicates the path of one particle, which passes through all periodic boundaries.

Fig. 2. Simulation box with open boundaries in the longitudinal direction (x axis) and periodic boundaries in the transverse one (y axis). Beam particles are injected into the system through the left edge. They cause a return current of plasma electrons and leave the system through the right border. The flux of all of the plasma particles is maintained in a self-consistent manner on both sides.

The necessity to use a simulation box with open boundary conditions and a continuously injected beam in proper simulations of processes in space plasma was clear dozens of years ago (Büchner et al. 2003), but for a long time such simulations were limited to different simplified cases. For example, in work by Mandrake et al. (2000) a 2.5D electrostatic code was used, and a 1D spatial geometry for dense beams was used by Lizunov et al. (2002). The simulation of low dense electron beams was hampered by the much larger amount of computational resources required for the model with beam injection. In this case, we have to use a simulation box long enough to fit the full relaxation length, L_R ∝ v_b/Γ, where v_b is the mean beam velocity, and Γ the growth rate of the most unstable plasma mode. The lower the Γ, the longer the system we have to use, and therefore the more computational resources are needed. This is the case for the low-density beam or the hot temperature of the electrons in the system. On the other hand, when we use the periodic boundary conditions to simulate beam-plasma interactions, a simulation box of even one plasma wavelength is enough to observe beam-plasma instability and the formation of plasma waves. A further increase in the length of the simulation box in the direction of beam propagation will only increase the spectrum of waves for which existence is possible in the system. This is of course important if we are interested in investigating non-linear wave transformations in long-term dynamics.

In recent years, we have observed a rapid growth in available computation power, which is opening up huge possibilities for simulating plasma processes. In particular, at the Budker Institute of Nuclear Physics, we were able to switch to PIC simulations using the model with open boundaries and continuously injected beams after developing a code for then-new general-purpose graphical processing units (GPGPUs) (Lindholm et al. 2008), which turned out to be a perfect solution for highly parallelised tasks like plasma simulations. Over the next decade, we investigated a lot of beam plasma systems in 2D spatial geometry, mainly in the context of plasma emission tasks (Annenkov et al. 2016a, 2018, 2019; Annenkov et al. 2020, 2024), electron acceleration in an open magnetic trap (Timofeev et al. 2022), and examinations of the influence of plasma density gradients (Annenkov & Volchok 2023). The possibility of simulating large-scale gradients is one more fundamental advantage of the open-boundary model over the infinite one. Other recent examples of successful usage of the model with open boundaries include the work by Nishikawa et al. (2023), in which relativistic jet injection into a plasma was simulated using the TRISTAN code (Spitkovsky 2005) on times up to 900 ${\omega }_{p0}^{-1}$ in full 3D spatial geometry, and the work by Kumar et al. (2022) in which EM emission of colliding electron beams was simulated by the EPOCH code (Arber et al. 2015) in 2D geometry. However, all of the studies mentioned above deal with relatively dense and mostly high-energy electron beams.

This work pursues several objectives. The first one is to investigate the relaxation process of a weak electron beam in a solar wind plasma using PIC simulations in a realistic model with open boundary conditions. For this, the numerical model and approaches being approbated in the context of laboratory experiments were used. The second objective is to demonstrate and discuss main differences of the obtained results with the simulations using infinite model. For this, the system parameters were chosen according to the parameters actively used for PIC simulations of electron beams in solar wind plasmas in the infinite plasma approximation by Krafft & Savoini (2021, 2022, 2023, 2024). The final objective is to briefly discuss other possible approaches to numerical simulations of beam-plasma systems in the context of solar-related problems.

2. Numerical model and parameters

To simulate a beam-plasma system, the 2D3V Cartesian PIC code with a standard numerical scheme was used. It is based on the Yee (1966) solver of Maxwell equations for EM fields, the Boris (1970) scheme for solving the equation of motion for collisionless macro-particles with a parabolic form factor, and the charge-conserving Esirkepov (2001) scheme for calculations of currents.

The proper open boundary conditions should support several features: they have to first be able to remove from the simulation box outgoing particles and EM waves; second, imitate the presence of real particles and fields beyond the edge of the simulation plasma slab; and third, realise the correct inflow of particles into the simulation area. The former feature is crucial not only for sustaining the thermal flow of plasma particles, but also to maintain the correct return current, which should compensate for the current of the beam. Moreover, this compensation should be realised in a self-consistent manner, because beam distribution and mean velocity can be changed dramatically and unpredictably during propagation in plasma. In the code in use, the self-consistent plasma buffers are realised. They fit these requirements and are described in Annenkov et al. (2018).

The simulations were performed using four Nvidia A100 GPGPU. The scheme of the simulation box is presented in Fig. 2. Open boundary conditions were imposed in the longitudinal direction (x axis).

The injected electron beam has a relative density of n_b/n₀ = 5 ⋅ 10⁻⁴, where n₀ is the density of the unperturbed plasma. This beam density is an order or two higher than the typical values measured in the solar wind. However, when such a density is used in simulations with periodic boundaries the resulting level of Langmuir turbulence becomes comparable to in situ measurements. This is caused by the limited amount of energy in the simulated beam and absence of fresh particles coming into the system, as is the case in real beam-plasma systems. The long-term non-linear dynamic of plasma waves in an infinite plasma system should also differ, because plasma waves do not have the possibility to interact with fresh beam particles and their associated oscillation modes. The mean velocity of the beam is v_b = 0.25c (≈17 keV), where c is the speed of light. All particles in the system initially have a Maxwell distribution:

$$\begin{array}{c}\hfill f(p)exp[-\frac{{(p-p_0)}^2}{\mathrm{\Delta }{p}^2}],\end{array}$$(1)

where p is the momentum of the particles in units of m_ec, m_e is the rest mass of an electron (511 keV), p₀ is the mean momentum, and $\mathrm{\Delta }p=\sqrt{2\widehat{m}T/m_e{c}^2}$, where $\widehat{m}=m/m_e$. For electrons in the plasma and the beam, the initial temperature is T_e = 320 eV, which corresponds to a thermal velocity of v_T ≈ 0.1v_b. For ions in the plasma, the temperature T_i = 0.1T_e was used. Such a ratio for the temperatures of plasma components will prevent damping of ion acoustic waves. In contrast to the infinite plasma model, it is not necessary to artificially compensate for the beam current by giving the plasma electrons an average velocity in the direction opposite to the beam direction. The external magnetic field was set to zero. The calculation of the full unstable spectrum in this beam-plasma system using the DispLib code (Annenkov et al. 2021) demonstrates (Fig. 3) the domination of purely longitudinal modes. The maximum growth rate, Γ/ω_p0 ≈ 1.4 ⋅ 10⁻², corresponds to the wave number, k_∥ = 4.29 ω_p0/c, where ∥ indicates a direction along the beam axis.

Fig. 3. a) Growth rate map for the beam-plasma instability, Γ(k∥, k⊥). The green line, k⊥ = k⊥(k∥), marks the local maximal growth rate. b) Growth rate for longitudinal modes, Γ(k∥, 0). c) Growth rate on the line of the maximum, Γ(k⊥).

The relaxation length of electron beams with low relative density imposes a requirement to use a relatively long simulation box. Moreover, it is also imposes a restriction on the minimal simulation time. It has to be long enough to fit the whole linear stage of the beam-plasma instability and give the ability to study the non-linear stage. Fulfilling both these conditions requires additional computation power. Since the most unstable beam-plasma modes are purely longitudinal, it is possible to restrict first studies to 1D regime. So, in this work, the simulations are limited to the quasi-1D case. The simulation box has N_y = 4 cells and periodic boundary conditions in the transverse direction. Such a choice allows us to study the main features of the beam-plasma instability with high resolution in space and time on large spatial scales and timescales, but prohibits the investigation of non-linear process, which involve transverse modes. Due to these restrictions, the present simulations do not allow for deep studies of the EM emission generation, but make possible an investigation of the evolution of the distribution functions of a plasma and a beam, as well as the growth and evolution of longitudinal plasma oscillations.

The longitudinal size of the simulation box is N_x = 3600 cells, which corresponds to a length of L_x = 72 c/ω_p0 = 2880 λ_D, where λ_D = v_T/ω_p0 is the Debye length. It should be noted that for simulations with open boundary conditions one should choose the size of the area with an extra space, because processes near the injector (e.g. the formation of strong density gradients Annenkov et al. 2019) can shift the area of intensive beam-plasma interactions further from the injector. For each particle species, N_p = 2025 macro-particles per cell were used. The size of cells is Δx = Δy = 0.02 c/ω_p0 and the time step is Δt = 0.01 ${\omega }_{p0}^{-1}$. The total simulation time is 10 000 ${\omega }_{p0}^{-1}$. Such small spatial and time steps enable all plasma oscillations to be investigated with high precision, while the long simulation time makes it possible to catch the process caused by the ion dynamic. At the initial time step, all particles of the plasma electrons and ions with the same densities, n_e = n_i = n₀, are distributed across the simulation box with the initial thermal spread.

Later in the text, all quantities are presented in a dimensionless form. Plasma and beam particle densities were calculated in units of n₀; all frequencies in ${\omega }_{p0}=\sqrt{4\pi n_0{e}^2/m_e}$; wave numbers in ω_p0/c; lengths in c/ω_p0; time in ${\omega }_{p0}^{-1}$; EM fields measured in units of m_ecω_p0/c; and particle velocities in the speed of light, c. Since the conversion was based on the initial density of plasma, n₀, all of the simulation results can be applied to different real systems with various plasma densities. For example, for coronal plasma with n₀ = 10¹⁰ cm⁻³, the unit of length is 1 c/ω_p0 ≈ 5.3 cm and the unit of time is $1{\omega }_{p0}^{-1}\approx 0.17$ ns, while for solar wind plasma with n₀ = 10⁵ cm⁻³, the unit of length is 1 c/ω_p0 ≈ 17 m and the unit of time is $1{\omega }_{p0}^{-1}\approx 55$ ns. For convenience, the conversion factors from dimensionless to dimensional units for various densities are available in Fig. A.1 in the appendix. However, the results of the current simulation are applicable only to systems with temperatures similar to ones we have chosen.

3. Simulation results

First, we shall discuss the process qualitatively. In Fig. 4, snapshots of the simulation results are presented. At the initial stage (Fig. 4 top), the electron beam undergoes two-stream instability and builds up plasma waves with the phase velocity of the beam particles. The spatially localised region with high frequency field oscillations is formed in the vicinity of the beam injection site (Fig. 5). The presence of continuously injected fresh beam electrons leads to the rise in the wave amplitude. It turns out that this amplitude along with the other system parameters is enough to make plasma waves unstable against the modulation instability (Vedenov & Rudakov 1964; Zakharov 1972; Thornhill & ter Haar 1978). Later time moments demonstrate the formation of ion density wells and the capturing of the plasma waves inside them. Such trapped oscillations fall out from the resonance with the electron beam and gradually expend their energy on well deepening. After they have burned out, regions with high density gradients are formed. Such gradients, in turn, prevent the development of beam-plasma instabilities. After some time, these density perturbations can be smoothed out by thermal ion dynamics, and conditions for instability development return. This process is described in more detail in the work Annenkov et al. (2019), in which it is also demonstrated that the formation of such a longitudinal modulation provides the necessary conditions for the beam-plasma antenna mechanism (Annenkov et al. 2016b; Timofeev et al. 2016) and the generation of EM emission near the plasma frequency and its second harmonic. Moreover, as can be seen from phase portraits in Fig. 4, the observed waves provide appropriate conditions for acceleration of the beam electrons to supra-thermal velocities.

Fig. 4. Simulation results at several moments in time. Top: Density of the ions (red line), different phases of one plasma oscillation (blue lines). Bottom: Phase space, f(vx, x), of the beam. The corresponding video is available online.

Fig. 5. Top: Evolution of the field, Ex, maximum in space and time. Bottom: Time dependence of the field, Ex, maximum (blue) and its position (brown). The horizontal blue line marks the average value of the longitudinal electric field, Ex, in the plasma wave during the non-linear stage.

Now we shall look at the results in more detail. In Fig. 6, the growth rate of the beam instability in simulations is shown. It was calculated using the energy, H(t), of the electric field in the whole area as follows:

Fig. 6. Dependence of the instability growth rate, Γ, on time. The horizontal solid black line indicates a prediction of the linear theory.

$$\begin{array}{c}\hfill {\mathrm{\Gamma }}^{\mathrm{PIC}}(t)=\frac{H\prime (t)}{2H(t)}·\end{array}$$

This graph allows us to define the boundary between two stages of the system evolution. At time moment t₀ = 900 ${\omega }_{p0}^{-1}$, the growth rate dependence shift from rapid changes to a flat-like regime. Before t₀, the regime of beam-plasma instability is almost linear and after it is non-linear. As we shall see further, this boundary also appears on the time dependencies of other integral system parameters. One more important difference between simulations in the infinite plasma model and the one with open boundaries should also be noted. In the first model, we initially placed in the simulation box all beam particles with the initial distribution function. Therefore, all beam electrons participate in the instability development at once, in good agreement with the predictions of the linear theory (Fig. 3). In the second model, the realistic situation when the fastest particles from the initial velocity distribution function (VDF) reach new plasma layers takes place. In this situation, the head of the beam build up plasma waves not in agreement with the predictions of linear theory for the whole initial distribution and the core of the beam enters in plasma layers with already existent plasma waves. This makes even the linear stage of the beam-plasma instability more complicated.

In Fig. 7, the distribution function, f(v_x), of all electrons (plasma and beam) at different time moments and in various areas of the simulation box are presented. Plasma electrons are heated, but their distribution remains Maxwellian. At the end of the linear stage (t ⋅ ω_p0 = 900), all distributions have areas with a positive slope of ∂f/∂v_x, and therefore continue to be unstable. At later moments, this is only the case for the first area, which is located near the beam injector site. However, at later moments, we can see significant acceleration of some beam particles to velocities far above the thermal speed. At the beginning, 20% of the beam energy is contained in supra-thermal particles (Fig. 8). During the simulation, that percentage increases to 50%. The average value over the entire non-linear stage in regions 1 and 2 turns out to be near 35%. This result exceeds by 5% the same value from the simulation in the infinite model, even for the initially inhomogeneous plasma density. It is also worth noting that in the model with open boundaries we obtained such an acceleration during the short time in which the beam particle was present in the system before leaving it, while in the model with periodic boundaries each particle remains in the system for the entire simulation time and has many chances to be accelerated. Another parameter that can be compared with the simulations by Krafft & Savoini (2023) is the fraction of beam energy that was transferred to waves. In the case of infinite plasma with an inhomogeneous initial density, it turns out to be ≈10%. The same value appears to be in the present simulations (Fig. 9). However, 10% is only the average value. We also can see the number of time moments when beam particles leaving the system significantly increase their energy in comparison with the time of injection.

Fig. 7. Velocity distribution functions, f(vx), of all electrons at different moments of time. The dashed gray line indicates the initial distributions. Different lines (0, 1, 2) corresponding to particles in the areas are indicated as 0, 1 and 2 in Figure 4. The corresponding video is available online.

Fig. 8. Fraction of the beam energy that is being carried out by supra-thermal particles (with v > vT) for each diagnostic area. The dashed black line indicates the average level over the time from t = 900 ${\omega }_{p0}^{-1}$ to the end of simulation for the second area.

Fig. 9. Power of the beam at the end of the system normalised on its power at the injection moment. The dashed black line indicates the average level over the time from t = 900 ${\omega }_{p0}^{-1}$ to the end of the simulation.

Turning to the spectrum of oscillations in the system, according to the linear theory (Fig. 3), the most unstable mode is one with k_∥ ≈ 4.3. In simulations, it is close to this value and appears to be near k_∥ ≈ 4.5 at the initial stage (Fig. 10). On the non-linear stage spectrum, there are various numbers of significant modes with wave numbers up to k_∥ = 7. In Fig. 12, the (ω, k_∥) spectrum over the entire simulation time is presented. Aside from the plasma oscillations near the plasma frequency, there is a noticeable amount for the second harmonic and evidence for the third one, but higher modes are not visible. According to the (ω, x) spectrum (Fig. 11), all high modes of plasma oscillations are localised in area 1, mainly near x = 10 c/ω_p0, where the amplitude of plasma waves are the highest and non-linear processes are the most intensive (see Fig. 5). In addition to plasma oscillations, we can see ion acoustic waves with a wide spectrum as well as a small number of EM waves in the B_z spectrum.

Fig. 10. Spectrum of the longitudinal electric field, Ex, over the simulation time. The horizontal pink line indicates a prediction of the linear theory. The vertical white line marks the time moment $t·{\omega }_{p0}^{-1}=900$.

Fig. 11. Spectrum of plasma oscillations over the system length.

Fig. 12. Spectrum Bz(ω, k∥) (top right) and Ex(ω, k∥) (others) over the entire simulation time and along the entire length of the system. The vertical red lines correspond to ω = |k|⋅vb.

4. Discussion and conclusion

The present simulation in the open boundary model with a continuously injected beam shows a higher number of supra-thermal particles appearing than in the model of infinite plasma. During the propagation through the plasma slab, beam particles lose on average 10% of their energy. We also can see some time moments when beam particles, on the contrary, increase their average energy up to 125% of the initial value. However, the most important results are the following. For the considered relative beam density, n_b/n₀ = 5 ⋅ 10⁻⁴, we see that the relaxation length appears to be only several dozens of c/ω_p0. Within this length, the level of beam-plasma instability is high enough to be followed by modulation instability and plasma wave capturing in the density wells. After passing through this distance, the electron beam distribution function appears to be almost flat-topped. Such a beam becomes stable against two-stream instability and passes forward without significant changes.

These results show the importance of using a realistic numerical model with open boundaries and a continuously injected beam when investigating beam-plasma systems in space conditions. However, to accurately reproduce the actual level of plasma turbulence in the solar wind, the relative density of the beam should be reduced by one or two orders of magnitude. Although the use of this model is more computationally demanding than the infinite plasma one, modern computational systems are sufficient for high-resolution studies on large spatial and temporal scales.

In addition, several more approaches for realistic simulations of beam-plasma systems can be suggested. We see that after some time the system reaches some sort of steady state regime. So, instead of continuous injection of a beam, we can consider a case with a limited injection time. This would allow us to study the evolution of the waves in the system that was built up by the beam, but without beam particles. Such a problem statement is relevant for the interpretation of real beam-plasma systems, but it is impossible in the infinite plasma approximation. In the present work, the injected beam has a Maxwellian initial velocity distribution. But, as we can see from simulations, the beam distribution rapidly becomes flat-topped. So, if our goal, for example, is to investigate solar radio bursts at a significant distance from the acceleration site, then we can simulate the injection of a beam with a flat-topped VDF. This allows us to investigate the influence of realistic regular and random gradients of plasma density and the external magnetic field on beam-plasma instability, which is necessary for the appearance of plasma emission mechanisms. On the other hand, if we are mainly interested in the evolution of a beam’s VDF over a long propagation distance, then it seems appropriate to use the model of a moving frame. In such a model, we can start from an electron beam with the initial distribution on a generation site (e.g. near the Sun) and investigate its evolution during the propagation through the solar plasma, taking into account all variations of external parameters (background plasma density and magnetic field). Nowadays, this model is mainly used for plasma accelerators (Vay 2020) with a frame length of several plasma lengths. As we can see from the present simulations, for solar plasma systems with low density beams, this length should be significantly (probably hundreds of times) longer. This is necessary to develop main beam-plasma instability at the time of beam propagation. Because of the significant demand for computational resources, most likely only 1D simulations using this model are possible nowadays. However, the result of such simulations should provide an unprecedented possibility for a direct comparison with the satellite measurements of beam VDFs.

Movies

Movie 1 Access here

Movie 2 Access here

Data availability

Movies associated to Figs. 4 and 7 are available at https://www.aanda.org

Acknowledgments

The work was supported by the Ministry of Science and Higher Education of the Russian Federation. The author would like to thank Irkutsk Supercomputer Center of SB RAS for providing the access to HPC-cluster “Akademik V.M. Matrosov”.

References

Appendix A: The conversion factors

The conversion factors from dimensionless to dimensional units for various densities are presented in Fig. A.1. By choosing the background plasma density we determine unambiguously all parameters in the simulations. This makes it possible to apply the results obtained in the simulations to various physical systems for which the collisionless plasma approximation is applicable.

Fig. A.1. Correspondence between dimensional and dimensionless quantities on a logarithmic scale. a) The unit of time ${\omega }_{p0}^{-1}$. b) The unit of length c/ωp0. c) The linear frequency of radiation at the second harmonic of the plasma frequency. d) The characteristic amplitude of the electric field of plasma waves.

All Figures

	Fig. 1. Schematic representation of a simulation box for infinite plasma. The cyan arrow indicates the path of one particle, which passes through all periodic boundaries.
In the text

	Fig. 2. Simulation box with open boundaries in the longitudinal direction (x axis) and periodic boundaries in the transverse one (y axis). Beam particles are injected into the system through the left edge. They cause a return current of plasma electrons and leave the system through the right border. The flux of all of the plasma particles is maintained in a self-consistent manner on both sides.
In the text

	Fig. 3. a) Growth rate map for the beam-plasma instability, Γ(k∥, k⊥). The green line, k⊥ = k⊥(k∥), marks the local maximal growth rate. b) Growth rate for longitudinal modes, Γ(k∥, 0). c) Growth rate on the line of the maximum, Γ(k⊥).
In the text

	Fig. 4. Simulation results at several moments in time. Top: Density of the ions (red line), different phases of one plasma oscillation (blue lines). Bottom: Phase space, f(vx, x), of the beam. The corresponding video is available online.
In the text

	Fig. 5. Top: Evolution of the field, Ex, maximum in space and time. Bottom: Time dependence of the field, Ex, maximum (blue) and its position (brown). The horizontal blue line marks the average value of the longitudinal electric field, Ex, in the plasma wave during the non-linear stage.
In the text

	Fig. 6. Dependence of the instability growth rate, Γ, on time. The horizontal solid black line indicates a prediction of the linear theory.
In the text

	Fig. 7. Velocity distribution functions, f(vx), of all electrons at different moments of time. The dashed gray line indicates the initial distributions. Different lines (0, 1, 2) corresponding to particles in the areas are indicated as 0, 1 and 2 in Figure 4. The corresponding video is available online.
In the text

	Fig. 8. Fraction of the beam energy that is being carried out by supra-thermal particles (with v > vT) for each diagnostic area. The dashed black line indicates the average level over the time from t = 900 ${\omega }_{p0}^{-1}$ to the end of simulation for the second area.
In the text

	Fig. 9. Power of the beam at the end of the system normalised on its power at the injection moment. The dashed black line indicates the average level over the time from t = 900 ${\omega }_{p0}^{-1}$ to the end of the simulation.
In the text

	Fig. 10. Spectrum of the longitudinal electric field, Ex, over the simulation time. The horizontal pink line indicates a prediction of the linear theory. The vertical white line marks the time moment $t·{\omega }_{p0}^{-1}=900$.
In the text

	Fig. 11. Spectrum of plasma oscillations over the system length.
In the text

	Fig. 12. Spectrum Bz(ω, k∥) (top right) and Ex(ω, k∥) (others) over the entire simulation time and along the entire length of the system. The vertical red lines correspond to ω = |k|⋅vb.
In the text

	Fig. A.1. Correspondence between dimensional and dimensionless quantities on a logarithmic scale. a) The unit of time ${\omega }_{p0}^{-1}$. b) The unit of length c/ωp0. c) The linear frequency of radiation at the second harmonic of the plasma frequency. d) The characteristic amplitude of the electric field of plasma waves.
In the text