On the ultrarelativistic two-stream instability, electrostatic turbulence and Brownian motion

Experimental evidence indicates that bulk plasma flow at ultrarelativistic speeds is common in astrophysical settings, e.g. the collimated jets of active galactic nuclei and gamma ray bursts. The low-plasma density of such flows implies their collisionless relaxation by means of wave-particle interactions. Such processes are not well understood in the ultrarelativistic regime. The thermalization of two interpenetrating equally dense electron–proton (e−p) beams in the absence of a magnetic field is examined here by means of 1.5D particle-in-cell simulations. The relative beam speeds correspond to Lorentz factors in the range 200–1000. The constraint to one spatial simulation dimension, which is aligned with the beam velocity vectors, implies that only the two-stream (TS) instability and the Weibel-type instability can grow, while filamentation instabilities are excluded. With this constraint and for our plasma parameters, the TS instability dominates. The electrostatic waves grow, saturate by the trapping of electrons, and collapse. The interaction of the electrons with the electric fields after the wave collapse represents a relativistic Wiener process. In response, the electrons are rapidly thermalized. The final electron distribution can be interpreted as a relativistic Maxwellian distribution with a high-energy tail arising from ultrarelativistic phase space holes.

These flows cannot relax via binary collisions due to the long mean-free path of the plasma particles, and instabilities develop like the two-stream (TS) instability [2]- [4] and the Weibel instability [5]- [8]. This relaxation is probably a source of cosmic rays and the intense release of electromagnetic radiation by GRBs. The TS instability, a microinstability driven by non-thermal particle distributions in plasmas, is thought to be an important relaxation mechanism for relativistic flows. It can cause the growth of waves that have a wide range of angles between their wavevectors k and the plasma flow direction. The electrostatic waves that are driven by non-relativistic or mildly relativistic beams of charged particles, have a k that is practically aligned with the beam flow direction. Such beams and waves are observed, for example, close to the Earth bow shock [9]. The dynamics of such systems in the absence of a magnetic field can be well approximated by a one-dimensional model, since only wavevectors parallel to the beam flow direction have to be represented. One-dimensional particle-in-cell (PIC) simulations [10] have examined the TS instability for non-relativistic and for mildly relativistic beam speeds in various settings [11]- [17].
As the beam speed is increased to relativistic values, the wavevector of the fastest growing mode is rotated away from the plasma flow direction, as has been shown for cold electron beams in an unmagnetized plasma [3,4,18]. The modulation of the plasma by this wave implies that the plasma rearranges itself into filaments with currents that vary in the plane orthogonal to the plasma flow direction and the process becomes three-dimensional. Its modelling is a challenge for PIC simulations that can be met if one considers leptonic plasma [19]- [22], since then the dynamics is determined only by the electrons and positrons. The ions in GRB jet plasma and in the medium into which the jets expand introduce a second timescale, and the electron and ion dynamics in multiple dimensions currently cannot be resolved if the full particle mass ratio is retained. Recent simulations [23,24] thus examine systems in which the proton to electron mass ratio is reduced to the range 20-40. Alternatively, the three-dimensional dynamics can 3 Institute of Physics ⌽ DEUTSCHE PHYSIKALISCHE GESELLSCHAFT be separated into lower dimensional case studies that can be examined at a high-simulation resolution and at the full mass ratio, each of which may reveal important aspects of relaxation mechanisms for relativistic plasma beams.
Here, we model a one-dimensional system that is aligned with the plasma flow direction. In such a system, the electrostatic TS instability competes with the Weibel instability that leads to electromagnetic waves with k parallel to the plasma flow direction. The Weibel instability has been introduced for non-relativistic beam speeds in [5] and examined for a wide range of physical settings, e.g. for cosmological magnetic field generation and magnetic field generation at GRBs [6,8]. We can, however, show that the Weibel instability does not play a significant role in our simulations.
We examine the development of the electrostatic TS instability in a spatially homogeneous plasma. By modelling only one spatial dimension we can retain the full proton to electron mass ratio and assess quantitatively the involvement of the ions in the development of the instability, in the mechanisms by which ultra-relativistic plasma flows thermalize and in the resulting electron distribution. The plasma flow at GRBs can probably be considered homogeneous only over short spatio-temporal scales, if this is possible at all. Any inhomogeneous plasma flow through an ambient medium will yield interpenetrating streams of plasma at varying density ratios. The largest growth rate of the TS instability is obtained for equal beam densities and it should have a larger probability to develop compared to the case of large beam density ratios.
We thus consider here two counter propagating electron-proton beams that have equal densities. We expand the previous work [25] to cover a wider range of flow speeds v b , ranging from 200 γ(v b ) 10 3 . The growing electrostatic waves have a mean wavelength close to that expected for the TS instability, but they grow over a wide range of wavenumbers. The wave initially saturates by trapping electrons [26,27], after which the wave modes rapidly collapse. The resulting electron distributions are well approximated by the relativistic Maxwell distribution [28] at low energies, which is identical for all beam speeds except for a linear scaling with the beam Lorentz factor. The likely origin of this thermal distribution is the potentially relativistic Brownian motion [29] of the electrons in the electrostatic wave field. Thermal bulk distributions together with non-thermal particle populations are thought to exist at GRBs [30,31]. We observe such a high-energy tail in the electron distribution that reaches, for increasing beam speeds, to larger peak energies. For γ(v b ) = 10 3 , this tail has its cut-off at about 10 10 eV at the simulation's end and in the reference frame of a beam. We show that this high-energy tail is due to electron jets [32] or phase space holes [33,34] that accelerate out of the thermal component.

Equations and initial plasma consideration
We consider here two electron-proton beams with equal densities. The non-relativistic electron plasma frequency is ω p,e = (n e e 2 / 0 m e ) 1/2 , where n e , e, 0 and m e are the electron number density, the magnitude of the electron charge, the dielectric constant and the electron mass. The non-relativistic proton plasma frequency is ω p,p = (n p e 2 / 0 m p ) 1/2 and the proton mass is m p . We consider charge neutral beams by setting n e = n p . These beams flow at ultrarelativistic speeds relative to each other. We consider only one-spatial direction, the plasma flow direction, which we align with the x-direction or x. Beam These equations, together with ρ x = i q i ∞ −∞ f i (x, p, t) dp and j = i q i ∞ −∞ p f i (x, p, t) dp, determine the time evolution of the initial phase space distribution. We define it in the laboratory frame, which we will use throughout this work, and in which the total plasma momentum vanishes It consists of four shifted Maxwell-Boltzmann-Jüttner distributions [28], or relativistic Maxwellian distributions. The symbol S i represents a multiplication with +1 for species 1 and 2 and a multiplication by −1 for species 3 and 4. We define the temperature of species i as T i . The normalization constants C i give the plasma frequencies ω p,e or ω p,p in the rest frame of the respective species. The beam speed V b and the Lorentz factor γ(V b ) in the laboratory frame are connected to the quantities in the beam frame by

Linear dispersion relation
The linear plasma dispersion properties of this system are determined by the dispersion relations discussed in [35]. By considering only one spatial direction, we exclude the filamentation instability, which is probably dominant for our chosen plasma parameters [4]. The filamentation instability requires the resolution of waves with wavevectors other than k x, as discussed in [3,4]. Since, apart from noise, there is initially no magnetic field present in our model, only two modes grow. These are the electrostatic TS instability and the electromagnetic Weibel instability. Their large mass implies that the protons do not contribute initially to the growth of waves. The cold TS instability for two equally dense electron beams, that determine the initial dynamics of the system and that move at the non-relativistic relative speed V d , is described in the beam frame by the dispersion relation [36] 5 The solution of (7) gives the most unstable wave frequencyω u /ω p,e = √ 3/2, the corresponding growth rateω i = ω p,e /2, and the wavenumberk u = √ 3 ω p,e /V d . The phase speed of the wave isω u /k u = V d /2. The Lorentz factor of the wave differs strongly from the beam Lorentz factor, and we can thus not introduce the same Lorentz factor for the beam and the wave into (7) as in [2]. We solve the dispersion relation (7) in the reference frame of the electrostatic wave. We transform this equation into the rest frame of the most unstable wave and define In this frame, the real part of the most unstable frequency isω u = 0 andk u = √ 3 ω p,e /2V b . Since we have applied a Galilei transform, the growth rate and the most unstable wavenumber are unchanged, i.e.ω i =ω i andk u =k u . We introduce relativistic beam speeds in this frame. Let both beams move at the speed V b that corresponds to the Lorentz factor γ(V b ) and equation (4) takes the form We substitute p,e = ω p,e /γ(V b ) and obtain from equation (9) 1 − 2 p,e which is equivalent to equation (8). Solving equation (10) gives the maximum growth rate in the laboratory frame and the wavelength is defined as The second type of waves we have to consider in our plasma setup, the electromagnetic waves excited by the ultrarelativistic Weibel instability, is considerably more difficult to describe in terms of the growth rate and wavenumber spectrum. It has previously been examined in the framework of a water bag model, where the electrons were distributed over a wide range of momenta in the x-direction [6,7]. This case does not apply here. More appropriate is the work in [3], where non-relativistic cold electron beams have been considered that move through a warm plasma, resulting in the wave growth at all k along the x-direction. The growth rate of the instability increases for increasing k. We may thus expect the growth of electromagnetic waves at large k, if they grow at all.

Simulation setup
The PIC code [10] solves the Maxwell equations (3) and (4) x and the momentum of particle j by Unless stated otherwise, we resolve the x-direction by N g = 5 × 10 4 simulation cells, each with the fixed length x. Initially, the plasma is spatially homogeneous and we set all fields to zero, and the waves thus grow out of noise [37]. We employ periodic boundary conditions that give us the fixed minimum wavenumber k min = 2π/N g x. As shown above, we obtain a scaling Both beams move at the speed modulus V b in opposite directions and the relative Lorentz We vary the beam Lorentz factor and we model the cases in table 1. The normalized thermal speed of the electrons is v th,e = (k B T e /m e ) 1/2 = 30 −1 c, and that of the protons is v th,p = (m e /m p ) 1/2 v th,e . We set the simulation timestep t c such that x/ t c ≈ 2.8c and π/ω p,e t c ≈ 170. The simulation cell size and time step in simulation (f) are x/2 and t c /2, respectively. The rate at which we sample the electric field data is ω s = 2.6 ω p,e , and the sampling timestep is t s . For simulations (a)-(e) the plasma species 1 and 3, the electrons, are represented by N e = 2.45 × 10 7 simulation particles each and the species 2 and 4, the protons, by N p = 8 × 10 6 simulation particles respectively. Simulation (f) employs three times those particle numbers.

Simulation results
In what follows, we discuss our simulation results in the laboratory frame, unless stated otherwise.

Wave spectra
The TS instability gives rise to waves with an electrostatic field amplitude that grows at the exponential rate ω i ≈ ω p,e /2γ(V b ). In the absence of other waves in the simulation, the electrostatic wave energy density 0 E 2 x (t) should grow at the rate 2ω i . To verify this, we compare the growth of the wave energy in all simulations. We make the wave growth comparable for the simulations (a)-(f) by normalizing 0 E 2 x (t) to the energy density E 0 = m e c 2 n e [γ(V b ) − 1] of one respective electron beam at the time t = 0. Figure 1 shows the time evolution of the electrostatic wave energy and confirms that all simulations show the correct exponential growth rate. The wave growth in the simulation (f) is delayed by lower initial noise levels, arising from the increased number of particles per cell [37].
The growth of the electrostatic wave energy is confined to a broad k-interval centred around k u . To show this, we first note that for the simulations (a)-(e) x = j x, t = l t s  and k = m k min are defined on a discrete grid and j, l and m are integers. The wave power is and where E x (x, t) is given in units of ω p,e c m e /e. Figure 2 exhibits that initially the electrostatic wave grows at the wavenumber k ≈ k u . At t ω p,e /γ(V b ) ≈ 20, the wave at low k saturates and collapses for the case γ(v b ) = 10 3 . The mean power E 2 x (k, t) is approximately constant up to large k after the collapse. It may thus be considered as a white noise along the space axis.
The normalized magnetic field energy density (B 2 y (t) + B 2 z (t))/µ 0 E 0 , where µ 0 is the permeability of free space, is displayed in figure 3. The simulation with γ(v b ) = 10 3 shows the strongest growth in the magnetic field energy for the simulation (e). By comparison, the simulation (f) shows much weaker magnetic fields, and it follows qualitatively the cases with γ(v b ) < 10 3 . The lower absolute value is due to the larger number of particles per cell, which results in lower statistical fluctuations of the currents. Since the curves representing the simulation (e) and the simulation (f) diverge at an early stage, the origin is likely to be associated with the particle initialization. The random initialization of the particle velocities may have introduced a deviation from a thermal equilibrium. Throughout the simulations, the electrostatic wave energy is four orders of magnitude larger than the electromagnetic energy. The impact of electromagnetic instabilities is thus small and probably negligible.

Particle distributions
The electrostatic and electromagnetic waves saturate at an energy level that is well below that of E 0 . The electrons and the protons also do not exchange substantial energy and the energy   densities they carry remain practically unchanged, as figure 4 shows. As in [25], the initial wave energy is provided by the electrons, while the proton beams increase their energy somewhat. The electrostatic wave saturates by trapping electrons, even for the high γ(v b ) that are considered here. Evidence is presented in figure 5 that shows the presence of periodic structures in the electron distribution. The period of these phase space structures is comparable to the wavelength x 0 = 2π/k u of the most unstable electrostatic wave.
Most energy is lost by the electrons after their initial trapping, that is during times when the broadband fluctuations in figure 2 have replaced the initial wave. A representative electron phase space distribution is shown in figure 6 for the simulation (e) with γ(v b ) = 10 3 . It shows electron jets similar to those discussed in [25]. This electron distribution further evolves under the influence of the broadband electrostatic noise that shows an almost white noise power spectrum in k. Since most electrons move at the speed v j ≈ c, they cross a distance x in x ω p,e /2πc ≈ 10 −2 , where ω −1 p,e is the characteristic fluctuation timescale. The electrons thus experience the electrostatic field, in their own frame, as a time series of white noise. In response to this, all simulations apparently develop electron momentum spectra that resemble relativistic Maxwellian distributions at low energies, as figure 7 shows. The electron momentum distributions agree with  high-energy tail expands to considerably higher momenta. Simulation (e) with γ(V b ) = 10 3 in figure 7 shows the most pronounced tail.
In figure 8, we compare the full momentum distribution, obtained from the simulation (f), with the relativistic Maxwellian distribution at two simulation times. The electron momentum distribution is a relativistic Maxwellian for both times in the interval −5 < p x /m e cγ(V b ) < 5. Outside this interval, we find an accelerating tail, in particular for p x > 0. This tail is formed by   a phase space structure resembling the ultrarelativistic equivalent of an electron hole [33,34] or an electron phase space jet [32], as figure 9 shows. This electron structure is the most energetic found in the simulation and it is still well represented by the simulation (f) at the time tω p,e /γ(V b ) = 46. We note the onset of small-scale momentum oscillations at the relative positions 10-40 in figure 9.A further evolution in time of these oscillations, which have a period comparable to the length x/2 of one grid cell in the simulation (f) is likely to introduce numerical artefacts. The simulation is thus stopped at this time. Such artefacts could not be found in the low-resolution simulations (a)-(e). The data from these simulations shows that the electron representation is not sufficient for a resolution of such low-density phase space holes.
The strong electron fluctuations in figure 5 and the ultrarelativistic electron phase space structures in figure 9 are associated with strong electrostatic fields. These fields, in turn, modulate also the protons as figure 10 shows for the simulation (f). We find momentum oscillations of the beam protons on all scales that are resolved by the simulation. The charge density fluctuations reach 20% of the mean density of the proton beam and show a broad spectrum N p (k) Note that the simulation (f) uses 2N g simulation cells. Here, ρ 0 is the average charge density of one proton beam and ρ 2 (x) is the charge density of the proton beam moving with p x > 0. This beam moves at a speed close to c and the Doppler shifted broadband fluctuation spectrum should constitute a key component of the white noise spectrum shown in figure 2. The charge density fluctuation spectrum, shown by (c) and normalized to the maximum value, is broad and covers the entire k-spectrum of the simulation.

Relativistic Brownian motion
In this section, we want to assess if the apparent relativistic Maxwellian distribution at low momenta shown in figure 7 is coincidental. A process that can produce such a distribution is the relativistic Brownian motion [29]. The relativistic Brownian motion requires that the system can be described by the relativistic Langevin equation where µ is the friction coefficient and w(t) is a stochastic increment that is equivalent to a Wiener process. For non-relativistic particle speeds a Wiener process implies that the stochastic increments w(t) in equation (14) are distributed according to the probability distribution with the parameter D. As discussed in [29], the equivalent of the distribution (15) for relativistic particle speeds is where γ(v) is the Lorentz factor of the considered particle. Equation (14) describes the timeevolution of a particle, for example one out of the 2N e electrons in our simulations (a)-(e). Plasma particles do not collide mechanically in the Vlasov-Boltzmann approach which the PIC code uses. The particles interact instead via the macroscopic electrostatic and electromagnetic fields. From figures 1 and 3, we find that for all simulations, except the simulation (e), the electromagnetic field energy is much less than the electrostatic one. This observation, together with the fact that the simulation (e) and the simulation (f) develop relativistic Maxwellian distributions despite the different electromagnetic energy densities, imply a negligible effect of the electromagnetic waves on the simulation results. We thus consider further only the electrostatic waves.
The simulation equivalent to the friction µ in equation (14) is the loss of electron energy to the charge density fluctuations associated with the electron motion. The increments w(t) in equation (14) are determined by the interaction of the particles with the electrostatic field E x (x, t), that is sampled by the PIC code in a fixed coordinate system, the laboratory frame in [29]. The increments are thus aligned with the x-axis and the one-dimensional approach in [29] should provide an adequate description of this system. Figure 2 further shows that the electrostatic fluctuation spectrum is that of white noise over a wide range of wavenumbers. This is a requirement for a relativistic Brownian motion of electrons.
The momentum increments w(t) for the electrons in equation (14) are associated with the fluctuations of E x (x, t) by the electrostatic term of the relativistic Lorentz force w(t)/dt = −eE x (x, t)/m e .We thus rewrite the Langevin equation to a form suitable for our simulations, i.e. dp x ( The Langevin equation can be transformed into an evolution equation for a phase space distribution, that is a Fokker-Planck equation Here, f(t, p x ) is the phase space distribution, in our case the distribution formed by all 2N e electrons in the simulations (a)-(e), and the term j(t, p x ) is the probability current. Using the Hänggi-Klimontovich approach from [29] to express the probability current, we obtain the relativistic Maxwell distribution given by equation (13) as the equilibrium distribution.
The changes in the average electron energy, computed by the PIC simulations, are quite low, as we see from figure 4. After t ω p,e /γ(V b ) ≈ 40, an equilibrium between the particle energy and the electrostatic field energy develops, as for thermal noise [37]. The electron momentum distributions, as shown in figure 7 at the time when the equilibrium between electrons and electrostatic fields is reached, resemble relativistic Maxwell distributions at low p x /m e c γ(V b ). If these are due to relativistic Brownian motion, we should expect E x (x, t) to be a Wiener process.
We obtain from E x (x, t) in physical units the unitless fieldẼ x (x, t) = E x (x, t)( 0 E 0 ) 1/2 . We calculate the normalized electric field derivative simulation (f) due to its different time step. By the division by γ 1/2 (V b ), we remove the explicit dependence of P r [w(t)] in equation (16) on γ(v b ). It is equivalent to the transformation w(t)/γ 1/2 = y(t) in [29].
The slowest wave growth occurs in the simulation with γ(v b ) = 10 3 , and here we measure the fields in the interval 48 < tω p,e /γ(V b ) < 54. That is well after the wave saturation at t ω p,e /γ(V b ) ≈ 20, and the magnitude and the mean of the electric field fluctuations are time stationary. For all other simulations the system has considerably more time to relax to an equilibrium. We count, for each γ(v b ) separately and for each position x, the number of values Ẽ x (x, t), that fall into specific intervals. We assume that P r [w(t)] is independent of x. We add the distributions of all N g simulation cells. We normalize the heights of the resulting curve for each γ(v b ) to its peak value. If the electric field increments represent a Wiener process with the probability distribution P r [w(t)], all curves should match. This is confirmed by figure 11, which shows that all curves are well approximated by a Gaussian envelope.

Discussion
In this paper, we have examined the growth and the saturation of electrostatic and electromagnetic waves driven by relativistically moving beams of electrons in an unmagnetized plasma. We have chosen two equally dense charge-and current neutral electron-proton plasmas that move at a relative Lorentz factor that has varied between 200 < γ(v b ) < 10 3 . Such beam speeds are thought to be representative for the expanding blast shell of GRBs. While the GRB shell is thought to be much denser than the ambient plasma, we may obtain a low-beam density ratio at the front end of the blast shell, where faster and thin shell fragments outrun the slower bulk plasma. Since the TS instability grows fastest if the colliding plasma species have a low density ratio, this scenario might be important.
We have modelled the plasma by means of a relativistic and electromagnetic PIC simulation. We have only resolved space along the beam flow direction by which we have excluded waves that develop obliquely to the plasma flow direction. This is a strong constraint that may not apply at GRBs, since oblique filamentation modes grow, in a homogeneous plasma, typically much faster than the TS modes that are aligned with the plasma flow [3,4]. The purpose of this work can thus not be to give an accurate and realistic account of the dynamics of plasma colliding at relativistic speeds, which requires multiple spatial dimensions. Instead, our ambition has been to examine the growth, saturation and collapse of waves with wavevectors that are aligned with a relativistic spatially homogeneous flow and this at a good simulation resolution.
In our simulations, we have not found evidence for the growth of electromagnetic waves via the Weibel instability. Probably the growth rates of this instability are too low for the chosen plasma parameters. We have found that the initial exponential rate of the TS modes is wellapproximated by the growth rate obtained from the linear dispersion relation. The electrostatic waves grow, however, over a wide band of wavenumbers. The TS waves eventually saturate by the trapping of electrons in the electrostatic wave potential. These non-linearly saturated modes are not stable and the TS waves collapse. We have analysed the power spectrum of the electrostatic turbulence that has developed out of the wave collapse, and have found it to be close to white noise. White electrostatic noise is a pre-requisite for the Brownian motion of particles in the wave fields. Our restriction to one spatial dimension and the dominance of electrostatic over electromagnetic fluctuations has then allowed us to compare a model of relativistic one-dimensional Brownian motion [29] with our simulation results. By doing this so, we have found that the distribution of electric field fluctuations represents a Wiener process, a second requirement for the Brownian motion which we have thus identified as the origin of the final relativistic Maxwellian electron distributions in all our simulations. The electrostatic fields are maintained by fluctuations of the electron and proton charge densities.
The restriction to one spatial simulation dimension has allowed us to use over 10 3 electrons per cell for 2N g cells in the simulation (f). We could thus follow this Maxwellian to high electron momenta. At high momenta, we have found a hot electron tail at the simulation's end. This electron tail is due to ultrarelativistic phase space holes or jets [32]- [34].
With respect to astrophysics we note that the dominant contribution to the electron acceleration has originated from the random walk of electrons in the stochastic wave fields, rather than by the interaction of electrons with monochromatic wave fields that are driven by the TS instability. If oblique modes are unstable and their collapse yields electrostatic turbulence with a spectrum similar to that found here, we would expect GRBs to consist of a bulk electron distribution that can be represented by a relativistic Maxwellian distribution plus, may be, a high-energy tail due to localized ultrarelativistic electron phase space structures.
Thermal particle distributions arising from the relaxation of ultrarelativistic plasma flows may have been observed at GRBs [30,31]. The simulation results obtained here may raise the question whether this thermal population could be due to a random walk of particles in electrostatic turbulence, similar to the relativistic Brownian motion. While the constraints imposed on our simulations are clearly too stringent to give a correct account of the physics at GRBs, the simulation results could have identified a process that can rapidly thermalize the relativistic plasma flow at GRBs.

Institute of Physics ⌽ DEUTSCHE PHYSIKALISCHE GESELLSCHAFT
In addition, the simulation results related to the high-energy tail may also serve as an orderof-magnitude estimate of the energies that electrons can reach by their interactions with the stochastic fields. Earth-bound observers typically measure the emissions of GRBs in a reference frame that is moving at mildly relativistic speeds relative to the ambient plasma in which the GRB source is immersed. Satellite experiments are thus likely to measure radiation in the beam frame rather than in the laboratory frame. The simulation (f) has found electrons with peak Lorentz factors ∼6 × 10 4 equivalent to electron energies of 30 GeV.
The electrostatic energy density in our simulations has reached the significant fraction 10 −2 of the electron kinetic energy density for the entire range of considered γ(v b ). These fields may give rise to electrostatic Bremsstrahlung processes by which some spectral components of the electromagnetic emissions by GRBs may be explained [38].
Beyond these findings, the evolution of the relativistic TS instability and the relaxing plasma has turned out to be a self-consistent potential test-bed for research into the motion of particles in stochastic fields and this might make its examination useful in areas outside plasma astrophysics, e.g. relativistic diffusion theory.
Future work has to examine stochastic wave-particle interactions for more realistic settings, i.e. incorporating multiple spatial dimensions as well as assessing the impact of spatial inhomogeneities of the plasma.