Electron surfing acceleration by mildly relativistic beams: wave magnetic field effects

Electron surfing acceleration (ESA) is based on the trapping of electrons by a wave and the transport of the trapped electrons across a perpendicular magnetic field. ESA can accelerate electrons to relativistic speeds and it may thus produce hot electrons in plasmas supporting fast ion beams, like close to astrophysical shocks. One-dimensional (1D) particle-in-cell (PIC) simulations have demonstrated that trapped electron structures (phase space holes) are stabilized by relativistic phase speeds of the waves, by which ESA can accelerate electrons to ultrarelativistic speeds. The 2(1/2)D electromagnetic and relativistic PIC simulations performed in the present paper model proton beam driven instabilities in the presence of a magnetic field perpendicular to the simulation plane. This configuration represents the partially electromagnetic mixed modes and the filamentation modes, in addition to the Buneman waves. The waves are found to become predominantly electromagnetic and nonplanar for beam speeds that would result in stable trapped electron structures. The relativistic boost of ESA reported previously is cancelled by this effect. For proton beam speeds of 0.6 and 0.8c, the electrons reach only million electron volt energies. The system with the slower beam is followed sufficiently long in time to reveal the development of a secondary filamentation instability. The instability forms a channel in the simulation domain that is void of any magnetic field. Proton beams may thereby cross perpendicular magnetic fields for distances beyond their gyroradius.


Introduction
Accretion discs form due to angular momentum conservation, as material is gravitationally attracted by a compact object [1]- [4], a neutron star or a black hole. The discs often subdivide into a cool component and a hot corona or wind [5]. Accreting compact objects are thought to be the origin of the jets of microquasars (MQs) [6,7]. The compact objects powering the MQs in our galactic neighbourhood contain about 10M s (solar masses) [8]. Their steady jets are ejected at the same time when a reduction of the total disc emission is observed [8]. The jets expand with subrelativistic speeds and the transient jets reach Lorentz factors of at least ≈ 2 [6,7,9].
While observational evidence of hot coronae and relativistic jets is available, the physical mechanisms behind their generation and the jet-disc coupling are still not clearly identified, except that they are due to radiation pressure and turbulence in the inner accretion disc [10]. In a steady state, the individual volume elements of the ionized inner disc probably follow regular orbits with velocities that are a function of the orbit radius. The presence of a magnetic field, as well as temperature gradients and radial velocity gradients, implies gyro-viscosity [11]. Viscous heating may explain the blackbody radiation of the steady state disc. Magnetohydrodynamic (MHD) instabilities during the unsteady phase of the disc result in a large-scale deformation of magnetic fields in the disc and in shock waves [12,13]. It has already been proposed [13] that the shock regions can be sources of the jet flow. The shock waves and the changing magnetic fields alter the orbits of neighbouring plasma volume elements, e.g. relative to the normal of a shock or a turbulent boundary layer [14]- [17], forcing some plasma volume elements to collide. In the absence of binary collisions, a thermalization of the colliding plasma volume elements can be achieved through particle scattering by waves. On an appropriate spatio-temporal scale, we can consider the electrons as magnetized, implying gyrotropic distributions, and the ions as unmagnetized beams. Large-amplitude waves may then be driven by ion beam instabilities.
In what follows, we consider only protons and assess how efficiently electrons can be accelerated by wave-particle interactions. A simple, yet powerful, electron acceleration mechanism is the electron surfing acceleration (ESA) [18]- [20], a type of plasma wave accelerator [21]. It has been introduced as a mechanism that can accelerate electrons in 3 laser-plasma interactions, if a spatially uniform constant magnetic field B 0 is present. Two electromagnetic laser pulses couple energy through a wave beat to an electrostatic wave with a wavevector k, which can trap electrons if its amplitude is sufficiently high. ESA assumes that B 0 ⊥ k. The trapped electrons are transported across B 0 with the phase speed of the wave. The cross-field transport continuously accelerates the electrons perpendicularly to k and to B 0 . This is a faster process than stochastic wave-particle interactions [22] or scattering by field fluctuations [23]- [25]. In principle, ESA could accelerate electrons to unlimited energies if the wave amplitude is sufficiently high so that the electrons cannot de-trap [20].
An electrostatic wave can also be driven by a proton beam that crosses a magnetic field at an angle of 90 • . Quasi-electrostatic upper hybrid waves grow, if the beam speed is nonrelativistic and they saturate by the trapping of electrons and beam protons [26]- [29]. The phase speed of the wave is close to the proton beam speed. The ESA mechanism can thus also occur in astrophysical settings, which we consider here. The characteristic wavelength of upper hybrid waves, and thus the scale size of the accelerator, is the electron skin depth if the proton beam is fast. Since the hot inner accretion discs of compact objects can be moderately warm and quite dense, e.g. 2.4 × 10 6 K and 10 12 cm −3 for the MQ GRS 1915-105 [30], ESA can operate on centimetre scales. ESA can thus serve as a small-scale energy dissipation and particle acceleration mechanism. The ESA and, more generally, kinetic instabilities may thus add a kinematic viscosity to the large-scale MHD evolution of the disc.
The original study of unlimited ESA in [20] and the test particle calculations in [31] employ an idealized model. Firstly, these studies assume a planar electrostatic wave. Secondly, the wave has a constant phase speed. Thirdly, the wave amplitude is constant and sufficiently high to trap electrons. The robustness of ESA against these assumptions has been tested. The impact of nonplanar wave fronts has been examined analytically in [32,33] and it has been found that a curved electrostatic potential can limit the peak energy accessible to the electrons. The limited electron acceleration by nonplanar quasi-electrostatic potentials has been confirmed by twodimensional (2D) nonrelativistic particle-in-cell (PIC) simulations [34,35]. The simulations also show that the electron motion across the time-dependent electric field gradients destroys the reversibility of the energy transfer from the wave to the electrons [26,34]. Another constraint arises, if the wave potential is not stationary in its rest frame, due to the wave damping and the nonlinear changes of the phase speed of the wave caused by the accelerating electrons [36,37]. Electron bunches, which have escaped the time-dependent wave potential, have been observed in [38,39], where the 1D geometry of the simulation boxes has excluded effects due to nonplanar wave fronts. This de-trapping occurs even if the coalescence instability [27] and the sideband instability [40] are suppressed [38]. The latter two instabilities destabilize the wave potential to the extent that ESA ceases to work even in an idealized 1D geometry, as long as the phase speed of the wave is nonrelativistic [38]. This is a possible reason for the absence of strong ESA at the Earth's bow shock, where proton beams move at nonrelativistic speeds across magnetic fields [41].
If the proton beams move at mildly relativistic speeds, which is realistic for the proton beams that are reflected by fast accretion disc shocks [12], the coalescence and the sideband instability become ineffective. In a 1D geometry, relativistic beam speeds v b > c/2 stabilize in PIC simulations the saturated upper hybrid wave and the trapped electron islands against secondary instabilities [42,43]. This enhanced stability, together with the relativistically strong v b × B 0 force, can result in a formidable ESA efficiency. The electrons in [43] have reached Lorentz factors of 10-100, and ESA can thus contribute relativistically fast electrons to accretion 4 disc jets and to the hot disc halo population. Hence, ESA merits a further examination in more realistic settings.
The purpose of this paper is to expand the 1D PIC simulation study in [43] to a 2D geometry, which can take into account the filamentation and mixed mode instabilities [44]- [48] that compete with the quasi-electrostatic upper hybrid mode, resulting in a higher-dimensional electric field modulation as well as in growing magnetic fields. To the best of our knowledge, the full 3D k-spectrum of unstable waves, which are driven by proton beams moving orthogonally across B 0 and through a background plasma consisting of electrons and protons, has not yet been explored. In section 2, we thus examine with a fluid model the interplay of the upper hybrid instability with the magnetized oblique and filamentation mode instabilities in the 3D k-space. We examine in section 3 with 2(1/2)D PIC simulations (two spatial and three momentum components) the nonlinear evolution of the competing instabilities. We consider beam speeds of v b = 0.6 and 0.8c. The beam speed is 2-3 times higher than that in the 2(1/2)D simulation in [34], allowing us to examine the impact of the relativistically boosted efficiency of ESA. These speeds, the beam density and the magnetic field strength are close to the values used in [43] and the selected speeds are the largest that do not trigger a strong electromagnetic finite grid instability (FGI) [49] with the feasible simulation resolution. We examine only these two case studies due to the high computational cost of the PIC simulations.
We find a new limiting mechanism for the ESA that constrains to a few million electron volt the maximum energy the electrons reach. This previously unreported and important limiting mechanism arises from the electromagnetic field components of the growing waves and of their harmonics. These fields quickly grow to an amplitude that equals |B 0 |, by which the ESA is cancelled and which permits the proton beam to breach the magnetic barrier imposed by the initial perpendicular magnetic field. In section 4, we discuss the results and their astrophysical relevance.

Initial simulation conditions
We consider a spatially homogeneous plasma of infinite extent, consisting of two proton beams moving along the x-direction through a background plasma (see figure 1). The latter consists of mobile electrons and protons. All species have single-Maxwellian velocity distributions in their respective rest frames. The electron plasma frequency is e = (n e e 2 /m e 0 ) 1/2 , where n e and 0 are the electron number density and the dielectric constant, respectively. The plasma frequency of the background protons and of each of the two proton beams is p = e m e /3m p in the rest frame of the electrons. The total proton density equals the electron density. The electron thermal speed v e = (K B T e /m e ) 1/2 = 10 7 m s −1 and the electron temperature is thus T e = 6.6 × 10 6 K (570 eV). The background protons and the proton beam moving in the positive x-direction have the temperature T p = T e . This temperature is comparable to that of the blackbody radiation some MQ accretion discs emit [6,7]. The second proton beam that moves in the negative x-direction cancels the current of the cool beam. It has the temperature T = 10 4 T p , which often introduces an asymmetry in the wave growth. Such an asymmetry is realistic, because typically ESA is driven by a single wave. However, it is not a priori known which wave saturates first. While cooler beams result in larger growth rates of the instability, the hotter beam provides higher noise levels [50] and thus a larger initial amplitude for the unstable wave, The considered geometry: the proton beams initially move along the x-axis in opposite directions. The bulk plasma is at rest. The magnetic field points along the z-direction. The red arrow denotes a general wavevector in 3D space. For a certain interval of the wavevector space, the associated wave is unstable and grows. The linear dispersion relation is solved in the full 3D wavevector space to identify this interval. In the simulations, the spectrum of unstable waves is restricted to the x-y-plane.
shortening the growth time to its nonlinear saturation. The magnetic field B 0 = (0, 0, B 0 ) with eB 0 /m e = e /10 equals that in [34,43]. The electrons and background protons are initially at rest. The two counter-propagating proton beams move with the mean speed modulus v b . We consider the cases v b = 0.6 and 0.8c.

Linear dispersion relation
Due to the orientation of the various components of the initial plasma configuration (see figure 1) and the lack of symmetry, the 3D problem cannot be simplified to any 2D one. Due to the analytic complexity involved in considering the full 3D wavevector space, we perform a linear analysis in the cold fluid approximation. This approach should describe properly the linear dispersive properties of our plasma, because the thermal velocity spreads are nonrelativistic and small, albeit not very small, compared to v b . Starting from the conservation and force equations of a linearized, relativistic cold fluid and the Maxwell equations, the full 3 × 3 dielectric tensor is symbolically calculated through a straightforward adaptation of the Mathematica notebook described in [51] and appended to this paper in terms of the dimensionless variables  Figure 2 displays the growth rates for v b = 0.6 and 0.8c and for Z x > 0. Apparently, a weak B 0 does not modify the linear growth rates in the y-z-plane. Unstable modes occur with wavevectors that have Z x ∼1 and Z 2 y + Z 2 z 0. We can estimate the growth rate of these oblique modes by extracting the x x component of the tensor and by neglecting the weak magnetic field with B = 1/10. This approximation yields the largest growth rate in units of e along the Z x -axis as The maximum value δ M of the growth rate in the oblique direction is then derived by multiplying this result by 2/3 [52], With R = 1/1836, we find δ M (β = 0.6) ∼ 0.0361 and δ M (β = 0.8) ∼ 0.0328, both in excellent agreement with the numerical values δ M (β = 0.6) = 0.0357 and δ M (β = 0.8) = 0.0327. Figure 3 displays for various v b the value of the growth rate for Z y = 0 and Z x = 1 in terms of Z z . The growth rate reaches its asymptotic value given by equation (3) as soon as Z z > 2. We recover here an effect that has already been identified for electron beams [52]. Due to their increased relativistic inertia along the beam direction, the particles oscillate more easily in oblique directions and, thus, the most unstable Z are not parallel to v b . Even for the fairly dense beams with n b = n e /3, these modes are found for ω ∼ e , so that the component of their phase velocity along the beam axis is almost equal to the beam velocity. These are quasielectrostatic Buneman-type modes [53] that arise from the interaction of the proton beams with the electrons of the bulk plasma. The recent linear kinetic (or warm fluid) model of oblique modes [44,54] predicts that thermal effects turn the continuum into a single oblique wavevector by damping the growth of waves with large Z . These predictions have been verified by PIC simulations [49,55]. We know that these modes are quasi-electrostatic and result in a filamentation of the beam [44]. Concerning the present investigation, we restrict ourselves to the case of a cold plasma and leave the identification of the most unstable mode to the PIC simulations. Regardless of the obliqueness of the fastest growing mode, it is quasi-electrostatic with Z x ∼ 1, so that one component of the ESA mechanism we consider is assured. The second component, the magnetic field B 0 , is also present.
The largest growth rate along the Z y -and Z z -axes can be investigated noting that they saturate at large Z . By developing the dispersion equation along these axes for large wavevectors, we find the maximum growth rate of one absolute unstable mode along the magnetic field axis z, It is worth noting that this quantity is exact and independent of B . The plasma filamentation along the magnetic field direction is unaffected by B . Turning now to the y-axis, we find a polynomial equation which cannot be solved in the general magnetized case. Nevertheless, it is easily checked that for B = 0 the growth rate (4) is retrieved, which was expected since the nonmagnetized system must be invariant by rotation around the beam axis x. Furthermore, it can also be determined that the growth rate along y vanishes exactly for With R = 1/1836, this gives limit values 29 and 60 for B and for β = 0.6 and 0.8 respectively, well beyond the value considered here.
Each CP corresponds to a phase space element and, thus, to many particles. While the mass and charge of a CP are larger than those of the physical particles it represents, their ratio is identical. Our PIC scheme [57] fulfills ∇ · E = ρ/ 0 and ∇ · B = 0 to round-off precision. The 3D electromagnetic fields are evolved on a 2D grid with periodic boundary conditions in all directions. A cell size x = y in the x-and y-directions is selected. The width of the single cell in the z-direction is also x .

Simulation set-up
We employ in both simulations a rectangular box with the side lengths L x = L y and L x = 6πv b / e in the x-and y-directions. This resolves three full wavelengths of the Buneman-type wave and thus three of the electron phase space holes that constitute the nonlinearly saturated wave in the x-direction. Even more oscillation periods of the wave in the y-direction are feasible, since the most unstable waves have a Z y 2. The presence of multiple electron phase space structures allows for the development of sideband and coalescence instabilities, the key limiting mechanisms for ESA [38]. The previous section has shown that the spectrum of unstable waves is 3D and our 2(1/2)D simulation is thus not completely realistic. However, since the magnetic field is parallel to the z-direction, the particle acceleration takes place in the x-y-plane. The electron speed along the z-direction remains low and the field gradients in this direction are less important than those in the x-y-plane.
To suppress the FGI [49], we set a cell size x = λ e /57.5 in simulation (a) with its v b = 0.6c and x = λ e /80 in simulation (b) with v b = 0.8c. The electron skin depth λ e = c/ e is identical in both simulations. We require a finer resolution in simulation (b), because the strength of the electromagnetic FGI increases with v b . We resolve the simulation planes with 650 2 cells and with 1200 2 cells in simulations (a) and (b), respectively. The electrons are represented by 200 particles per cell (PPC) in simulation (a) and by 128 PPC in simulation (b). Each of the proton species is modelled by 128 PPC in (a) and 72 PPC in (b). We employ the mass ratio m p /m e = 1836. By restricting the simulation duration to T s m p /eB 0 , we can neglect the proton magnetization and assume a spatially homogeneous directed proton beam. On longer timescales, a ring distribution would be more appropriate [58].
Insight into the plasma dynamics is given by the box-averaged energy densities of certain quantities. The electron acceleration results in a growing E K = L −1 , where j is the Lorentz factor of the jth computational electron with the mass m CP . We introduce E 0 as the initial (thermal) electron kinetic energy density. The total number N p of CPs is N p = 8.45 × 10 7 in simulation (a) and N p ≈ 1.85 × 10 8 in simulation (b). The electrostatic component of the waves 650 and 1200 for simulations (a) and (b), respectively. The total electric field energy density incorporates also the electromagnetic fields of the mixed/filamentation modes and the electric fields due to the currents of the accelerated electrons. The contribution of E y to E E exceeds that of E z in the simulations by 3-4 orders of magnitude. The contribution of E z to E E is negligible, at least in the considered 2(1/2)D geometry. The mixed/filamentation modes also have a magnetic field component, which results The energy densities in the simulation with its total runtime e t ≈ 900 are shown in figure 4.
The E Ex and E E develop similarly until e t ≈ 600. The growth of E E prior to e t ≈ 200 is mainly due to the obliquely propagating electromagnetic waves driven by the FGI [49]. A second growth phase is evident after e t ≈ 200, which can be attributed to the physical streaming instability. The E E grows at 0.8 times the rate 2δ M / e = 0.072, where δ M corresponds to the fastest-growing wave estimated in the previous section. Note that the energy density grows twice as fast as the δ M of the wave amplitude. The growth rate reduction is firstly due to the averaging over a wide range of waves and secondly due to the nonzero temperature of our plasma.
The larger E E indicates that the waves have an electromagnetic E y component, which is characteristic of mixed modes. This also explains the simultaneous growth of E B . The growth of electromagnetic fields will modify the ESA, that is based in its original form on electrostatic waves. However, the exponential wave growth during 200 < e t < 300 and the abrupt saturation of E Ex and E E has been evidenced also by a previous 2D simulation [34] and is indicative of the two-stream instability and electron trapping by electrostatic potentials. If the plasma dynamics is quasi-electrostatic, ESA may still work. The electron trapping yields the fast growth of E K at e t ≈ 250.
At the time e t = 300 the contribution to E E of the high-k oscillations due to the FGI is 25% of that of the physical oscillations. The oscillations driven by the FGI are thus not negligible. Their broad k-spectrum [49] and their short wavelengths imply, however, that they will probably not interact resonantly with electrons, in contrast to the long physical waves. The phase speed of these high-frequency electromagnetic waves is just below the speed of light due to the selected numerical approximation of the Maxwell's equations [49] and they will scatter some of the relativistic electrons. After e t ≈ 300, the E K continues to grow. Initially, this could be due to ESA, since the E K grows similarly to the electron energy in [43]. The scattering by the high-k oscillations may also contribute. After e t ≈ 600 the curves for E Ex and E E diverge, which suggests the development of a secondary instability.
We examine the field distributions calculated by the PIC simulations at the times of the initial wave saturation and after the secondary instability has saturated, in order to compare these self-consistently computed wave fields with the electrostatic plane waves, on which ESA is based. Figure 5 displays the field distributions at e t = 300. The x-and y-positions are normalized to L u = 2 π v b / e , the wavelength of the most unstable electrostatic two-stream mode. The low-pass filter removes waves with wavelengths less than 3L u /50. The strengths of the E x and E y components are comparable. The electromagnetic component of the mixed modes is thus not negligible, as the comparison of E Ex with E E has already shown. All waves propagate on the branch with k x /k 0 = 1 and k 0 = 2π/L u . They are thus not filamentation modes, which would have a wavevector strictly perpendicular to the beam flow direction and thus k x ≈ 0. The amplitude of the electric structures is below cB 0 and unlimited ESA [20] is thus not possible. The power spectrum of the electric fields reveals a wide band of waves along the k y -direction. The wide k y band of the waves with k x ≈ k 0 has been predicted in the previous section. The presence of these waves and their spatial harmonics explains the patchy structures of all fields along the y-direction. Note, however, that for any given y we obtain three oscillation periods in E x , E y and B z . Figure 5 also reveals oscillations of the B z component with an amplitude comparable to B 0 . In addition to the electromagnetic wave component, the spatial harmonics constitute a modification of the assumptions, from which ESA has been derived in [20]. The modulations of the fields along the y-direction imply that the ESA mechanism based on plane waves, as proposed by Katsouleas and Dawson [20], will not work in this form in our simulation. Since the fields are homogeneous over small spatial intervals, the ESA may nevertheless accelerate electrons for limited durations. This constraint is similar to curvature effects [32,33].
The fields in figure 6 at e t = 800 reveal that the late growth of E E in figure 4 has changed qualitatively the field structure. The electrostatic E x component still shows waves show the E x , E y and B z components, respectively. All fields are normalized to cB 0 and space is normalized to L u . Panel (d) shows the 10logarithm of the spatial power spectrum of E x + iE y , normalized to its peak value.
with the length L u . Both the E y and B z components now display oscillations along the y-direction. The power spectrum of the electric field peaks at k x = 0. The growth of E E is thus due to the filamentation instability. Note that the solution of the linear dispersion relation in the previous section does not predict the growth of the filamentation instability. It is thus a secondary wave, which is destabilized by the plasma distribution after the initial wave saturation. The filamentation instability spatially modulates along y the charge densities of the bulk protons, the cool proton beam and the electrons. The number density distributions of the latter two species, integrated over x, and the B z , which has been averaged over x, are displayed in figure 7. The magnetic pressure gradient triggered by the spatially inhomogeneous B z gives rise to the E y field, which acts as to restore quasi-neutrality [59]. The filamentation instability has practically removed the B z component in some y-intervals. The protons can flow ballistically in these fieldfree channels, in which the ESA mechanism will not work due to the absent perpendicular magnetic field.

Beam speed
The increase of v b from 0.6 to 0.8c results in a stronger relativistic reduction of the growth rate of the purely electrostatic mode with k y = 0 compared to the partially electromagnetic mixed modes. We thus expect that E E is relatively stronger compared to E Ex . Figure 8 confirms this. Until e t ≈ 120, the curves for E Ex and E E grow in unison, with E E being dominant. We Figure 6. The low-pass filtered fields for v b = 0.6c at e t = 800: the panels (a)-(c) show the E x , E y and B z components, respectively. All fields are normalized to cB 0 and space is normalized to L u . Panel (d) shows the 10logarithm of the spatial power spectrum of E x + iE y , normalized to its peak value.
attribute this to the FGI as well as to the rapidly growing electromagnetic modes. In contrast to figure 4, no clear exponential growth of E Ex or E E is detectable after this time except, possibly, in the interval 230 < e t < 330. This is probably due to a higher absolute strength of the FGI. The change in the growth rate at e t ≈ 200 in figure 4 took place at an energy density well below E 0 . In figure 8, we find a change in the growth rate at e t ≈ 150, at which E E is comparable to E 0 . The relative contribution of 30% of the FGI to E E at e t = 330 is barely stronger than in the simulation with v b = 0.6c. The physical instabilities outgrow the FGI at energies, for which nonlinear processes have already started. The field energy densities are thus not a good measure for the wave growth during the linear stage in the simulation with v b = 0.8c, at least in comparison to that with v b = 0.6c.
The curve for E Ex in figure 8 shows a cusp at e t ≈ 330. At around this time, the electron kinetic energy E K experiences its strongest growth. We do not, however, find the same simultaneous drop of E Ex and E E after the saturation, as we did in figure 4. The relativistic stabilization of the islands of trapped electrons [43] may be responsible. This cannot, however, be the only cause because E E also converges to E B , which has been attributed to the development of a filamentation instability in the simulation with v b = 0.6c. We thus suspect that the nonlinear saturation mechanism changes from being quasi-electrostatic for v b = 0.6c to being quasi-electromagnetic for v b = 0.8c. Figure 9 provides further insight by showing the low-pass filtered electric and magnetic fields. The positions are normalized to L u , which has changed compared to the previous  show the E x , E y and B z components, respectively. All fields are normalized to cB 0 and space is normalized to L u . Panel (d) shows the 10-logarithm of the spatial power spectrum of E x + iE y , normalized to its peak value. subsection due to the higher v b . Waves with wavelengths below 3L u /50 are filtered out. The E x field shows the presence of at least partially electrostatic structures. This is, provided that the flow-aligned E x field is related at this time to charge density modulations, as it is the case during the linear growth phase of the waves. Figure 9 further reveals the presence of strong oscillations also in the E y and B z fields. The power spectrum of the complex electric field evidences waves with k x /k 0 = 1, as well as with k x = 0, and their harmonics. The electric field amplitudes exceed cB 0 . Unlimited ESA may, however, not work because spatial oscillations of all field components and strong electromagnetic wave components have not been considered by [20].
The comparison of figure 5 with figure 9 reveals no striking qualitative differences between the fields in both simulations. It is thus not obvious why the nonlinear saturation mechanism differs for both v b as suggested by figures 4 and 8. The structures in figure 9 have a smaller size in terms of L u and a larger amplitude in units of cB 0 compared to figure 5. On one hand, the accelerating electrons cross the electrostatic structures more rapidly along the y-direction, which decreases the probability for them to get trapped. Trapping in the wave potential is the electrostatic saturation mechanism, which is also one cause for ESA. On the other hand, the larger E y fields may deflect them more easily into current channels, which is the electromagnetic saturation mechanism.
The time evolution of the energy densities in figure 8 after e t = 330 reveals no drastic change. We thus do not anticipate major qualitative changes in the field structure after this time. show the E x , E y and B z components, respectively. All fields are normalized to cB 0 and space is normalized to L u . Panel (d) shows the 10logarithm of the spatial power spectrum of E x + iE y , normalized to its peak value.
The low-pass filtered fields in figure 10 support this claim. We find similar patchy structures in figures 9 and 10, which contrasts the drastic changes observed in figures 5 and 6. The electric fields in figure 10 are generally weaker than those in figure 9 and the size along the y-direction of the magnetic fluctuations is larger, as well as their amplitude. We find that the convergence of E E and E B is thus not an unique indicator for a development of the filamentation instability, because the final states in figure 6 and 10 are different. It is actually possible that the filamentation instability will also grow in the simulation with v b = 0.8c. To follow the instability further is, however, prohibitively expensive in terms of computing time.

Electron acceleration
The previous section has revealed significant differences between the electromagnetic fields calculated self-consistently by PIC simulations and the planar electrostatic wave fields, which have resulted in efficient electron acceleration in [20,43]. We thus anticipate differences also in the electron acceleration, which we quantify here. Figures 4 and 8 demonstrate that the electron energies E K grow to quite similar values during the simulations. The quantity E K is ambiguous, since it is not known how many electrons are accelerated. We thus sample the electron energy distributions at the times when the fields were measured. The kinetic energy of the computational electron with index j is E j = (1 − v 2 j /c 2 ) −1/2 − 1 in units of m cp c 2 . We evaluate the kinetic energy using the full velocity vector (v x , v y and v z ). The energy of each CP  Figure 11 shows the electron energy distributions. The typical energies that the electrons reach in both simulations are comparable and limited to about 1.5 MeV in physical units. The electrons are accelerated to mildly relativistic energies when the waves saturate and gain only moderately in energy thereafter. These energies are comparable to what electron trapping by electrostatic waves alone can achieve [42]. We could also not find in the particle data the phase space cones that are characteristic for ESA [20], and which are formed by the trapped electron oscillations in the x, v x -plane and the v y -component that increases due to the electron transport across B 0 . The ESA mechanism is thus not efficient for the considered plasma configuration. The final electron distribution measured in the simulation with v b = 0.6c is decreasing faster than exponential. This is the characteristic final distribution of the filamentation instability [59]. The final electron distribution in the simulation with v b = 0.8c is better approximated by an exponential. It is likely that it will also evolve into a more rapidly decreasing function, if the filamentation instability is triggered.

Discussion
In this work, we have continued previous studies of the ESA mechanism, originally proposed by [18,19] and expanded to relativistic regimes by [20], to incorporate mildly relativistic proton beam speeds and a second spatial dimension. The ESA mechanism is a type of plasma wave accelerator [21] that has been considered, for example, in the context of electron acceleration by laser-driven plasma waves [20] and also in astrophysical particle acceleration [38]. The ESA mechanism is based on a strong electrostatic wave that has trapped electrons. The electrons are transported by the wave with its phase speed across a perpendicular magnetic field. The cross-field transport accelerates electrons.
The motivation for this work has been to test whether the high phase speeds of the waves, driven by mildly relativistic proton beams, also result in stable electron phase space holes [42] in 2D. Such a stabilization boosts the ESA efficiency [43] as 1D PIC simulations have evidenced. The second dimension takes into account curvature effects of the wave fronts [32,33] and the richer wave spectrum [44], which has been considered for nonrelativistic beam speeds in previous 2D PIC simulations [34,35]. We have selected the proton beam speed v b = 0.6c. It results in waves with a phase speed that is high enough to show a weak relativistic stabilization of the electron phase space holes in 1D [42,43]. The second simulation has considered a beam speed v b = 0.8c, because for this v b the electron phase space holes have been quite stable in 1D [42].
Our simulations have demonstrated that the spatial structuring of the wave fronts prevents a strong electron acceleration by ESA. This limitation is similar to that discussed in [32,33] that took into account a curvature of the electrostatic potential in a second dimension. The field changes we observe here occur on scales less than an electron skin depth and involve both the electric and the magnetic fields. The PIC simulation demonstrated that the wave evolution is similar to the (quasi-)electrostatic case, if v b = 0.6c [34]. The proton beam-driven waves undergo an exponential growth phase, they saturate by the trapping of electrons and then they collapse. The wave collapse de-traps the electrons and ESA ceases to work. The saturated waves should be stabilized in the second PIC simulation by the beam speed v b = 0.8c, but here the wave spectrum has been dominated by electromagnetic waves. This electromagnetic component apparently cancels the ESA, which is based on electrostatic waves. We conclude that ESA has been inefficient in the considered case study also for relativistic beam speeds, which severely limits the validity of the conclusions of Dieckmann et al [43].
Our case study has only examined waves that have been driven by counterpropagating cool proton beams, which saturate by the trapping of the bulk electrons. Our conclusions are thus not generally valid. The importance of the electromagnetic wave components relative to the electrostatic wave components could be reduced by introducing a high-proton beam temperature perpendicular to the beam's flow velocity vector. A stronger background magnetic field can also suppress the filamentation and mixed mode instabilities, which have outgrown in our simulations the purely electrostatic Buneman wave. Our linear dispersion analysis has shown that the plasma and wave structuring can, in principle, be suppressed in the plane perpendicular to the magnetic field vector, which is the most relevant one for ESA. The required field strengths are, however, quite high and at least unlimited ESA that requires electric fields stronger than cB 0 may not be possible. Future studies have to examine how hot electron sub-populations influence the wave growth and stability and thus the ESA efficiency. The waves may also not originate in proton beam instabilities. They could, for example, be produced by a beat of other waves [20]. However, our simulations clearly demonstrate that ESA cannot be considered to be a plasma wave accelerator that always accelerates electrons to highly relativistic speeds, but which may do so in some cases that have yet to be explored.
While the final electron energies have been low compared to those in [43], the electron Lorentz factors of a few imply that ion-beam driven kinetic instabilities and, possibly, ESA could still be a significant energy dissipation mechanism in relativistically colliding plasmas. Shocks and turbulent boundary layers that move at mildly relativistic speeds through accretion discs [12] should give rise to a hot electron population with million electron volt temperatures. These electrons could form a component of the jet plasma and of the hot disc halo [13]. The electrons could further contribute through synchrotron emissions to the electromagnetic disc radiation [6]. The scale at which these beam instabilities work is about an electron skin depth, which is less than one metre at the inner accretion disc of the MQ GRS 1915-105. Such accelerators thus provide energy dissipation and particle acceleration on a scale that is practically infinitesimally small compared to the accretion disc size. Wave-particle interaction might thus complement magnetic turbulence as the viscosity term in the large-scale disc dynamics.
Our simulations have revealed a constraint on ESA other than the well-known instability of electron phase space holes [27] and the presence of modes other than the Buneman waves during the linear growth phase of the instability. We could follow in the simulation with the beam speed v b = 0.6c the time-evolution of the system to the point when a secondary filamentation instability sets in. At this time, the electrons are already relativistically hot and the system differs from the initial configuration, which has been dominated by the mixed mode waves. The plasma and wave structuring perpendicular to the proton beam velocity vector due to this secondary filamentation instablity will stop ESA. The secondary instability is strong. We have observed flow channels, in which the strength of the perpendicular magnetic field has been reduced to zero. At least mildly relativistic proton beams can thus drill themselves into magnetic field barriers. A consequence of this experimental finding is, for example, that the model of planar, nonrelativistic and perpendicular plasma shocks [15,16] may not work for the high plasma collision speeds observed in many astrophysical regimes. This model assumes that a perpendicular shock undergoes cycles, due to the spatial confinement of the charged particles. The structure is convected by the shock reflected ion beam. Each cycle results in a spatial propagation of the shock front by about an ion gyroradius. If ions can overcome with the help of the filamentation instability the barrier imposed by the perpendicular magnetic field, the shock dynamics would be modified. The secondary filamentation instability has further accelerated electrons. The acceleration efficiency has been weak and the final electron energy distribution decreased faster than exponential, in line with previous findings [59].
For the plasma parameters we have used here, the energy of the waves driven by the FGI has not been negligible. It has amounted to about 20-30% of the energy contained in the waves originating from the beam instability, when the latter have saturated. The wavelengths unstable to the FGI are comparable to the grid cell size, which is the electron Debye length. This short wavelength, the broad spectrum of unstable waves [49] and their large amplitudes imply that the FGI probably introduces artificial electromagnetic diffusion into the simulation, just like electrostatic diffusion is introduced by charge density fluctuations [23,24]. The short scale waves are pumped by the proton beams. Their high frequency and short scale implies that they will not act back on the proton beam, in contrast to the long waves that are driven by the physical beam instability. Only the latter can thus couple the electrons to the protons, which is required for a strong acceleration of the electrons. However, the short-scale waves can scatter electrons in and out of wave potentials, which will influence the evolution of the system. These waves may also yield a slow growth of the electron energies, like that we have observed at late times in the simulations. More likely is, however, that this energy increase is due to short ESA events or stochastic interactions between the electrons and either the electromagnetic fields [22] of the long physical waves or the magnetic structures found in particular for v b = 0.8c. Future work will, however, have to examine in more detail and quantify the impact of the FGI.