Self-consistent Monte Carlo model of particle acceleration by relativistic shocks

We present a self-consistent Monte Carlo model of particle acceleration by relativistic shock waves. The model includes the magnetic field amplification in the shock upstream by cosmic ray driven plasma instabilities. The parameters of the Monte Carlo model are obtained based on PIC calculations. We present the spectra of accelerated particles simulated in the frame of the model.


Introduction
The processes of particle acceleration and the magnetic field amplification by shocks that occur when a powerful relativistic current collides with the environment are appropriate for many astrophysical objects. Such flows are formed, for example, in the propagation of a pulsar wind, the explosion of relativistic supernovae, and the formation of gamma-ray bursts. These objects are powerful sources of non-thermal radiation. Relativistic shock waves in astrophysical objects are often collisionless. The particles are accelerated to high energies by such shock waves according to the first-order Fermi mechanism, when they are scattered by magnetic fluctuations. The distribution function of accelerated particles has a significant anisotropy in the upstream of a relativistic collisionless shock wave. This anisotropy leads to the development of plasma instabilities that amplify magnetic fields [1][2][3]. In the same way, the pressure of accelerated particles modifies the prefront of the shock wave [4][5][6][7][8][9]. Thus, the problem of particle acceleration by a relativistic collisionless shock wave is essentially nonlinear, since the propagation of particles is determined by magnetic fields, and the amplification of fields is a function of the distribution of accelerated particles.
We develop a stationary self-consistent Monte Carlo model of particle acceleration by relativistic shock waves, which takes into account the amplification of magnetic fields by plasma instabilities. Model allows us to construct the particle distribution function, plasma flow profiles, and magnetic fluctuation spectra in an iterative scheme that implements the laws of conservation of momentum and energy near the shock wave front. PIC (particle in cell) calculations [10][11][12][13][14][15][16][17] based on the first principles allow us to solve this problem more accurately. PIC calculations have very significant limitations on the size of the calculation space due to the limited computing power. We compare the results of PIC calculations and the results of our Monte Carlo simulation,

Description of the Monte Carlo model
The Monte Carlo model of particle acceleration by an unmagnetized relativistic stationary planeparallel shock presented here is self-consistent. The model in an iterative scheme constructs the shock flow structure which obeys the momentum and energy conservation laws. The iteration scheme in the relativistic Monte Carlo model is similar to that used in the modeling of nonrelativistic shocks with the magnetic field amplification developed earlier [18][19][20][21]. The plasma in the model is divided between the background and the accelerated particles. The accelerated particles are those that have crossed the shock wave front at least once from the downstream to the upstream. In most non-relativistic calculations, the background plasma is described hydrodynamically in most of the upstream, which makes it possible to speed up greatly the operation of the numerical scheme. In the calculations of this paper, the entire plasma is described as particles in the entire computational domain.
The pressure gradient of the accelerated particles modifies the upstream momentum fluxes of the background plasma. The distribution function of accelerated particles has a significant anisotropy in the upstream, which leads to the development of plasma instabilities. In this paper, we describe the magnetic field amplification in the upstream by the following equation where u w (x) and γ w (x) are velocity and Lorentz factor of the scattering centre rest frame relative to the shock rest frame accordingly, W (x) and G (x) are turbulence energy density and the growth rate of plasma instabilities in the rest frame of scattering centers accordingly. The model assumes that the magnetic field amplification occurs in a narrow range of wave numbers. The second term of the left side of the equation determines the adiabatic amplification of magnetic fluctuations with a change in the velocity of the scattering centers. Thus, the parameter α w determines the adiabatic amplification of the magnetic field, in the calculations α w = 2. The x coordinate is calculated from the shock front, along the normal to the front. Negative values x correspond to the upstream, positive values of x correspond to the downstream. In the far upstream, there is a free-escape boundary located at x = x F EB . Any accelerated particle which is crossing the boundary leaves the simulation domain. Near the free-escape boundary, the Lorentz factors of the background plasma flow and the rest frame of scattering centers coincide and are equal to the Lorentz factor of the shock γ sh . The scattering of particles occurs in the rest system of the scattering centers and in this system they are isotropic and elastic. The free path length of particles in the scattering center rest frame in this paper where L c (x) is the turbulence correlation length, which is determined by the range of wave numbers in which the magnetic fields are amplified, the particle gyroradius is where c is the speed of light, e is the particle charge, p is the momentum modulus of a particle in the scattering center rest frame, International  The growth rate of plasma instabilities G (x) and the turbulence correlation length L c (x) were obtained by comparing the simulation results with the results of PIC calculations [15][16][17]. The blue curves correspond to γ sh = 100. r g0 = m e c 2 eB 0 ≈ 5.7 · 10 8 cm.

Results
Based on our model in this paper we perform a simulation of particle acceleration and magnetic field amplification by shocks with Lorentz factors γ sh = 10 and γ sh = 100. The electron-positron plasma is considered. The value of the other calculation parameters is the same: x F EB = 10 3 r g0 , where r g0 = m e c 2 eB 0 , m e is the electron mass, B 0 = 3µG is the root-mean-square value of the magnetic field in the far upstream in the background plasma rest system, n 0 = 4 · 10 −3 cm −3 is the total concentration of electrons and positrons in the far upstream in the background plasma rest system. Figure 1 is given to demonstrate the convergence of the iterative scheme of our Monte Carlo model. This figure shows the energy flux near the shock for two cases: a converged nonlinear solution of the problem and an unmodified calculation. The unmodified calculation does not include self-consistent magnetic field amplification and modification of the upstream by the momentum flow of accelerated particles. As can be seen from Figure 1, in the unmodified case, the energy flow near the shock wave is not conserved.
The results of the self-consistent simulation are shown in the Figures 2 and 3. The particle distribution functions in the shock rest system are shown on Figure 2. As can be seen from this Figure for the values of the particle momentum ≤ γ sh m e c, the particle distribution function has the form of a thermal spectrum. For the values of the particle momentum > γ sh m e c, the particle distribution function has the form of a power spectrum. The spectrum slope in the power region is close to -4.22. This slope of the spectrum is obtained both in calculations for particle acceleration by relativistic shocks and in the test-particle regime. Since the spatial domain of calculations in this paper is an order of magnitude larger than that achievable in PIC modeling the value of the maximum particle momentum is several times greater than in PIC modeling for similar parameters of the problem.
The simulated magnetic field profiles are shown on Figure 3. A sharp increase in the turbulent magnetic field occurs near the free-escape boundary due to Weibel and filamentation plasma instabilities [14,17,[22][23][24][25][26]. Closer to the shock front, the magnetic field amplification occurs mainly adiabatically due to the second term in the left part of the equation (1) field amplification effect is essential to explain the observed non-thermal synchrotron emission from the shocks driven by relativistic outflows.