Particle acceleration at magnetized, relativistic turbulent shock fronts

The efficiency of particle acceleration at shock waves in relativistic, magnetized astrophysical outflows is a debated topic with far-reaching implications. Here, for the first time, we study the impact of turbulence in the pre-shock plasma. Our simulations demonstrate that, for a mildly relativistic, magnetized pair shock (Lorentz factor $\gamma_{\rm sh} \simeq 2.7$, magnetization level $\sigma \simeq 0.01$), strong turbulence can revive particle acceleration in a superluminal configuration that otherwise prohibits it. Depending on the initial plasma temperature and magnetization, stochastic-shock-drift or diffusive-type acceleration governs particle energization, producing powerlaw spectra $\mathrm{d}N/\mathrm{d}\gamma \propto \gamma^{-s}$ with $s \sim 2.5-3.5$. At larger magnetization levels, stochastic acceleration within the pre-shock turbulence becomes competitive and can even take over shock acceleration.

One limitation of previous numerical studies based on fully PIC simulations, though, is to consider laminar inflow conditions, i.e., nonturbulent, homogeneous background plasmas of uniform magnetization. Notable exceptions are, to our knowledge, Ref. [21], which shows that the interaction of a shock with a monochromatic linear eigenmode of the upstream plasma leads to particle acceleration in the resultant downstream turbulence, Refs. [22,23] which examine the impact of a superposition of magnetostatic plane waves seeded upstream, 1 We define the magnetization level as σ ≡ B 2 /4π , in terms of mean-squared magnetic field B 2 and the energy density = n γ mec 2 as measured in the simulation (downstream) frame; n represents the total apparent density. When measured in the comoving turbulence frame, the resulting value of σ may be a factor ∼ 2 larger. and finally Ref. [24] that considers the influence of an anisotropic transverse upstream magnetic profile. The presence of a strong turbulence upstream of a fast shock may change the picture in various ways: it may preaccelerate the plasma particles via a stochastic Fermi process [7][8][9][10]25], just as it may corrugate the shock front so that the turbulence does not transform trivially through the shock [26][27][28][29][30][31][32], possibly unlocking particles from the field lines and enabling their acceleration. In light of these considerations, the paradigm of inefficient relativistic magnetized shocks as particle accelerators needs to be revisited in the likely common case of turbulent environments.
To this goal, we here report on the first PIC simulations of relativistic shocks propagating in turbulent, magnetized pair plasmas. We demonstrate that, despite a substantial magnetization (σ ∼ 0.01), shock acceleration is manifest, and that the particle spectrum develops a powerlaw tail extending in time, which is absent without turbulence. This drastic change involves a significant (δB > B 0 ), but not too strong upstream turbulence, otherwise its own contribution to particle acceleration can supersede that of the shock. The discussion is organized as follows: we detail the simulation technique in Sec. II, then investigate the acceleration processes at play in Sec. III, before concluding and summarizing our results in Sec. IV.

II. NUMERICAL METHOD
We perform such simulations using the fully electromagnetic and relativistic calder code [33][34][35]. Turbulence is excited close to the right-hand side of the domain, in a pair plasma continually injected along −x at a relativistic velocity v ∞ = −0.87 c (Lorentz factor γ ∞ = 2). The flow is left to propagate across the domain until the turbulence hits its left-hand side. Switching the local boundary condition from open to reflective at that time arXiv:2303.11394v2 [astro-ph.HE] 3 Jul 2023 triggers a rightward propagating shock wave [5], which sweeps the incoming turbulent plasma. The simulation frame thus corresponds to the downstream rest frame of the shock. Due to physical constraints discussed thereafter, we restrict ourselves to 2D3V geometry (2D in space, 3D in momentum). A uniform magnetic guide field B 0 is applied along the (out-of-plane) z direction with corresponding magnetization level σ 0 .
The turbulence is driven in the rest frame of the drifting plasma, in a finite region covering a few stirring length scales c . Elsewhere, the energy injection in the system is halted and the turbulence is let to develop freely, thereby initiating the cascade before it interacts with the shock. Following [7], we aim at exciting turbulence on a scale c as large as possible compared to the kinetic scale c/ω p , in order to simulate an inertial range under near-magnetohydrodynamics conditions where the turbulence cascades down to the dissipative range. Those simulations are demanding, because the interaction between the turbulence and the shock must be followed over a long enough timespan, the transverse dimension must accommodate 1 − 2 c , and also because the need to stir turbulence in the plasma rest frame brings in further constraints due to time dilation effects.
A relativistically hot initial plasma as in S2 and S3 could describe internal shocks inside a strongly turbulent jet. We have run ancillary simulations, in particular S2a, similar to S2 albeit deprived of turbulence, S2b which retains open boundary conditions and thus models drifting turbulence without a shock, and finally S4, for which k B T /m e c 2 = 0.1 as in S1, but with larger σ 0 and σ δB , as in S2. Figure 1 displays (from top to bottom) the spatial distributions of plasma positron density, mean Lorentz factor and magnetization level at the final simulation time t 12 400 ω −1 p for simulation S2. The rightward-moving shock front has then reached x 2 400 c/ω p . Once swept up by the shock, the plasma is compressed by a factor of 3.5 and the mean kinetic energy per particle slightly increases, in good agreement with the shock-crossing conditions [37] [Figs. 1(a,b)]. Magnetic fluctuations are at their highest near the right boundary where turbulence is continually excited [ Fig. 1(c)]. The ∼ 2 000 c/ω p distance between the shock and the boundary of forced turbulence is then just below the minimum distance cτ nl needed for nonlinear evolution of the turbulence. This guarantees that the shock-turbulence interaction is not affected by the stirring procedure in the right part of the domain.

III. RESULTS
As shown in Fig. 1(d), the turbulence profile reaches an approximately steady state by the time it encounters the shock. We have checked that the spatial power spectrum of magnetic fluctuations in the transverse y direction, which extends over three orders of magnitude, shows a general scaling ∝ k −5/3 y at large scales k y 15 −1 c , and a steeper behavior at kinetic scales, consistent with previous PIC studies of nondrifting decaying turbulence [8].
Further upstream, corresponding to an earlier stage in the turbulence evolution, stronger fluctuations are observed, as expected. The eddies are compressed when transiting across the shock and continue interacting until the turbulence eventually relaxes further downstream. This general picture resembles that observed in MHD simulations of the interaction of a monochromatic, linear plasma eigenmode with a relativistic shock front [30].
The shock moves at velocity 0.4 c in the simulation (downstream) frame (0.36 c is predicted by the shockcrossing conditions [37]), and consequently at v sh 0.93 c in the upstream frame, corresponding to a shock Lorentz factor γ sh 2.7. Ahead of the shock, the transversely averaged magnetic field strength is δB/(m e ω p c/e) 0.8, so that particles with Lorentz factor γ 100 − 300 have a gyroradius r g 120 − 360 c/ω p . Figure 2 plots the time evolution of the particle energy spectra γ 2 dN/dγ (per log-interval of energy) in each of our simulations, as integrated over a moving window centered on the shock front position x sh (located from the plasma density map) and with a 200 c/ω p half-width along x. The peaks and widths of the spectra -for those simulations with shock -are consistent with shock dis- x + δB 2 y (normalized to 4πn0mec 2 ) for simulation S2 at t 12 400 ω −1 p . In panel (d), the longitudinal profile of the turbulent magnetization σ δB = (δB 2 x + δB 2 y )/(4π α nα γ αmec 2 ), with nα the total (apparent) density of the plasma, averaged over the transverse dimension (y). The left and right vertical lines locate the shock front and the boundary of forced turbulence.
FIG. 2. Particle energy spectra γ 2 dN/dγ at different times (from light to dark solid blue) in simulations S1, S2 and S3, from top to bottom. Middle panel: in dotted line, spectrum from S2a, i.e., a shock interacting with a non-turbulent plasma, and in dashed line, spectrum from S2b, i.e., a drifting turbulence without a shock, both in conditions otherwise similar to S2. In the lower panel, the light orange band delineates the range of spectra measured in a simulation similar to S3 albeit without a shock, as extracted at various places and times. sipation, as predicted by the shock-crossing conditions. Remarkably, a suprathermal tail develops in all cases, with approximate powerlaw index s 2.5 → 3.5 (as defined through dN/dγ ∝ γ −s ), providing manifest evidence of particle acceleration. In S1 and S2, and unlike in S3, the maximal energy is seen to increase with time. The middle panel of Fig. 2 presents additional spectra from the turbulence-free simulation S2a, thus where the shock forms immediately in the external field B 0 (dotted line), and from shock-free simulation S2b, which contains only drifting turbulence (dashed line). Clearly, the suprathermal tail only arises when the shock interacts with the turbulence. The absence of particle acceleration in S2a can be attributed to the background magnetization level σ ≡ σ 0 10 −3 , large enough to inhibit Fermi cycles around the shock. In S2b, the magnetic fluctuations are too slow to accelerate particles, as the characteristic acceleration timescale t acc ∼ γ ∞ c c /v 2 A ∼ 10 4 ω −1 p indeed exceeds the time needed for the plasma to cross the domain. Note that the spectrum of S2b coincides with that of S2 at t = 5940 ω −1 p , because the plasma has then just hit the reflective wall, and the shock has not formed yet.

A. Particle acceleration
Standard theory depicts acceleration at a magnetized shock front as the result of diffusion back and forth across the shock, or of shock-drift motion along the mean advected electric field [38]. In the relativistic limit, shock acceleration is ineffective [12,13] unless intense turbulence can unlock particles off the field lines [17], and thereby trigger diffusive-type acceleration [39][40][41], a form of shock-drift process [42,43], or a combination of both, meaning orbits in the regular field upstream of the shock, diffusive orbits downstream [15]. Let us stress here that what we mean by "diffusive-type" does not correspond to the standard "spatial diffusion" at play in subrelativistic shocks, but rather to "diffusion in pitch-angle", which ensures that particles can return to the shock. It is indeed known that in relativistic shocks, spatial diffusion does not have time to set in properly, because of the high advection velocity downstream of the shock and because the upstream particles can be caught back by the shock front just by barely changing their propagation direction [39,41]. Similarly, the shock-drift type process that we will refer to in the following is sustained by pitch-angle scattering in the turbulence, which allows particles to remain close to the shock surface. It could be termed "stochastic shock-drift" in analogy with Ref. [43], although those authors considered a subluminal configuration, not a superluminal one as in the present case.
To probe the acceleration process at work, we have tracked a large number ∼ O(10 5 ) of particles sampled in various (initial) energy intervals (see Appendix B). Figure 3 shows the trajectories and energy histories of four particles in S2, representative of the population able to circulate around the shock for an extended period of time. The Lorentz factor of some particles (e.g. orange and cyan in that figure) undergoes sizable oscillations before they encounter the shock; this results from their gyromotion along the fast-moving magnetic field lines [44], not from acceleration per se. Notwithstanding this effect, the energization of the particles traveling in the vicinity of the shock is evident.
We discriminate the acceleration processes in simulations S1 and S2 using the following argument. In a shock-drift process, the energy gained by a particle of velocity v equals the amount of work performed by the mean motional electric field E 0 = v ∞ B 0ŷ , i.e., W (E 0 ) = q dt v · E 0 , whereas for diffusive-type acceleration, the energy gain is rather related to the work W (δE z ) = q dt v z δE z performed by the motional turbulent electric field. The latter is mostly directed alonĝ z because the plasma, which flows along −x, carries essentially (δB x , δB y ) magnetic fluctuations; hence δE z −v ∞ δB y . For each tracked particle with initial Lorentz factor at the onset of the powerlaw tail, i.e. γ ≥ 20 (in S1) and γ ≥ 200 (S2), we have thus recorded W (E 0 ) and W (δE z ) during the time interval ∆t sh between the first and last encounters of the particle with the shock front. . In these plots, the dashed red line indicates the expected level of contribution from either shock-drift or diffusive-type acceleration. Figures 4 (c) and (d) reveal that a shock-drift process sustained by particle scattering along the shock front nicely accounts for particle energization in S2, whereas in S1, acceleration appears dominated by diffusive-type acceleration.
This general picture is further supported by the angular distribution of the suprathermal particle momenta (see Appendix B). In S2, this angular map presents a clear asymmetry between positrons and electrons, roughly polarized along E 0 , such that positrons (resp. electrons) appear to drift with negative (resp. positive) p y , as expected for E 0 < 0. By contrast, the angular map is significantly more isotropic in S1, as expected if diffusive-type acceleration dominates. Additionally, simulation S2 displays a net linear correlation between ∆γ and ∆t sh , from which one can infer an acceleration timescale, t acc ≡ | ∆γ/γ | −1 ∆t sh ∼ 2 − 3 p/(eE 0 ), again consistent with that expected for particles drifting along E 0 at mildly relativistic speeds.
The difference in spectral index observed between sim-ulations S1 (s ∼ 2.5) and S2 (s ∼ 3.5) likely results from the distinct acceleration mechanisms at play. In particular, the spectral index for shock-drift acceleration -at least, in subrelativistic shocks -depends sensitively on the shock speed and on the ratio of the scattering frequency of particles to their gyrofrequency in the mean field [42]. At large scattering frequencies, the index approaches the canonical value of 2, associated with diffusive shock acceleration, whereas at small scattering frequencies, the spectrum steepens significantly, encompassing the value measured in S2.
Regarding the different acceleration processes in S1 and S2, we observe that the r g vs c ordering, which controls the scattering rate of particles, varies between those two simulations because of different initial temperatures and magnetizations: r g / c ∼ 0.1 at the onset of the powerlaw tail for S1, whereas r g / c ∼ 1 for S2. A detailed study of the ancillary simulation S4, which shares the same initial temperature as S1 and same magnetization as S2, reveals, however, that although a powerlaw index similar to that in S2 is obtained, shock-drift and diffusive-type processes now contribute in about equal amounts to particle energization; accordingly, the angular map is less anisotropic than in S2, more than in S1. Overall, this suggests that both plasma temperature and magnetization may influence the prevalence of one mechanism over the other. While a comprehensive explanation for this change of regime certainly deserves further investigation, we emphasize that the main result of the present work, namely, the formation of a nonthermal spectrum at relativistic, magnetized turbulent shocks, is a robust feature.
Let us finally address simulation S3, characterized by a relativistically hot initial plasma and a substantial magnetization σ δB ∼ 0.1 (Appendix A). This simulation probes a new regime in which stochastic Fermi acceleration inside the pre-shock turbulence controls the acceleration process, because the (stochastic) acceleration timescale t acc ∼ γ ∞ c /σ δB c now becomes short enough ( 10 3 ω −1 p ) to energize the freshly injected particles before they attain the shock. Accordingly, the spectrum shown in Fig. 2 (lower panel) does not vary with time because the turbulence inside the simulation box is stationary, up to its fluctuations. This figure also reveals that the particle distribution has undergone significant heating beyond the simple shock-crossing conditions, as can be seen by direct comparison to simulation S2 (middle panel). Furthermore, we have verified that the same simulation, albeit with open boundary conditions to prevent shock formation, yields a similar spectrum. Finally, the spectral index s ∼ 3.5 falls in line with that observed in PIC simulations of turbulence in the semi-relativistic regime [7,8,10].

IV. DISCUSSION
Interestingly, the range of magnetizations that we consider here, σ ∼ 0.01 → 0.1, and the range of spectral indices that we measure, s ∼ 2.5 → 3.5, appear rather typical of what is inferred from one-zone models of blazars [45] and gamma-ray bursts [46]. This supports the idea that mildly relativistic shocks interacting with magnetized turbulence can play a leading role in dissipation and particle acceleration in a broad range of relativistic high-energy sources, up to moderate magnetizations. Our study thus significantly extends the realm where relativistic shock acceleration can operate without turbulence (i.e. σ 10 −4 ). As stochastic turbulent acceleration is observed to take over shock acceleration at σ 0.1, one is tempted to sketch a picture in which, as the magnetization level rises, a shock, or a shock plus turbulence, then turbulence and eventually magnetic reconnection control dissipation and acceleration.
As noted, the present simulations are computationally expensive, which limits a broad parameter study. Future works should explore a larger parameter range, in particular larger dimensions (and dimensionalities) in order to examine how the particle spectrum changes with increasing c , to make better contact with phenomenology. This work was granted access to the HPC resources of TGCC under the allocations 2019-A0050407666, 2020-A0080411422, 2021-A0080411422 and 2022-A0130512993 made by GENCI. We thank X. Davoine for his assistance on particle-tracking diagnostics.

Appendix A: Simulation parameters and turbulence generation scheme
Our PIC simulations are conducted with the fully electromagnetic and relativistic calder code [33] which has been optimized to expunge beam-grid numerical instabilities, known to affect relativistic shock simulations [34,35]. We adopt a 2D3V geometry, with periodic boundary conditions in the transverse direction. In the longitudinal direction, particles are continually injected with mean velocity v ∞ = −0.87cx through the righthand side boundary. At the left-hand boundary, conditions are either open or reflective for fields and particles, as discussed below. The mesh size is ∆x = ∆y = 0.1 c/ω p , where ω p ≡ (4πn ∞ e 2 /m e ) 1/2 represents the nonrelativistic plasma frequency of the far-upstream (injected) pair plasma, with n ∞ = 2n 0 /γ ∞ as the total proper density, n 0 the apparent density of one species and γ ∞ = (1 − v 2 ∞ /c 2 ) −1/2 . The simulation domain has dimensions of 48 000 ∆x × 6 000 ∆y. The time step is ∆t = 0.99 ∆x/c. A uniform magnetic guide field B 0 is applied along the (out-of-plane) z direction. Each species (electrons or positrons) of the drifting plasma is initially represented by 10 particles per cell.
Immediately after injection, the drifting plasma is subject to turbulence forcing in its proper frame via a Langevin antenna scheme [36]. In detail, an external random current j z,ext = (c/4π)∇ 2 A z with A z = Nw i=1 a i (t )e ik i ·r excites external magnetic perturbations δB x and δB y . The coefficients a i (t ) obey the equation of a stochastically driven, damped harmonic oscillator. We use N w = 24 plane waves, with mean wavenumber k 2.9×2π/L y (primed quantities are measured in the comoving plasma frame and L y = 600 c/ω p denotes the transverse box size). Numerically, the excitation scheme is implemented so as to have it evolved in the simulation grid, while the antenna external vector potential is evaluated on refined comoving grids with x max = 4 800 c/ω p the right-hand boundary of the domain, where particles are injected. While the choice of δB · B 0 = 0 and k · B 0 = 0 points to the excitation of Alfvén modes, we stress that no velocity perturbations are excited externally. Furthermore, the fluctuations are initialized with a large amplitude, which places them in the nonlinear regime. The role of this stirring is thus to initialize the system off-pressure balance so that it evolves rapidly towards a turbulent state. Consequently, we expect the resulting turbulence cascade to comprise a significant fraction of compressive modes, as discussed in Ref. [7]. By construction, our reduced dimensionality excites an anisotropic turbulence with k k, where k = k z denotes the wavenumber component parallel to the mean field. Clearly, however, 3D simulations of the shock-turbulence interaction problem that we study remain prohibitive at the present time. Nevertheless, 3D and 2D simulations with an out-of-plane magnetic field have been shown to share the same characteristic in terms of turbulent cascade and particle acceleration, provided the turbulence level is large enough, as is the case here [10,47].
Turbulence is excited over a few coherence lengths c in the vicinity of the right-hand boundary, x max − 600c/ω p ≤ x ≤ x max . Initially, the left-hand boundary (x = 0) is left open to let the plasma exit freely. For all simulations but S2b, this boundary condition is turned to reflective once the turbulent part of the plasma has crossed the domain. This triggers a turbulent shock that mimics the interaction of two similar turbulent flows. As such, once the shock forms in the box, it propagates at a roughly constant speed. The numerical parameters characterizing our reference simulations S1, S2 and S3 and the ancillary ones S2a (no turbulence), S2b (no shock) and S4 are compiled in Table I. For S2b, the boundary conditions are left open  on the left-hand side at all times, so that the turbulent plasma can exit freely the simulation domain, which prevents shock formation.
We have also carried out additional runs in which the shock is triggered right at the onset of the simulation, in order to investigate the interaction of a shock with an ambient plasma that progressively becomes turbulent as time passes. While we do not observe a significant difference in terms of particle acceleration, the shock evolves in time, as it first interacts with a nonturbulent plasma, then with a turbulent plasma, which itself evolves downstream as it progressively fills the whole post-shock region.
Appendix B: Particle-based diagnostics Regarding particle tracking, we follow a large number (∼ 2×10 5 ) of particles, injected at different times and locations upstream of the shock, in order to study different histories of interaction with the upstream turbulence and the shock. We note that only a small fraction ( 1 %) of the tracked particles return to the shock after bouncing specularly on the reflective wall; such orbits do not therefore alter our results.
We extract the mean shock velocity from a x − t diagram of the plasma density in order to assign an analytical form to the shock front trajectory, which is thus averaged over the transverse dimension of the simulation box. The amount of time (∆t sh ) spent by each particle in the vicinity of the shock is obtained by comparing the trajectory of the particle with that of the shock front and locating the first and last shock-crossing times. In Fig. 4 of the main text, we present the correlation obtained between the energy gained (∆γ) by the tracked particles during ∆t sh with the theoretical estimate for shock-drift acceleration, namely, W (E 0 ) ≡ qE 0 ∆y, in terms of the mean electric field E 0 = v ∞ B 0ŷ and ∆y the displacement along y. We also provide a correlation with the work performed by the motional electric fields carried by the incoming turbulent fluctuations drifting at the mean bulk velocity, δE −v ∞ δB yẑ . These electric fluctuations are mostly directed alongẑ, because in our 2D3V simulations, the magnetic fluctuations mostly lie in the simulation plane while the plasma drifts alongx. W (E 0 ) accounts for shock-drift acceleration, while W (δE z ) characterizes diffusive-type Fermi acceleration. The latter mechanism must be distinguished, of course, from purely stochastic acceleration, for which the contributing electric field scales as −δv×δB, with δv the turbulent velocity component as measured in the upstream rest frame; those provide smaller contributions, as we have checked.  5. Normalized angular distributions of (left) positrons and (right) electrons, as obtained in simulation S1 at t = 8 700ω −1 p , for particles with Lorentz factor γ ≥ 20; see text for further details.
The amount of work performed by the various electric fields is measured over the time interval ∆t sh between the first and last shock crossings for each particle. The correlation is then built from the sample obtained. In simulation S1, the particles that we track have Lorentz factors γ ≥ 20, while in S2, γ ≥ 200.
Finally, as discussed in the main text, we have also recorded angular distributions of the accelerated particles to probe a possible anisotropy associated with the drift along the mean motional electric field E 0 ŷ. The angular distribution has been extracted at the final time of the simulation from the whole sample of particles (not only the tracked ones) that lie within ±250 c/ω p of the shock front, and whose Lorentz factor γ > 20 for simulation S1, and γ > 200 for simulation S2, i.e. in the powerlaw tail. Those angular maps are displayed in Fig. 5 for simulation S1 and Fig. 6 for S2, with for each a map for positrons (left) and one for electrons (right). The longitude is here defined as φ ≡ arctan2 (p y , p x ) ∈ [−π, π] and the latitude as θ ≡ π/2 − arccos (p z /p) ∈ [−π/2, π/2], in terms of the momentum components (p x , p y , p z ) and norm p. As discussed in the main text, the positron and electron angular distributions are strongly anisotropic in S2, with bright peaks at approximately opposite angles relative to the shock normal. In detail, positrons (electrons) are preferentially drifting with p y < 0 (resp. p y > 0), consistent with shock-drift acceleration (noting that E 0 < 0 here, since B 0 > 0 and v ∞ < 0). In S1, the suprathermal particles show a much broader angular distribution with weaker differences between electrons and positrons, consistent with the dominant mechanism of diffusive-type