Particle Acceleration and Magnetic Field Amplification by Relativistic Shocks in Inhomogeneous Media

Particle acceleration and magnetic field amplification in relativistic shocks propagating in inhomogeneous media are investigated by three-dimensional magnetohydrodynamical (MHD) simulations and test-particle simulations. The MHD simulations show that the interaction between the relativistic shock and dense clumps amplifies the downstream magnetic field to the value expected from observations of the gamma-ray burst. The test-particle simulations in the electromagnetic field given by the MHD simulation show that particles are accelerated by the downstream turbulence and the relativistic shock. We provide the injection energy to the shock acceleration in this system. If the amplitude of upstream density fluctuations is sufficiently large, low-energy particles are initially accelerated to the injection energy by the downstream turbulence and then rapidly accelerated to higher energies by the relativistic shock. Therefore, the density fluctuation significantly affects particle acceleration in the relativistic shock.


INTRODUCTION
High-energy charged particles called cosmic rays (CRs) can be mainly separated into two categories, galactic and extra-galactic CRs.The most promising acceleration mechanism is first-order Fermi acceleration (diffusive shock acceleration) (Axford et al. 1977;Krymskii 1977;Blandford & Ostriker 1978;Bell 1978).In this process, particles are accelerated by a shock, crossing the shock front many times.The galactic CRs are thought to be accelerated by this mechanism in supernova remnants in our galaxy.On the other hand, the origin of extra-galactic CRs has long been a mystery.Naively, the extra-galactic CRs are expected to be accelerated by shocks in the same way as the galactic CRs.Since objects that can potentially produce the extra-galactic CRs mostly have relativistic outflows (Hillas 1984;Batista et al. 2019), relativistic shock waves have been intensively studied as an accelerator of the extra-galactic CR (Kirk & Schneider 1987;Begelman & Kirk 1990;Gallant & Achterberg 1999;Achterberg et al. 2001 Lemoine et al. 2006;Niemiec et al. 2006;Sironi et al. 2013;Kirk et al. 2022;Huang et al. 2023).
The relativistic shock acceleration needs strong downstream magnetic turbulence to return particles to the upstream region.For weakly magnetized relativistic shocks, the Weibel instability generates the magnetic turbulence, so that the shock acceleration works (Spitkovsky 2008).However, the Weibel-mediated shock cannot accelerate CRs with 10 20 eV within the dynamical time scale of astrophysical relativistic shocks because the length scale of the magnetic turbulence is much smaller than the gyroradius of accelerated particles (Sironi et al. 2013).Particles can be rapidly accelerated by the relativistic shock if the shock generates magnetic turbulence with a larger length scale (Lemoine & Revenu 2006;Kirk et al. 2022;Huang et al. 2023).In addition, the larger magnetic turbulence in the relativistic shock is required to explain bright prompt or afterglow emissions from gamma-ray bursts (GRBs) (Gruzinov & Waxman 1999;Huang et al. 2022).Nevertheless, the generation mechanism of the magnetic turbulence is still an open issue (Tomita & Ohira 2016;Tomita et al. 2019).

Morikawa et al.
In this work, we demonstrate, by solving the generation of magnetic turbulence, that the large-scale magnetic turbulence is generated in the relativistic shock, and then particles are rapidly accelerated in the relativistic shock.We consider a relativistic shock propagating in an inhomogeneous medium.For non-relativistic or mildly relativistic shocks, it is shown that the interaction between the shock and a shock-upstream denser region can drive the strong magnetic turbulence in the shock-downstream region (Giacalone & Jokipii 2007;Inoue et al. 2011;Mizuno et al. 2014;Fraschetti 2013;Ohira 2016a,b;Slavin et al. 2017;Romansky et al. 2020;Hu et al. 2022;Demidem et al. 2023;Bresci et al. 2023;Fulat et al. 2023).Although there are few studies about a generation of downstream magnetic turbulence in the case of the relativistic shock (Sironi & Goodman 2007;Goodman & Macfadyen 2008;Tomita et al. 2022), none of them demonstrate the generation of the magnetic field turbulence that can accelerate CRs rapidly and explain emissions from GRBs.First of all, we perform three-dimensional relativistic magnetohydrodynamical (MHD) simulations to obtain the magnetic field structure in the downstream region of a relativistic shock propagating to an inhomogeneous medium.Then, we perform test-particle simulations, showing that particles are rapidly accelerated in the magnetic field obtained in the relativistic MHD simulations.Hereafter, we use the notation X ,U , X ,S , and X ,D to indicate the quantity, X, measured in the upstream, shock, and downstream rest frames, respectively.

Relativistic MHD simulations
In order for particle acceleration to work in the relativistic shock, particles have to cross magnetic field lines.To avoid the unphysical suppression of the particle diffusion perpendicular to the magnetic field line, we need a magnetic field structure in the realistic threedimensional space for the test-particle simulation (Giacalone & Jokipii 1994;Jokipii et al. 1993).We used the special relativistic three-dimensional MHD simulation code, SRCANS+ (Matsumoto et al. 2019), where the approximate Riemann solver is the HLL method, the MUSCL method in 2nd order is adopted for the spatial higher-order method, the limiter is the van Leer limiter, and the error of ∇ • B is removed by using the 9-wave method (Dedner et al. 2002;Mignone et al. 2010).We adopted the equation of state proposed in Mignone & McKinney (2007) for primitive recovery.The simulation box size is (L x /∆x, L y /∆y, L z /∆z) = (8000,200,200), where L s and ∆s = 0.1c∆t MHD are the box size and grid size in the s directions (s = x, y, z) and ∆t MHD is the time step.The x-direction is set to be perpendicular to the mean shock front.The boundary condition in the xdirection is the open boundary and the periodic boundary condition is imposed for the y-and z-directions.We set an initial condition so that the shock front averaged over the y-z plane appears stationary, that is, the simulation frame is the mean shock rest frame.
In the simulation frame, the upstream bulk Lorentz factor is Γ up,S = 5.2 and the upstream density structure is given as follows, where ρ 0,S is the constant density and the individual clumps are represented by where x ci , y ci , z ci , d, and L cl = 0.1L y are the center of the individual clump, distance from the center, and the half-width of the clumps, respectively.ρ cl,S is the density at the half width.We randomly put the clumps in the upstream rest frame.Then, the Lorentz transformation was performed to obtain the center of each clump in the simulation frame.The shape of the clumps is a sphere in the upstream rest frame.The volume filling factor of the density clumps is set to N cl V cl,S /L x L y L z = 0.2512 in this work, where N cl and V cl,S = 4πL 3 cl /3Γ up,S are the number of clumps in the simulation box and the volume of each clump, respectively.Initially, the magnetic field in the upstream rest frame and gas pressure are assumed to be uniform, B up,U and p up , and the magnetic field is pointing in the ydirection.The upstream magnetization parameter and plasma beta are set to σ = B 2 up,U /4πρ 0,U c 2 = 3.7 × 10 −6 and β = 8πp up /B 2 up,U = 5.4 × 10 3 , where c is the speed of light.In this work, we perform two MHD simulations for ρ cl,S /ρ 0,S = 2 and 10.

Test-particle simulations
Since the shock front is almost stationary in our simulation frame and the downstream turbulent velocity is about c ∼ 0.1, we use the snapshot data of the electromagnetic field obtained by the MHD simulations at t = 3000L cl /c to solve the equation of motion for charged particles.This stationary velocity field u ,S can accelerate particles through the electric field, E ,S = −(u ,S /c)×B ,S .In order to solve the equation of motion for the charged particles, we adopt the Buneman-Boris method (Birdsall & Langdon 1991) with the time step of ∆t TP = 0.02π∆x/c.10 6 simulation particles are uniformly injected in the upstream region.The initial momentum distribution is isotropic with the Lorentz factor of γ p,U = 10 in the upstream rest frame.We follow particle orbits until the particles reach the downstream boundary at x = 300 L cl .Because the shock front is located at x ∼ 240 L cl in our MHD simulation, the downstream size in the x-direction is L down ∼ 60 L cl .
To understand the effect of clump size on particle acceleration, we perform two test-particle simulations (R up,U /L cl = 26 and 0.26) for each of the two MHD simulations, where R up,U is the initial gyroradius in the upstream region in the upstream rest frame.

RESULTS
Figure 1 shows the MHD simulation results at t = 3000 L cl /c.The top two panels are the two-dimensional distribution of magnetic field strength in the z = 0 plane for ρ cl,S /ρ 0,S = 2 and 10, and the bottom two panels are the one-dimensional distribution of magnetic energy density and electric field strength averaged over the yz plane.The solid and dashed lines correspond to the cases of ρ cl,S /ρ 0,S = 2 and 10.The left region of the shock front is the upstream region.The magnetic field is amplified by the shock compression just behind the shock front and by the turbulence in more distant regions.The coherent length scale of the downstream turbulence is about L cl .Interestingly, the saturation value of ϵ B = B 2 ,S /4πΓ up,S ⟨ρ up,S ⟩c 2 is on the order of 10 −4 , which does not strongly depends on the density of clumps, ρ cl,S /ρ 0,S .The value of ϵ B ∼ 4×10 −4 is a highly anticipated result in the GRB community (Gruzinov & Waxman 1999;Huang et al. 2022).
Figure 2 is the energy spectra of particles in the test-particle simulation when they arrive at the downstream boundary.The particle energy is normalized by E cl = qB up,U L cl , where q is the particle charge.The maximum energy in this system is limited by the distance between the shock front and the downstream boundary, E max,S ∼ 10 3 E cl .The results with and without density variations are shown in the solid and dotted histograms, and the results for R up,U /L cl = 0.26 and 26 are shown in the red and blue histograms, respectively.For R up,U /L cl = 26, some particles are accelerated but the peak of the energy spectrum does not change significantly compared to the case without density fluctuations.On the other hand, for R up,U /L cl = 0.26, the position of the spectral peak significantly moves to the higher energy region.Comparing the top (ρ cl,S /ρ 0,S = 2) and bottom (ρ cl,S /ρ 0,S = 10) panels, the particles are accelerated more when the amplitude of density fluctuations is larger.
To understand the acceleration mechanism, we plot trajectories of the highest-energy particles for four runs in Figure 3, where the horizontal and vertical axes show the particles' position in the mean shock normal direction and the particle energy.The gray region represents the shock front.The solid and dotted lines are the trajectories for R up,U /L cl = 0.26 and = 26, respectively.For both cases of ρ cl,S /ρ 0,S = 2 and 10 with R up,U /L cl = 0.26, the particles do not go back to the upstream region, but are accelerated in the downstream region in the range of E ,S ≲ 100E cl .This acceleration by the downstream turbulence is theoretically expected in early studies (Ohira 2013;Asano & Terasawa 2015;Pohl et al. 2015;Yokoyama & Ohira 2020).The particle with the small gyroradius in the case of ρ cl,S /ρ 0,S = 10 is effi- ciently accelerated just behind the shock front compared with one for ρ cl,S /ρ 0,S = 2 (solid lines).The fourth panel in Figure 1 shows the electric field strength associated with the downstream turbulence, δE ,S /E up,S , where δE ,S = ⟨|−(δu ,S /c)×B ,S |⟩ yz , E up,S = (u up,S /c)B up,S ∼ B up,S , and δu ,S = u ,S − ⟨u ,S ⟩ yz is the fluid velocity associated with turbulence.⟨Q⟩ yz is a quantity Q averaged over the y-z plane.u up,S is the upstream flow velocity.The turbulent electric field for ρ cl,S /ρ 0,S = 10 is stronger than one for ρ cl,S /ρ 0,S = 2.The interaction with a denser clump generates stronger turbulence, resulting in faster magnetic field amplification, so that particles with the small gyroradius for ρ cl,S /ρ 0,S = 10 are rapidly accelerated by the downstream turbulence.For both cases of ρ cl,S /ρ 0,S = 2 and 10 with R up,U /L cl = 26 and the case of ρ cl,S /ρ 0,S = 2 with R up,U /L cl = 0.26 in the range of E ,S ≳ 100 E cl , particles are accelerated by crossing the shock front many times.Our simulations show that the relativistic shock acceleration works if there are particles with a sufficiently large gyroradius in the vicinity of the shock front.This can be understood as follows.It takes the eddy turnover time, T eddy,S ∼ L cl /δu ,S , for the shocked magnetic field to be disturbed.Hence, the magnetic field is disturbed at the distance of T eddy,S × u down,S away from the shock front, but closer to the shock front than that, the magnetic field is almost uniform and perpendicular to the shock normal, where u down,S = u 1,S /r ≈ c/r is the downstream flow velocity in the shock rest frame, and r is the shock compression ratio.To go back to the upstream from the downstream, the downstream gyroradius, E ,S /qrB up,S , has to be larger than T eddy,S × u 2,S , which gives the injection energy for the shock acceleration, (4) Since our MHD simulations show c/δu ,S ∼ 10 in the downstream region, E inj,S ∼ 50 E cl , which is consistent with results in our test-particle simulations.
For the case of ρ cl,S /ρ 0,S = 2 with the initial small gyroradius, although some particles are accelerated to more than E inj,S in the far downstream region, they can not go back to the upstream region because they have already passed the diffusion length.In contrast, for ρ cl,S /ρ 0,S = 10 with the initial small gyroradius, the particle is quickly accelerated to E inj,S around the shock front, so that the particle can experience the relativistic shock acceleration.To show that most particles are similarly accelerated, first, we extract parts of the trajectories of the top 10 3 most energetic particles moving from upstream to downstream and back upstream again.Then, we sample the incident energy from upstream to downstream, E ini,D , and the downstream energy gain, ∆E ,D , in the downstream rest frame.For the shock acceleration without turbulent acceleration, ∆E ,D = 0 and the mean energy gain during the one cycle is about E ini,D .Figure 4 shows ∆E ,D /E ini,D ≳ 1 for Therefore, for ρ cl,S /ρ 0,S = 10 with the initial small gyroradius, most particles are accelerated first by the downstream turbulence and then by the relativistic shock.The injection energy to the shock acceleration is consistent with Equation ( 4).This injection process to the shock acceleration has never been observed in previous studies.The detailed mechanism of the downstream acceleration will be addressed in future work.
Finally, we show the mean evolution of the energy for the top 1000 energetic particles in Figure 5, where the solid and dashed curves are the simulation results and theoretical curve of the Bohm rate in the upstream magnetic field strength, (dE/dt) ,S = qB up,S c.Except for the top left panel, the simulation results are very close to the theoretical curve in the late phase when particles are accelerated by the relativistic shock.This is the first simulation that reproduces the ideal fast acceleration of CRs in the relativistic shock by calculating the magnetic field amplification and particle trajectory.In the relativistic shock acceleration, particles double their energy for every round trip of a relativistic shock (Achterberg et al. 2001).The downstream residence time is negligible because the downstream magnetic field is strongly amplified in our simulation.In this case, the relativistic shock acceleration rate is given by the Bohm rate in the upstream magnetic field strength (Kamijima et al. 2020).

DISCUSSION AND CONCLUSION
The evolution of the magnetic field amplification and downstream turbulence depends on the resolution of MHD simulations (e.g.Ji et al. 2016).Our MHD simulations did not have enough resolution to converge all the results.Nevertheless, we can theoretically explain our results concerning the injection to acceleration by relativistic shocks with Equation (4).Therefore, most results on particle acceleration are qualitatively robust although dependency of the clump density and initial gyroradius may change quantitatively.
In Figure 5, we showed that the ideal rapid acceleration works in the relativistic shock.If the downstream gyroradius is larger than the coherence length scale of  Mean downstream energy gain in the downstream rest frame for the top 10 3 energetic particles in the case of ρ cl,S /ρ0,S = 10 and Rup,U/L cl = 0.26.The horizontal axis is the incident energy from upstream to downstream, Eini,D.magnetic field turbulence, L cl , the scattering efficiency is smaller than the Bohm rate in the downstream magnetic field strength.In this case, the downstream diffusion coefficient, κ, is given by where r down,S is the downstream gyroradius.Then, the downstream residence time of accelerated particles is roughly given by t d,S ∼ κ/c 2 (e.g.Sironi et al. 2013;Lemoine et al. 2019;Vanthieghem et al. 2020).On the other hand, the upstream residence time of accelerated particles is given by t u,S ∼ Ω −1 up,S (Achterberg et al. 2001).If the upstream residence time is longer than the downstream residence time, the acceleration time scale is given by the upstream residence time, that is, the acceleration rate is given by the Bohm rate in the upstream magnetic field strength (Kamijima et al. 2020).From the condition, t u,S > t d,S , the energy range in which the ideal acceleration at the Bohm rate works is given by where Bdown,S ∼ 20 B up,S is the average downstream magnetic field in our MHD simulation.If the particle energy is larger than the above range, the acceleration rate is slower than the Bohm rate in the upstream magnetic field.We showed in Figure 3 that some particles are accelerated by the relativistic shock, but Figure 2 does not show a clear power-law spectrum.This is probably due to the limitation of the acceleration region.The maximum energy is at least limited by the finite size of the downstream region, L down , and given by the condition, κ/c ∼ L down .Then, from Equation ( 5), the sizelimited maximum energy is estimated to E max,S /E cl ∼ 800 (L down /60L cl ) 1/2 (Γ up,S /5)( Bdown,S /20B up,S ).We cannot expect a clear power-law spectrum above this energy, which is consistent with our results.If the upstream density clumps have different larger length scales, higher energy particles would be efficiently confined in the downstream region.A clear power-law energy spectrum could then be observed in the test particle simulation.To confirm this point, we need to perform simulations with a larger simulation box.In addition, we have to investigate dependences of physical parameters (σ, β, Γ up,S ), and other realistic density structures (Hu et al. 2022).
Our simulations showed that particles with the energy of 0.26 E cl are accelerated in the downstream region, where E cl ∼ 10 12 eV (B up,U /1 G) L cl /10 10 cm .The clump size of 10 10 cm in the stellar wind is suggested by several studies (Sironi & Goodman 2007;Moens et al. 2022;Chené et al. 2020).Since the energy scale of E cl is much higher than the thermal energy of upstream particles, there is a huge gap between thermal particles and our simulation particles.The time scale of the downstream acceleration is an increasing function of the particle energy (see Figure 5).We naively expect a similar energy dependence even in a smaller energy range, but our simulations cannot investigate whether thermal particles are accelerated mainly by the same mechanism or not.On the other hand, it has been shown by particle-in-cell simulations that the Weibel-mediated collisionless shock can accelerate thermal particles to higher energies (Spitkovsky 2008;Sironi et al. 2013).In addition, the Weibel turbulence would suppress particle streaming along the magnetic field line.Particle streaming makes the high-density region smoother, resulting in weaker downstream turbulence (Tomita et al. 2022).The Weibel turbulence, which cannot be resolved by MHD simulations, therefore plays an important role in bridging the gap between the thermal particles and our test-particle simulations, and in keeping the amplitude of density fluctuations large.
In summary, using three-dimensional relativistic MHD simulations and test-particle simulations, we have shown the following.
1. Relativistic shocks propagating to inhomogeneous media amplify the downstream magnetic field to ϵ B ∼ 4 × 10 −4 .2. Particles are rapidly accelerated by the relativistic shock and downstream turbulence in this system.3.If the amplitude of the upstream inhomogeneity is sufficiently large, low energy particles are initially accelerated by the downstream turbulence and then accelerated by the relativistic shock.4. The injection energy to the shock acceleration in this system is given by Equation (4).
These results have long been awaited to explain emissions from relativistic shocks such as GRBs and active galactic nuclei.Calculating electromagnetic waves, neutrinos, and CRs from this system and comparing them with observations will help us to understand the environment surrounding the high-energy objects, leading to a better understanding of the origin of extra-galactic CRs and enigmatic high-energy objects.Moreover, this work can be applied to non-relativistic perpendicular shocks realized in supernova remnants of massive stars.The injection to the shock acceleration at the perpendicular shock would be enhanced by the upstream density inhomogeneity.

Figure 1 .
Figure1.MHD simulation results at t = 3000 L cl /c.The top two panels show the two-dimensional distribution of the magnetic field strength in the z = 0 plane.The third and last panels show one-dimensional distributions of magnetic energy density and electric field strength associated with the downstream turbulence, respectively.The bottom two panels show quantities averaged over the y-z plane.

Figure 2 .
Figure 2. Energy spectra of test-particle simulations.The top and bottom panels are for ρ cl,S /ρ0,S = 2 and 10, respectively.The red and blue solid histograms are the energy spectra for Rup,U/L cl = 0.26 and = 26, respectively.The dotted histograms show the energy spectra without density fluctuations.

Figure 3 .
Figure 3.The representative particle trajectory for all the simulations.The top and bottom panels are for ρ cl,S /ρ0,S = 2 and 10, respectively.The gray region represents the shock front.
Figure 4.Mean downstream energy gain in the downstream rest frame for the top 10 3 energetic particles in the case of ρ cl,S /ρ0,S = 10 and Rup,U/L cl = 0.26.The horizontal axis is the incident energy from upstream to downstream, Eini,D.

Figure 5 .
Figure5.Mean time evolution of the particle energy for the four runs.t = 0 is defined by the time when each particle crosses the shock front for the first time.The solid curves are the simulation results averaged for the top 1000 particles.The dashed lines are the theoretical line of the Bohm rate in the upstream magnetic field strength, where the initial energy was ignored.