Simulation study of surfing acceleration in magnetized space plasmas

We present a numerical study of the surfing mechanism in which electrons are trapped in Bernstein–Greene–Kruskal (BGK) modes, and are accelerated across the magnetic field direction by the Lorentz force in magnetized space plasmas. The BGK modes are the product of an ion-beam Buneman instability that excites large-amplitude electrostatic upper-hybrid waves in the plasma. Our study, which is performed with particle-in-cell (PIC) and Vlasov codes, reveals the stability of the BGK mode as a function of the magnetic field strength and the ion beam speed. It is found that the surfing acceleration is more effective for a weaker magnetic field owing to the longer lifetime of the BGK modes. The importance of our investigation to electron acceleration in astrophysical environments has been emphasized.


3
Institute of Physics ⌽ DEUTSCHE PHYSIKALISCHE GESELLSCHAFT fashion. We describe the numerical setup and physical parameters in section 2. In section 3 we discuss numerical results from the Vlasov and PIC simulations. Conclusions and applications of our results are discussed in section 4.

Physical model and numerical setup
Collisionless shocks exist both at the bow shock of the Earth, where the surrounding solar wind streams with supersonic speed and collides with the geomagnetic field, and in the vicinity of supernovae where the plasma is ejected and propagates with high speed into the interstellar medium. Ions can be reflected either by the electrostatic potential of the shock or by the compressed magnetic field downstream the shock, to form a beam of ions in front of the shock. The beam may then be rotated in the upstream magnetic field and returned into the shock as a counterpropagating beam [24]. The ion beams can interact with electron waves and give rise to a Buneman instability where waves that have the same phase speed as the beam grow exponentially. The nonlinear saturation of the electrostatic wave occurs when electrons get trapped at the potential maxima of the wave potential, and subsequently BGK modes are formed.
In our numerical model, we use as an initial condition, two counterpropagating ion beams, which are propagating relative to the bulk plasma with the speeds ±v b . The electrons are assumed initially to form a homogeneous Maxwell-distributed background with a temperature T e and a zero mean velocity. The ion beam propagating in the positive x direction with speed v b (relative to the mean speed of the electrons) is assumed to have a shifted Maxwellian particle distribution with the temperature 10 T e , while the counterpropagating ion beam (with speed −v b ) has the temperature 100 T e . The bulk ions have a temperature of 10 T e . We assume that both the bulk and beam ions consist of protons with mass m i = 1836m e , where m e is the electron mass. The particle number densities of the two ion beams are each taken to be n b = n 0 /6, where n 0 denotes the unperturbed electron number density, while the bulk protons have the density n bulk = 2n 0 /3 so that the unperturbed plasma is charge neutral. In our simulations, we have set spatial box length L = 2πv b /ω pe , which corresponds approximately to the wavelength of the most unstable Buneman mode. Here, ω pe = (n 0 e 2 / 0 m e ) 1/2 is the electron plasma frequency, e the magnitude of the electron charge and 0 the electric permittivity in vacuum. The electrostatic and nonrelativistic Vlasov simulations solve the Vlasov-Poisson system where B 0 is the constant external magnetic field, and where the charge density is ρ = e(n i − n e ) and the particle number densities are obtained from the continuous particle distribution functions f j (x, v x , v y , t). Here, j equals e for electrons and i for ions, q e = −e and q i = e, where e is the magnitude of the electron charge, andx,ŷ andẑ denote the unit vectors in the x, y and z directions, respectively. The values of f j are represented on a fixed grid, and the numerical method is based on a Fourier method in velocity space, a pseudo-spectral method in space and a fourth-order Runge-Kutta scheme in time [21,22]. The PIC simulation code represents the particle distribution function with a large number of discrete superparticles, where each superparticle represents a large number of real electrons or ions. The electromagnetic field is represented on a fixed equidistant grid in space, and the value of the electric and magnetic fields are interpolated to the position of each particle. Then the momentum and position of particle k is advanced via the equations of motion dp (k) where where Q j and M j are the charge and mass, respectively, of the superparticles. From the positions and velocities of each particle, the charge and current densities ρ and j are calculated via an interpolating scheme and entered into Maxwell's system of equations and which is advanced in time with the initial conditions B = B 0ẑ and ∇ · E = ρ/ 0 for the magnetic and electric fields. In our one-dimensional case, we have ∇ =x∂/∂x. The PIC code is based on a Lagrangian scheme [20] which has a relatively low numerical noise level. In the physical problem, the box length and particle velocities are relatively small so that relativistic and electromagnetic effects will play only a minor role, thus making the two codes comparable. The Vlasov code uses a maximum representation of velocity space of v max = 20πv Te ≈ 63v Te , which is resolved with 400 intervals, where v Te = (k B T e /m e ) 1/2 is the electron thermal speed and k B is the Boltzmann constant. The space is resolved with 50 intervals, and the time step is adjusted dynamically to maintain numerical stability. The PIC simulations use 180 intervals to resolve the space, 2.7 million particles for the electrons and 432 000 particles for each proton species (the bulk protons and the two beams).

Simulation results
We first study the linear growth rates of the exponentially growing ∝ exp(γ t) amplitude of the Buneman instability. The results are presented in figure 1, where one can see that the Buneman instability is strongly dependent on the beam velocity. The growth rate was not large enough to be measurable in the numerical experiments for beam speeds below v b ≈ 7v Te , while for larger speeds the growth rate approached the theoretical value for a cold ion beam in a cold electron plasma, as given by [25] where ω pb = (n b e 2 / 0 m b ) 1/2 denotes the plasma frequency of the beam. Here, we have ω 2 pb /ω 2 pe = (n b /n 0 )(m e /m b ), where n b /n 0 = 1/6 and m e /m b = 1/1836, yielding the value γ cold /ω pe ≈ 0.03, indicated by the dashed line in figure 1. We observe that the linear growth rates obtained from the Vlasov and PIC simulations are in excellent agreement.
We next turn to the nonlinear saturation of the Buneman instability. When the electrostatic wave reaches a large enough amplitude, electrons get trapped in the potential maxima of the wave, and at this time the transfer of energy from the ion beam to the wave becomes less effective, and the growth of the wave amplitude gets reduced or halted. In figure 2, we have compared the amplitude of the electric field as a function of time for different beam speeds and magnetic fields, latter represented by different values of the electron gyrofrequency eB 0 /m e . We clearly see the initial exponential growth phase of the instability, followed by the saturation where the electric field reaches a maximum value, and then a more or less rapid decrease of the electric field. In general, the growth of the wave starts to saturate at a somewhat earlier stage (at a smaller amplitude) in the Vlasov simulation than in the PIC simulation. This is most clearly seen in the lower panels for the weak magnetic field case where the electric field in the Vlasov  Figure 2. The peak electric field as a function of time, obtained from the Vlasov (left panels) and PIC simulations (right panels). Parameters are: v b = 15v Te and ω pe = 30ω ce (upper panels), v b = 10v Te and ω pe = 30ω ce (middle panels) and v b = 15v Te and ω pe = 100ω ce (lower panels).
simulation shows a clear decrease of the growth rate at ω pe t ≈ 400, while in the PIC simulation the electric field rapidly reaches its peak value after which it decreases. We have previously identified this effect as a difference between the Vlasov and PIC codes in the representation of the electron distribution function, where the Vlasov code can represent the phase density tail of the Maxwellian distribution to lower values than the PIC code [11]. This leads to an earlier trapping of electrons and the formation of electron phase space holes in the Vlasov simulations compared to the PIC simulations, and to an earlier saturation in the Vlasov simulations than in the PIC simulations. The electric field reaches the largest amplitude and also shows the most rapid decrease for the largest values of the beam speed and the magnetic field, displayed in the upper 7 Institute of Physics ⌽ DEUTSCHE PHYSIKALISCHE GESELLSCHAFT panels of figure 2. Compared to this case, the peak electric field is almost halved for the case with somewhat slower beam speed in the middle panels. Here the electric field amplitude does not decrease as rapidly as in the higher beam speed case, but shows a complicated oscillatory behaviour. In the weak magnetic field case, displayed in the lower panels of figure 2, we observe that after the peak of the electric field, the electric field amplitude decreases relatively slowly compared to the case with the stronger magnetic field in the upper panels, and for the Vlasov simulation (left panels), there is a more distinct peak of the electric field in the strong magnetic field case than in the weak magnetic field case. This can be explained in that the fastest electrons, which are trapped early in the wave potential, start to experience surfing acceleration and are accelerated in the v y direction, after which they are decelerated and detrapped from the wave, so that they do not contribute to the saturation of the growing wave. In this way, the electric field can grow faster for a stronger magnetic field than in the weaker magnetic field, where the trapped electrons contribute to the saturation of the wave, clearly seen in the lower left panel of figure 2.
In order to investigate the influence of the magnetic field strength on the surfing mechanism, we have plotted the distribution function for electrons at different times for the case v b = 15v Te in figures 3-6. In these plots, we have made a projection of the three-dimensional distribution function in (x, p x , p y ) space down to two dimensions by integrating over the third dimension. Accordingly, we display the function F(x, v x ) = f dv y in the left panels and the function F(v x , v y ) = L 0 f dx in the right panels of figures 3-6. Figures 3 and 4 show the distribution functions at different times, corresponding to the electric field in the upper left and right panels of figure 2, respectively. The upper panels of figures 3 and 4 exhibit an electron distribution function at times shortly before the electric field has reached its peak value. Here (in the upper left panels), we see that unmistakable BGK modes have been formed in both the Vlasov and PIC simulations. In the upper right panels, we see that electrons have been accelerated mainly in the v x direction by the wave potential, while there is only a slight acceleration of electrons in the v y direction by the magnetic field. In the second right panel from the top of figures 3 and 4, we can see that electrons are being accelerated in the v y direction also by the magnetic field. In the lower panels of figures 3 and 4, we see that the BGK mode has disappeared (left panels) and that the distribution of electrons has reached a more or less isothermal distribution rotating in the (v x , v y ) space. The fastest electrons here have reached speeds of around 30v Te , which corresponds to the fastest particles initially trapped in the BGK mode seen in the upper panels. Thus, there has not been any efficient surfing acceleration in this case.
We now turn to the weaker magnetic field case, displayed in figures 5 and 6. The initial stage with the trapping of electrons is here similar to the stronger magnetic field case, as depicted in figures 3 and 4; therefore in the following, we present the phase-space distribution at somewhat later times. We see in the left panels of figures 5 and 6 that the BGK modes survive for a much longer time than in the stronger magnetic field case, displayed in figures 3 and 5. In the right panels of figures 5 and 6, we see that the acceleration cone in the v y direction continues to grow, and in the final stage, we see in the lower right panel of figure 6 that the fastest electrons have attained speeds of approximately 50v Te , which well exceeds the speed of the initially trapped electrons. At this stage, a population of electrons get detrapped from the BGK wave potential and start to rotate counterclockwise in the (v x , v y ) space. The surfing acceleration mechanism has thus been more efficient in the weak magnetic field case than in the strong magnetic field case, owing to the stability of the BGK modes.    The electron phase-space function (10-logarithmic colour scale) obtained from the PIC simulation, integrated over v y as a function of v x and x (left panels) and integrated over x as a function of v x and v y (right panels) at different times, for the weak magnetic field case ω pe /ω ce = 100. The beam speed is v b = 15v Te .