Long-term Evolution of Relativistic Unmagnetized Collisionless Shocks

We study a relativistic collisionless electron–positron shock propagating into an unmagnetized ambient medium using 2D particle-in-cell simulations of unprecedented duration and size. The shock generates intermittent magnetic structures of increasingly larger size as the simulation progresses. Toward the end of our simulation, at around 26,000 plasma times, the magnetic coherence scale approaches λ ∼ 100 plasma skin depths, both ahead and behind the shock front. We anticipate a continued growth of λ beyond the time span of our simulation, as long as the shock accelerates particles to increasingly higher energies. The post-shock field is concentrated in localized patches, which maintain a local magnetic energy fraction ε B ∼ 0.1. Particles randomly sampling the downstream fields spend most of their time in low field regions (ε B ≪ 0.1) but emit a large fraction of the synchrotron power in the localized patches with strong fields (ε B ∼ 0.1). Our results have important implications for models of gamma-ray burst afterglows.

1. INTRODUCTION Gamma-ray bursts (GRBs) are powerful cosmic explosions driven by the collapse of a massive star or the merger of binary neutron stars (Mészáros 2002;Piran 2004;Kumar & Zhang 2015).The prompt emission stemming from the internal dissipation of the GRB jet is followed by an afterglow powered by the external relativistic collisionless shock, which expands into the surrounding medium.The afterglow unfolds from a few seconds up to decades after the explosion, as we are now witnessing for GRB 221009A, the "brightest of all time" (Burns et al. 2023;Williams et al. 2023;Laskar et al. 2023).The emission from GRB afterglows is distinctly nonthermal and broadband.It ranges from the radio to the gamma-ray band (Abdo et al. 2009;Ackermann et al. 2010;De Pasquale et al. 2010;Chandra & Frail 2012;Ackermann et al. 2014;Perley et al. 2014), and has been recently observed even at TeV energies (Abdalla et al. 2019;H. E. S. S. Collaboration et al. 2021;MAGIC Collaboration et al. 2019a,b;Huang et al. 2022a;LHAASO Collaboration et al. 2023).
GRB afterglow shocks propagate into a very weakly magnetized ambient medium.The ratio of the magnetic to the particle energy density of the cold ambient medium (i.e., the Corresponding author: Daniel Grošelj daniel.groselj@kuleuven.bemagnetization) is σ ≲ 10 −5 for massive stellar progenitors (Crowther 2007), and could be as low as σ ∼ 10 −9 if the shock expands into an interstellar-like medium (e.g., Sironi et al. 2013;Grošelj et al. 2022).A mechanism for magnetic field generation must be present at the shock in order to facilitate efficient particle scattering and their acceleration via the Fermi process (e.g., Blandford & Ostriker 1978;Bell 1978;Drury 1983;Blandford & Eichler 1987;Achterberg et al. 2001), and to explain the resulting broadband emission.Kinetic particle-in-cell (PIC) simulations of relativistic collisionless shocks (Spitkovsky 2008a,b;Martins et al. 2009;Nishikawa et al. 2009;Haugbølle 2011;Sironi et al. 2013) have shown that the Weibel filamentation instability (Fried 1959;Weibel 1959;Medvedev & Loeb 1999;Silva et al. 2003;Achterberg et al. 2007;Bret et al. 2014;Takamoto et al. 2018;Lemoine et al. 2019a) generates the required fields, which serve as the scattering agents for the Fermi process.However, Weibel-generated fields grow on plasma microscales and decay quickly downstream from the shock (Gruzinov 2001;Chang et al. 2008;Lemoine 2015), which is at odds with the relatively high field strengths inferred for some GRB afterglows (e.g., Wijers & Galama 1999;Panaitescu & Kumar 2002), including the recent GRB 221009A (Laskar et al. 2023).Moreover, particle scattering in microscale fields presents a challenge for the acceleration to very high energies (Kirk & Reville 2010;Sironi et al. 2013;Reville & Bell 2014;Asano et al. 2020;Huang et al. 2022b).The largest PIC simulations reported in the literature (Keshet et al. 2009) found that the acceleration of particles to higher energies drives the formation of magnetic fields on progressively larger scales over time, which slows down the decay of the post-shock magnetic fields.However, despite of the computational progress, previous as well as present PIC simulations need to be extrapolated to orders of magnitude longer time scales to make direct contact with observations, and so the long-term evolution of the shock remains a subject worth investigating.
In this Letter, we revisit the microphysics of relativistic unmagnetized collisionless shocks using PIC simulations of unprecedented duration and size.We observe the generation of intense magnetic structures in the upstream flow, reaching sizes up to ∼ 100 plasma skin depths.The bulk of the postshock plasma becomes magnetized, which leads to particle trapping in the magnetostatic turbulence, and prevents the fast decay of magnetic fields predicted for untrapped particles (Chang et al. 2008;Lemoine 2015).The downstream field is concentrated in localized patches, which occupy roughly one percent of the total volume but contain about half of the magnetic energy.Our results have direct implications for particle acceleration and emission from GRB afterglow shocks.

METHOD
We perform 2D simulations of relativistic collisionless shocks using the PIC code OSIRIS 4.0 (Fonseca et al. 2002;Fonseca et al. 2013).In order to evolve the simulations for as long as possible we employ for simplicity an electronpositron pair composition.As is routinely done (e.g., Sironi et al. 2013), the shock formation is driven by the reflection of the cold upstream flow (moving to the left) from the left simulation boundary at x = 0.The rest frame of the simulation coincides with the downstream frame and the shock propagates in the positive x direction.In order to reduce re-flecting boundary artifacts, the initial upstream flow Lorentz factor transitions linearly from unity to the far-upstream value γ 0 = 10 over a distance of 400 c/ω p from x = 0, where c/ω p = (γ 0 m e c 2 /4πn 0 e 2 ) is the plasma skin depth and n 0 is the simulation-frame density of the upstream electrons and positrons.At t = 0, the upstream plasma fills the transition layer between x = 0 and x = 400 c/ω p .As the simulation progresses, the region x > 400 c/ω p is gradually populated with upstream plasma introduced from a moving particle injector, which moves to the right at light speed c until it reaches the end of the box.
Our fiducial simulation has a box size of 20000×1248 c/ω p and is evolved up to tω p = 26100.The cell size is ∆x = 0.2 c/ω p and our time step ∆t = 0.5∆x/c.The upstream plasma is represented with 144 particles per cell per species with cubic spline shapes.The total number of particles at the end of the simulation exceeds 0.4 trillion.Particle noise is further reduced by low-pass filtering the electric currents at each step.The fields are advanced using a modified stencil (Blinne et al. 2018) with improved numerical dispersion of electromagnetic waves to suppress the numerical Cherenkov instability (see Grošelj et al. 2022).For improved computational performance, we employ a vectorized particle push and deposit, together with a hybrid MPI/OpenMP parallelization strategy, all of which are readily available in the OSIRIS code (Fonseca et al. 2013).
3. SHOCK EVOLUTION Figure 1 shows the evolution of the shock structure in our fiducial simulation, as observed through the snapshots of the out-of-plane magnetic field (top panels) and the electron density (bottom).The particles returning from the shock into the upstream flow drive plasma streaming instabilities, which grow and amplify magnetic fields.The upstream region populated with the self-generated fields defines a turbulent precur-sor, which expands with time as the returning particles propagate further into the upstream.At early times (tω p ≲ 4000; left panels in Fig. 1), the upstream plasma near the shock front develops a microscale filamentary structure, typical of shocks mediated by the Weibel instability (e.g., Lemoine et al. 2019a).In contrast, at later stages of the shock evolution (right panels in Fig. 1), the upstream Weibel filaments undergo a secondary nonlinear instability (Honda et al. 2000;Peterson et al. 2021Peterson et al. , 2022;;Grošelj et al. 2022), which results in large-scale magnetized plasma cavities. 1he structure of the precursor in the plasma cavitydominated phase is illustrated in Fig. 2. Cavity formation is triggered by the development of a net electric current and charge in the returning population of electrons and positrons (defined as having γβ x > 0), which stream against the incoming flow (Fig. 2(e)).For example, at the time of Fig. 2, the returning particles carry an excess positive charge.The asymmetry of the returning particles is then imprinted onto the current filaments of positive (J x > 0) or negative (J x < 0) polarity in the incoming background plasma (Fig. 2(d)), which tries to neutralize the returning-particle charge and current.The current filaments with the same polarity as the returning-particle current (J x > 0 in Fig. 2(d)) become unstable to the cavitation instability and expand in the transverse (y) direction; the filaments of the opposite polarity (J x < 0 in Fig. 2(d)) do not expand.In the rest frame of the simulation, the electric current inside the cavities is sustained by a low-density population of incoming particles of the opposite charge sign to the returning-particle charge density (Figs.2(c)-2(e)).The current inside the cavities is left essentially unscreened, which enables the generation of intense magnetic fields (Grošelj et al. 2022).The newly born plasma cavities then grow and merge as they propagate toward the shock (Fig. 2(b)), which allows the magnetic field to cascade to progressively larger scales (Fig. 2(a)).At the end of our simulation, the largest individual cavities reach transverse sizes of about 200 plasma skin depths (Fig. 1, right panels).
The current and charge of the returning particles, which drives cavity formation, exhibits oscillations in sign over time and in space.This gives rise to bursts (or cycles) of cavity generation with alternating polarity, as demonstrated in Fig. 3.The returning-particle charge fluctuations originate near the shock with a time-varying sign and spread into the upstream (Fig. 3(g)).The amplitude of the charge fluctuations amounts to a significant fraction of en r , where n r is the returning-particle number density (Fig. 3(h)).In response to the charge oscillations in the stream of returning particles, the polarity of the magnetized plasma cavities changes over time (Figs. 3(a)-3(e)).To quantify the symmetry breaking in the induced magnetic field structure we introduce an "order parameter" Q(x), defined as Q(x) = ⟨sgn(−A x )⟩ y , where A x is the longitudinal component of the magnetic vector potential (Fig. 3(f)).We solve for A x in the Coulomb gauge, so that −A x can be related via Ampere's law to the current as ∇ 2 (−A x ) = (4π/c)J x .Essentially, Q(x) measures the relative difference in the volume occupied by magnetic structures with A x < 0 versus those with A x > 0, and as such it provides a convenient diagnostic for tracking the changes in polarity over time and in space.In our fiducial run, the first cavity cycle develops a positive polarity (see Fig. 3(b) and the corresponding profile of Q(x)).It starts around tω p ≈ 4500 and lasts until tω p ≈ 11000 or so.The second cycle displays a negative polarity (Fig. 3(d)).It starts around tω p ≈ 13500, continues up to tω p ≈ 23000, and transitions into the third cycle, which features cavities of positive polarity (Fig. 3(e)).The third cavity cycle is still ongoing at the end of the run (Fig. 1, right panels).
The origin of the returning-particle charge and current oscillations is connected to the mechanism of particle reflection at Weibel mediated shocks, which was recently studied in Parsons et al. (2023).The authors found that returning particles of a given species (electrons or positrons) need to "find" current channels with a matching sign (J x < 0 for electrons and J x > 0 for positrons when the shock travels in the pos- In (h) we show the transversely averaged number density of the returning particles.Squares in (f) show the location of cavity birth (see Appendix A). itive x direction, as in our run) in order to efficiently move away from the shock into the upstream.This is supported by Figs.2(d)-2(e), which show that the local sign of the returning-particle charge density ρ r correlates with the local sign of the current J x , within a few hundred c/ω p from the shock.The injection of particles back into the upstream through the current channels provides a mechanism for spontaneous symmetry breaking.Namely, if the distribution of the background plasma current develops a seed asymmetry, one of the returning species finds current channels with a matching sign more easily and is more successful at moving away from the shock, and so the returning population as a whole develops a net current and charge (Fig. 2(e)).The returning particles of the species that sustains the excess charge then propagate further into the upstream and drive the formation of cavities, which amplifies the original seed asymmetry of the background plasma current, reinforcing the positive feedback loop.However, the continuous promotion of returning particles with a given charge sign generates a growing electrostatic potential between the shock and the far upstream.The electrostatic field pulls the species that provides the excess charge toward the shock and feeds back negatively on the asymmetry.Instead of continuous generation of cavities with a fixed polarity we therefore expect oscillations with a time-varying polarity, which are indeed observed in our simulations.Note that when the composition includes ions, the difference in inertia between electrons and ions acts as a natural mechanism of symmetry breaking (Lemoine & Pelletier 2011;Grošelj et al. 2022).Thus, the exact nature of cavity cycles in electron-ion shocks could differ from that in pair shocks.2 4. MAGNETIC FIELD STATISTICS In Fig. 4 we characterize the evolution of the magnetic field behind (left column) and ahead of (right column) the shock.Owing to cavity cycles (see Sec. 3), the evolution of magnetic field properties is not strictly monotonic.On the other hand, well-defined trends in the evolution can still be identified when the full time span of the simulation is considered.In particular, the magnetic field transverse coherence scale3 λ in front of the shock (Fig. 4(d)) grows by an order of magnitude over the duration of the simulation and approaches λ ∼ 100 c/ω p at late time.The corresponding late-time magnetic spectrum (Fig. 4(f)) is broadband; the magnetic energy uniformly fills the transverse wavenumber range 0.02 ≲ k y c/ω p ≲ 0.5.The size of structures in front of the shock correlates with the distance ∆x cav between the shock and the upstream location where the cavities arriving at the shock are born.In Fig. 4(d) we compare measurements of ∆x cav (gray dots) with the evolution of λ and find empirically an approximately linear scaling λ ≈ 0.012∆x cav . 4We discuss the possible implications of this result in Sec. 6.
Magnetic structures generated over the precursor are compressed at the shock and transmitted into the downstream.Structures of growing size ahead of the shock thus imply that more energy is transferred to longer wavelength modes in the downstream, which results in a slower overall decay of the magnetic energy behind the shock (Figs.4(a is the magnetic energy fraction.The mean of ε B , measured in the same downstream slab as the magnetic spectrum, is ⟨ε B ⟩ ≳ 10 −3 (Fig. 4(a)), except at early times (tω p ≲ 4000).For ⟨ε B ⟩ ∼ 10 −3 particles with γ ∼ γ 0 become marginally magnetized at a wavenumber k y ∼ π/R L (γ 0 ) ∼ 0.1 ω p /c, which is near the peak of the magnetic spectrum (Fig. 4(e)).This confirms that the bulk of the post-shock plasma becomes magnetized in the long-time regime.A similar finding was reported by Keshet et al. (2009) in a simulation evolved up to tω p ≈ 12600.
The decaying magnetic field behind the shock evolves into a patchy configuration composed of long-lived structures with strong fields (reaching values up to ε B ∼ 0.1) and a background with weaker fields (ε B ≪ 0.1).This is reflected in the measured probability distribution of ε B behind the shock (Fig. 4(g)).The distribution (with logarithmic bins for ε B ) rises as ε B f (ε B ) ∝ ε 1/2 B on the low energy side toward the peak around ε B ∼ 10 −3 .The 1/2 slope is consistent with a normal Gaussian distribution for the fluctuating field B z . 5Due to intermittent magnetic structures, the distribution does not cut off exponentially above the peak but instead transitions into a hard power-law tail with a slope around −1/2 ).The tail extends up to ε B ∼ 0.1 at late times. .Fraction of the volume (blue) and of the magnetic energy (orange) contained in downstream (x < x sh ) regions with εB larger than a given threshold εB.The value εB * = ⟨ε 2 B ⟩/⟨εB⟩ (dotted vertical line) is much greater than the mean ⟨εB⟩ (dashed line); εB * represents the typical local energy density of intense magnetic fluctuations, which make a significant contribution to the synchrotron emission (see text for details).The inset shows the magnetic energy content (arbitrary units) per logarithmic interval in εB.
The relatively hard power-law tail above the peak of the probability distribution of ε B in the downstream has an immediate implication for models of GRB afterglows.Namely, a significant fraction of the synchrotron power P syn ∝ γ 2 ε B from particles randomly sampling the downstream field will be emitted at the localized structures with ε B ∼ 0.1.This is demonstrated in Fig. 5, which shows the fraction of magnetic energy and the volume fraction contained in downstream regions with ε B larger than a given threshold εB .Roughly ∼ 50% of the total magnetic energy, and hence of the total synchrotron power, is contributed by structures with ε B ≳ 0.1, which occupy about ∼ 1% of the downstream volume.That a significant fraction of the emission is contributed by the high field regions is further supported by the downstream magnetic energy distribution (inset of Fig. 5, showing ε Our results motivate the design of two-zone synchrotron emission models (e.g., Kumar et al. 2012; Khangulyan et al.   5 2021), where particles spend most of their time in low field regions, but emit an order unity fraction of their synchrotron power intermittently at the locations with strong fields. 6In the scenario implied by our simulations, the average synchrotron cooling rate (for a given γ) is determined by the mean value ⟨ε B ⟩ ≪ 0.1 of the magnetic energy fraction in the emission region (see Fig. 5).In contrast, the photon energy E syn ∼ γ 2 ℏeB * /m e c at the peak of the emitted spectral energy distribution is set by the high field strength B * , corresponding to the peak of ε 2 B f (ε B ) at ε B * ∼ 0.1 (inset of Fig. 5), since it is the value at the peak of ε 2 B f (ε B ) that makes the single largest contribution to the emission.A convenient analytic estimate of ε B * can be given based on the synchrotron-power-weighted mean of ε B , which equals ⟨ε 2 B ⟩/⟨ε B ⟩ (for fixed γ).In our simulation, indeed ε B * ∼ ⟨ε 2 B ⟩/⟨ε B ⟩ ∼ 0.1.Using our analytic estimate, we can express E syn ∼ κ 1/2 γ 2 ℏeB rms /m e c, where κ = ⟨ε 2 B ⟩/⟨ε B ⟩ 2 is the kurtosis of the magnetic fluctuations and B rms ∝ ⟨ε B ⟩ 1/2 is the root-mean-square field strength.The downstream field at the end of our simulation has κ ≈ 25 on average, which is much greater than the value κ = 3 for normal Gaussian distributions.As a result, E syn is raised by a factor of ∼ 3 compared to the emission in a non-intermittent field with κ = 3.The factor could be increased further for more intermittent distributions with even larger κ.Note that the factor ∼ κ 1/2 applies to all the emitting particles, including the highest energy ones.If the highest-energy particles are accelerated at maximum efficiency (i.e., close to the Bohm regime with λ ∼ R L (γ)), the intermittency of the post-shock magnetic field enables the production of synchrotron photons beyond the classical burn-off limit (de Jager et al. 1996;Kumar et al. 2012;Khangulyan et al. 2021).
5. PARTICLE ACCELERATION Finally, we study particle acceleration at the relativistic shock.Fig. 6 shows the evolution of the particle energy spectrum (γ 2 dN/dγ) in a fixed slab behind the shock at x−x sh = [−1200, −800] c/ω p .Particles are accelerated to higher energies as the simulation progresses, consistent with earlier works (e.g., Spitkovsky 2008a).At select times, which correlate with bursts of plasma cavity generation (around tω p ∼ 10 4 and toward the end of the simulation; see Fig. 3), the particle maximum Lorentz factor in the slab grows faster than the canonical γ max (t) ∝ (tω p ) 1/2 scaling (Sironi et al. 2013;Plotnikov et al. 2018), expected for Weibel-mediated relativistic shocks (see inset of Fig. 6).However, in between the intervals of fast acceleration lies a quiescent phase with a very slow growth of γ max (t), so that the maximum Lorentz factor still traces out a ∼ (tω p ) 1/2 envelope on average.The rapid changes in the instantaneous rate of particle acceleration can be largely attributed to the evolution of the magnetic field properties, which set the particle scattering frequency as ν scatt ∼ λc/R 2 L ∼ ε B (λω p /c)(γ/γ 0 ) −2 ω p (e.g., Plotnikov et al. 2011).An integration of dγ/dt = 0.25 ν scatt γ (using 0.25 as an ad hoc prefactor) with a time-dependent ⟨ε B ⟩⟨λω p /c⟩, measured in the same downstream slab as the particle spectrum (see Fig. 4), gives a result broadly consistent with the measured γ max (t) (gray dotted curve in the inset of Fig. 6).A reasonable overall agreement is also obtained for a fixed ⟨ε B ⟩⟨λω p /c⟩ ≈ 0.4 (dashed black curve).A novel feature observed in the particle spectrum at late times is the formation of a distinct suprathermal component, which connects the Maxwellian peak (near γ ∼ γ 0 = 10) to the high-energy power law tail with an index close to −2 at late time (i.e., dN/dγ ∝ γ −2 at 100 ≲ γ ≲ 200).The fraction of kinetic energy contained in the suprathermal component grows with time.Toward the end of the run, the suprathermal component contains about 30% of the total kinetic energy, while the high-energy tail contains stably around 10%.The suprathermal component grows at the expense of the thermal (i.e., Maxwellian) part of the distribution, which loses energy over time.It is interesting to note that the start of the high-energy tail shifts with time to higher γ, in favor of the growing suprathermal component.In this regard, let us mention that standard treatments of the first order Fermi process at relativistic shocks (e.g., Achterberg et al. 2001) assume a discontinuous flow profile across the shock transition.However, in practice the shock develops a broad transition region over which the incoming flow gradually slows down (Lemoine et al. 2019a).The enhanced scattering due to large-scale structures and the strongly turbulent shock transition layer observed in our simulation imply that moderate energy particles (γ ∼ a few γ 0 ) returning into the upstream may deflect toward the downstream before being able to probe the full velocity difference between the far upstream and downstream.We speculate that these particles represent the dominant contribution to the suprathermal component, but we leave a more detailed investigation of this aspect for the future.
6. DISCUSSION AND CONCLUSIONS We studied the evolution of a relativistic pair plasma shock propagating into an unmagnetized medium using 2D PIC simulations of unprecedented duration and size.Toward the end of our fiducial run, at tω p ≈ 26100, the shock generates magnetic structures reaching sizes of roughly ∼ 100 plasma skin depths on both sides of the shock.The size of magnetic structures generated in our simulation may be sufficient to circumvent at least some of the outstanding issues related to the modeling of GRB afterglows.For instance, it has been argued (Huang et al. 2022b) that magnetic fields with a coherence scale of about ∼ 100 ion skin depths or larger are required to explain the X-ray and gamma-ray afterglow of GRB 190829A (H.E. S. S. Collaboration et al. 2021). 7It should be noted that further evolution of the shock structure is expected beyond the time span of our simulation, as discussed below.We also find that the bulk of the post-shock plasma becomes magnetized, which slows down the magnetic field decay, as compared to early times of the shock evolution (Sec.4).The probability distribution of the downstream magnetic energy fraction ε B is intermittent; it features a hard power-law tail which extends up to ε B ∼ 0.1.Particles randomly sampling the downstream field spend most of their time in low field regions (ε B ≪ 0.1), but emit roughly half of their synchrotron power at localized magnetic structures with ε B ∼ 0.1.
The shock exhibits rich spatiotemporal dynamics, characterized by bursts of magnetized plasma cavity generation and oscillations in the sign of their polarity (Sec.3).The cavitation instability enables the formation of large-scale magnetic structures over the turbulent shock precursor, giving rise to a broadband magnetic wavenumber spectrum at the end of the simulation (Fig. 4(f)).The spectrum extends from the large energy-containing scale, characterized by the transverse coherence scale λ, down to plasma microscales ∼ c/ω p .Empirically, we find that λ measured immediately ahead of the shock is proportional to ∆x cav , where ∆x cav is the distance between the shock and the upstream location of cavity birth (Fig. 4(d)).Beyond the duration of our simulation, an upper limit on λ(t) may be obtained by assuming that the distance ∆x cav scales linearly with the size of the precursor ℓ prec , which is set by the scattering length of the highest energy particles.The upstream scattering length ℓ scatt ∝ γ 2 (Lemoine et al. 2019b).This implies λ ∝ t, assuming the shock keeps 7 For application to GRB shocks propagating into an electron-proton ambient medium, time and length scales should be measured in units of the ion (i.e., proton) inverse plasma frequency ω −1 pi = (m i /me) 1/2 ω −1 p and skin depth c/ω pi = (m i /me) 1/2 c/ωp, where m i /me is the ion-electron mass ratio.In our simulation, m i /me = 1 for computational convenience.
accelerating particles to higher energy as γ max (t) ∝ t 1/2 (see Sec. 5).A linear growth rate for λ presents a potentially favorable scenario for particle acceleration; it suggests that the ratio λ/R L (γ max ) ∼ ε 1/2 B (λω p /c)(γ max /γ 0 ) −1 grows over time and converges toward the Bohm limit (λ ∼ R L (γ max )).While these results are encouraging, we caution the reader that extrapolations beyond the time span of our simulation are uncertain, and the scaling λ ∝ t is only an optimistic estimate.
The size of magnetic structures increases from early to late time also behind the shock front (Fig. 1).However, quantitatively the increase of the field coherence scale λ in the downstream is not large (see also Keshet et al. 2009).At different times and/or at different locations behind the shock, we measure typical values of the order of λ ∼ 100 c/ω p (Fig. 4(c)).This scale corresponds within factors of a few to the wavenumber k y ∼ 0.1 ω p /c at which the magnetic spectrum peaks (Fig. 4(e)), which also happens to be the scale at which the bulk of the post-shock plasma becomes marginally magnetized.Presumably, the demonstration of a large increase of λ in the downstream requires longer duration simulations with even wider boxes.In particular, we expect the coherence scale at a fixed distance behind the shock to grow when the size of magnetic structures arriving at the shock exceeds our "typical" downstream value of ∼ 100 c/ω p .Then, as the growing magnetic structures are transmitted into the downstream, the coherence scale will increase on both sides of the shock.
Even though the presented PIC simulation is the largest and longest of its kind, extrapolations are still required in order to make direct contact with GRB afterglows.The apparent duration of our simulation in the observer frame (Mészáros 2006) is T obs ≃ T (1 + z)/2γ 0 ≈ 0.02 (γ 0 /10) −1 (1 + z)(n u /1 cm −3 ) −1/2 (m i /m e ) 1/2 s, where T is the duration in the (downstream) simulation frame, n u is the upstreamframe number density of the ambient medium, z is the redshift, and m i /m e is the ion-electron mass ratio (m i /m e = 1 in our run, but for electron-proton shocks m i /m e ≈ 1836).The transverse size of the box is L ⊥ ≈ 7 × 10 8 (n u /1 cm −3 ) −1/2 (m i /m e ) 1/2 cm.Thus, our present simulation probes the shock physics over time and length scales intermediate between the plasma microscales and the global macroscales of the blast wave.In our fiducial run, spatial regions featuring cavities of a given polarity extend up to thousands of skin depths in the longitudinal direction, and across the full transverse dimension of the simulation box.In a macroscopically large domain, these long-range correlations will be ultimately limited by causality.In that case, we expect the shock to form cavities of different polarity in causally disconnected regions transverse to the shock normal.Our present simulation box is too narrow to cancel out the system-wide correlations.Strictly speaking, the exact times and locations of cavity generation observed in our fiducial run represent only a particular realization of our numerical experiment.Nevertheless, the general trend of large-scale field generation is still clearly evident when considering the full time span of the simulation.
Previous PIC simulations found that as much as ∼ 90% of the post-shock kinetic energy is contained in thermal (i.e., Maxwellian) particles (e.g., Spitkovsky 2008a,b;Martins et al. 2009).Radiative signatures of thermal electrons from GRB afterglows have been considered by a number of authors (e.g., Eichler & Waxman 2005;Giannios & Spitkovsky 2009;Ressler & Laskar 2017;Laskar et al. 2019;Warren et al. 2022).However, conclusive observational evidence for the presence of a prominent thermal component is presently lacking.In contrast to previous PIC simulations, we find that the fraction of energy contained in the thermal component drops over time.At the end of our fiducial run, the fraction of energy in the Maxwellian component drops from the "canonical" ∼ 90% to ∼ 60% (Sec.5), which reduces the tension with the widely used assumption of purely nonthermal electrons in afterglow models (e.g., Sari et al. 1998;Granot & Sari 2002).
In order to evolve the shock simulations for a maximum possible duration we focused on a pair plasma composition.Previous works showed that the physics of electron-ion relativistic unmagnetized shocks is similar to pair shocks, owing to efficient preheating of the incoming electrons (e.g., Spitkovsky 2008b;Martins et al. 2009;Sironi et al. 2013).Therefore, we expect qualitatively similar results for electronion shocks.A critical parameter affecting the shock evolution is the upstream magnetization of the ambient medium σ (see Sironi et al. 2013;Reville & Bell 2014;Plotnikov et al. 2018).The physical picture described here requires that the shock operates close to the σ = 0 limit.Presently, it is not well understood how small σ needs to be for the σ = 0 limit to be relevant.It has been shown that different levels of ambient magnetization can dramatically alter the behavior of collisionless shocks (e.g., Bret et al. 2017;Bret & Narayan 2018;Haggerty et al. 2022;Grošelj et al. 2022).Thus, it may be reasonably expected that only a fraction of GRBs take place in environments where a vanishingly small magnetization can be assumed, which might help explaining the relatively large scatter of best-fit modeling parameters with respect to different afterglow observations.

Figure 1 .
Figure 1.Structure of a collisionless relativistic pair plasma shock, propagating into an unmagnetized upstream medium with a bulk Lorentz factor γ0 = 10, at tωp = 4000 (left) and at tωp = 26100 (right).On the top we show the out-of-plane magnetic field Bz, and on the bottom the electron number density n e − .An animated version of this figure is available online.The animation lasts 40 s and shows the spatiotemporal evolution of Bz (top panel) and n e − (bottom panel) from the start (tωp = 0) to the end (tωp = 26100) of the simulation.

Figure 2 .
Figure 2. Structure of the shock precursor populated with magnetized plasma cavities at tωp = 7000.Only a fraction of the box is shown in the range 100 < yωp/c < 400 and 0 < [x−x sh (t)]ωp/c < 1300, where x sh (t) is the location of the propagating shock front.From top to bottom we show the z component of the magnetic field Bz (a), the particle number density n (b), the total electric charge density ρ (c), the x component of the electric current Jx (d), and the charge density ρr of the returning particles with γβx > 0 (e).The inset plots zoom in on the structure of a single cavity.

Figure 3 .
Figure 3. Evolution of the magnetic field and the emergence of symmetry breaking.Panels (a)-(e) show snapshots of Bz at five times in the simulation.In panel (f) we show the "order parameter" Q(x) = ⟨sgn(−Ax)⟩y.Values of Q(x) significantly above or below zero indicate the presence of cavities of a given magnetic polarity (see text for details).Panel (g) shows the transversely averaged upstream charge density of the returning particles with γβx > 0.In (h) we show the transversely averaged number density of the returning particles.Squares in (f) show the location of cavity birth (see Appendix A).

Figure 4 .
Figure 4. Evolution of the magnetic field behind (left) and ahead of (right) the shock.Panels (a)-(d) show the average magnetic energy fraction εB = B 2 z /8πγ0n0mec 2 and the field transverse coherence scale λ as a function of time in a fixed downstream (x − x sh = [−1200, −800] c/ωp) and upstream (x − x sh = [100, 500] c/ωp) slab.Panels (e)-(f) show the magnetic spectrum as a function of the transverse wavenumber ky in the downstream (e) and upstream (f) slab (colors represent time).The spectrum is compensated by ky to highlight the energy content per scale.In (g)-(h) we show the distribution of εB in the downstream (g) and upstream (h) slab.Finally, in (i) we show the transversely averaged profile of εB.Gray dots in panel (d) measure the distance ∆xcavωp/c (rescaled by a numeric prefactor ≈ 0.012) between the shock and the upstream location of cavity birth (see Appendix A for details).
Figure5.Fraction of the volume (blue) and of the magnetic energy (orange) contained in downstream (x < x sh ) regions with εB larger than a given threshold εB.The value εB * = ⟨ε 2 B ⟩/⟨εB⟩ (dotted vertical line) is much greater than the mean ⟨εB⟩ (dashed line); εB * represents the typical local energy density of intense magnetic fluctuations, which make a significant contribution to the synchrotron emission (see text for details).The inset shows the magnetic energy content (arbitrary units) per logarithmic interval in εB.

Figure 6 .
Figure6.Evolution of the particle energy spectrum in a fixed slab behind the shock (colors represent time).The dashed red curve shows a Maxwellian fit to the low-energy part of the spectrum at the end of the simulation.The inset shows the evolution of the maximum Lorentz factor γmax(t) (brown curve), and compares the measurement with theoretical expectations for particle scattering in a magnetic field with fixed ⟨εB⟩⟨λωp/c⟩ ≈ 0.4 (black dashed curve) versus scattering in a field where ⟨εB⟩⟨λωp/c⟩ is time-evolving (gray dotted curve; see text for details).The maximum Lorentz factor is determined by the value of γ where γdN/dγ drops by 10 5 below the peak.