The maximum energy of shock-accelerated electrons in a microturbulent magnetic field

Relativistic shocks propagating into a medium with low magnetization are generated and sustained by small-scale but very strong magnetic field turbulence. This so-called"microturbulence"modifies the typical shock acceleration process, and in particular that of electrons. In this work we perform Monte Carlo (MC) simulations of electrons encountering shocks with microturbulent fields. The simulations cover a three-dimensional parameter space in shock speed, acceleration efficiency, and peak magnetic field strength. From these, a Markov Chain Monte Carlo (MCMC) method was employed to estimate the maximum electron momentum from the MC-simulated electron spectra. Having estimated this quantity at many points well-distributed over an astrophysically relevant parameter space, an MCMC method was again used to estimate the parameters of an empirical formula that computes the maximum momentum of a Fermi-accelerated electron population anywhere in this parameter space. The maximum energy is well-approximated as a broken power-law in shock speed, with the break occurring when the shock decelerates to the point where electrons can begin to escape upstream from the shock.


INTRODUCTION
When particles are accelerated by the Fermi process at a collisionless shock, the distribution takes the form of a power law in momentum, with an exponential cutoff.While energy gains from acceleration exceed losses, the distribution is that of a power law; above some momentum, however, energy losses dominate and the number of particles drops off exponentially.The resulting function may be written dN dp = Kp −σ e −(p/pmax) χ . (1) In the above equation, K is a normalization constant, σ is the spectral index of the power law, p max is the location of the exponential cutoff, and χ controls the curvature of the turnover.In the case of Bohm diffusion limited by a free escape boundary, χ = 1 (Caprioli et al. 2009).If particle energies are instead limited by radiative losses such that ṗ ∝ p 2 , then χ = 2 (Zirakashvili & Aharonian 2007;Lefa et al. 2012).In the extreme limit χ 1, the exponential turnover behaves like a step function.An additional component is expected in physical distributions, that of the particles that crossed the shock but were not further accelerated.This "thermal" population is discussed elsewhere in the literature (Giannios & Spitkovsky 2009;Warren et al. 2015;Ressler & Laskar 2017), but since it occurs at the low-energy end of the particle distribution we ignore it henceforth.For this work we are not interested in how particles are injected into the acceleration process; we focus solely on what happens to the particles once they are being Fermi-accelerated.
Both σ and p max can be predicted analytically under certain simplifying conditions.If particle orientations diffuse isotropically, and if the shock is assumed to be discontinuous, then σ depends only on the shock Lorentz factor Γ 0 ; one can then recover σ = 2 in the non-relativistic limit Γ 0 → 1 and σ ≈ 2.2 in the ultra-relativistic limit Γ 0 → ∞ (Keshet & Waxman 2005).This result has been confirmed by many numerical studies both pre-and post-dating the analytical work (Ostrowski 1988;Ellison et al. 1990;Ellison & Double 2002;Baring 2004;Aoi et al. 2008;Ellison et al. 2013).
To determine p max one must assume an energy loss mechanism.If radiative losses such as synchrotron (or synchrotron self-Compton) are responsible, then one may set the acceleration time equal to the loss time (Gallant & Achterberg 1999;Achterberg et al. 2001;Aharonian et al. 2002;Lemoine & Pelletier 2003;Piran & Nakar 2010;Bykov et al. 2012;Lemoine 2013).The resultant formula for p max has many more dependencies than that for σ; as shown in the cited works p max potentially depends on large-scale hydrodynamical parameters (e.g.shock speed) and microphysical parameters (e.g.magnetic field strength, turbulence characteristics, or scattering model).
Particle-in-cell (PIC) simulations (Nishikawa et al. 2007;Spitkovsky 2008;Sironi et al. 2013;Ardaneh et al. 2015Ardaneh et al. , 2016) ) offer a first-principles method for determining the distributions particles take as they interact with a collisionless shock.By self-consistently treating both how particles self-generate magnetic field turbulence and how they scatter in it, PIC simulations act as the "gold standard" for numerical studies.However, they are extremely computationally intensive, and this severely limits their ability to simulate the times and volumes needed to discuss the value of p max in realistic settings. 1hile PIC simulations can not yet address the value of p max , they are in general agreement about the structure of the magnetic field around the relativistic shocks they simulate.Such shocks are susceptible to a wide variety of plasma instabilities (Medvedev & Loeb 1999;Bret et al. 2010;Lemoine & Pelletier 2011), and these instabilities drive small-scale "microturbulence" in the magnetic field that can exceed 10% of the local energy density (e.g.Sironi et al. 2013).This microturbulence decays rapidly behind (and also in front of) the shock, meaning particles probe magnetic fields of different strengths and turbulence spectra.The maximum momentum attainable by particles can still be determined by equating acceleration time and loss time, but both of those quantities are more complicated than they are in a zeroth-order model (Lemoine 2013).
Monte Carlo simulations (Ostrowski 1988;Ellison et al. 1990;Ellison & Double 2002;Lemoine & Pelletier 2003;Baring 2004;Niemiec & Ostrowski 2006;Ellison et al. 2013;Warren et al. 2015) trade additional assumptions about the shock microphysics for substantially lower computational costs.Instead of handling magnetic field turbulence and particle scattering self-consistently, Monte Carlo schemes typically fix one or both.The computational cost of Monte Carlo codes is sufficiently low that they can be used to probe p max , albeit sparsely, over a meaningful range of shock parameters and magnetic field configurations.
In this work, we apply an established Monte Carlo code (Ellison 1985;Ellison et al. 1990;Ellison & Double 2002;Vladimirov et al. 2008;Ellison et al. 2013;Warren et al. 2015Warren et al. , 2017) ) to the problem of particle acceleration in the microturbulent field.We simulate electron acceleration, and radiative losses, over a multidimensional physical parameter space.From the accelerated particle spectra thus obtained, we identify an empirical formula that quickly computes p max over a wide parameter space of astrophysical relevance.

NUMERICAL METHODS
The Monte Carlo (hereafter MC) model used here assumes a steady-state shock and explicitly defines how particles scatter in the magnetic turbulence, in effect solving the Boltzmann equation with collective scattering (Ellison & Eichler 1984;Ellison et al. 2013).This steady-state approximation is valid as long as the particle acceleration time is shorter than the dynamical time of the system.
Particle motions, energy losses, and scatterings are treated in a series of discrete time steps δt.The MC model assumes that particles scatter with a mean free path proportional to their gyroradius r g = γmvc/(qB) in the local magnetic field B, where η mfp characterizes the strength of the scattering; the case of η mfp = 1 is called Bohm diffusion.Bohm diffusion is commonly assumed because of its simplicity, but there is no conclusive evidence that η mfp has any particular value, or indeed a single value at all.Observations of rapid variability in bright knots of the supernova remnant RX J1713.7-3946suggest that particle acceleration takes place with η mfp ∼ 1 (Uchiyama et al. 2007).In the ultrarelativistic limit, interpreting GeV emission in the afterglows of gamma-ray bursts (GRBs) as due to synchrotron emission places strict limits of η mfp 1 on the acceleration process (Sagi & Nakar 2012).Simulations of particle scattering in the relativistic regime variously suggest that η mfp > ∼ 10 (Lemoine & Revenu 2006) or η mfp ∼ 1 (Riordan & Pe'er 2019).Indeed, the discussion surrounding η mfp is far from resolved (and the dependence of λ mfp on r g may not be linear in some regimes (Plotnikov et al. 2011)), and for simplicity Equation 2 is assumed to hold everywhere in the shock structures simulated.
The MC model divides each gyroperiod into N g smaller propagation and scattering steps which take place over a time δt.That is, δt ∝ r g /N g .The amount by which particle orientations are allowed to change at each scattering event depends on both η mfp and N g (Ellison et al. 2013): δθ max ∝ (η mfp N g ) −1/2 .In the continuous limit of δt → 0 (equivalent to N g → ∞), we approach pitch-angle diffusion as used in (Keshet & Waxman 2005).The value of N g used in the simulations was tested to confirm that larger values did not affect the results.
The only energy gain process used in the MC model is first-order Fermi acceleration, i.e. crossing the shock and scattering off of turbulence in the new rest frame.There are numerous other mechanisms by which particles in the vicinity of a shock could gain energy, including wake field acceleration (Tajima & Dawson 1979), second-order Fermi acceleration (Virtanen & Vainio 2005, and references therein), and cross-shock potentials (e.g., Tran & Sironi 2020).We assume that once injection into the first-order Fermi process occurs, those energy gains dominate those from other sources.
The code can treat the acceleration process of particles with any charge to mass ratio, and has in the past been applied to electrons, protons, and heavier ions.For this work we restrict ourselves to considering only electrons.Electrons carry significantly less momentum and energy flux than do ions, and their ability to influence the shock structure is similarly limited.Even in relativistic shocks with Γ 0 ∼ 30, nonlinear interactions between the shock and the accelerated particles can change the velocity profile of the shock (Ellison et al. 2013).At higher particle energies, especially in the limit of relativistic shocks, the induced precursor is of sufficiently small extent that it is ignorable.Since the primary objective of this work is to estimate the value of p max for the electron spectrum given in Equation 1, nonlinear effects can be reasonably neglected.
Another break from previous work with this code is the magnetic field profile used.Nonrelativistic versions of the code self-consistently handle generation and decay of magnetic field turbulence (Vladimirov et al. 2008(Vladimirov et al. , 2009;;Bykov et al. 2014), but the relativistic version assumed that the mean field was parallel to the plane of the shock and of uniform strength throughout the shock structure (both up-and downstream from the shock).In the work presented here we assume that the magnetic field is dominated by isotropic turbulence, consistent with the downstream region in PIC simulations.The magnetic field upstream is more ordered due to the filamentation instability, but in relativistic shocks electrons spend the bulk of their time scattering in the downstream medium.The strength of the magnetic field is related to the local energy density by the parameter B .The peak value of both magnetic field and B occurs at the shock, and may be expressed in terms of bulk hydrodynamical parameters, where B is the mean local magnetic field, Γ 0 is the Lorentz factor of the shock, and n 0 is the restframe ambient density the shock is presently encountering (cf. the definition of proper energy density in Blandford & McKee 1976;Granot & Sari 2002).Upstream and downstream from the shock, the magnetic field used is consistent with PIC simulations (Keshet et al. 2009;Sironi et al. 2013) and previous analytical work (Lemoine 2013): in which the distance x from the shock is measured in ion skin depths, λ sd = (Γ 0 m p c2 /4πe 2 n 0 ) 1/2 .This functional form has a plateau within 50 plasma skin depths of the shock, then decays upstream (negative positions) and downstream (positive positions).The plateau has been used in previous analytical treatments of microturbulence in GRB afterglows (Lemoine 2013), but is a simplification of the sharper, exponential-like peak in B found in PIC simulations.The downstream decay rate of x −1 is expected from PIC simulations (Chang et al. 2008) and interpretations of GRB afterglows with microturbulence (Lemoine et al. 2013).The upstream decay rate matches the precursors presented in Sironi et al. (2013, Fig. 7); we have not seen an extended discussion of this decay rate in the literature, so we do not make any claims of universality here.
The oblique instability responsible for the microturbulence is expected to disappear when the shock speed drops below a certain critical value, Γ 0 ∼ 10 for reasonable choices of physical parameters (Lemoine & Pelletier 2011;Pelletier et al. 2017).The calculation assumes a relativistic shock, however, and it is less clear what instabilities are present in the mildly relativistic regime.At low Lorentz factors both analytical and numerical approaches encounter difficulty, and there are comparatively fewer studies in the literature.Numerical work on mildly relativistic shocks suggests that short-wavelength turbulence exists in the shock precursor at Lorentz factors as low as Γ 0 = 1.5 (Romansky et al. 2018).We therefore assume here that amplified microturbulence is present at all Lorentz factors considered, and remain agnostic as to the precise source of the turbulence.
As particles scatter in the above-mentioned magnetic field turbulence, they experience radiative losses.We consider only synchrotron losses in this work.Electrons should additionally radiate away energy via the inverse Compton process, and in the strong magnetic fields present near the shock synchrotron self-Compton (SSC) can be significant (Sari & Esin 2001).The rate of SSC losses depends on the synchrotron spectrum at an electron's location within the shock structure, which in turn depends on the local electron distribution.The full calculation of synchrotron-SSC cooling is thus highly nonlinear.Schlickeiser et al. (2010) discuss the problem of accelerating electrons that are simultaneously suffering synchrotron and SSC losses, but their work covers only interactions in the Thomson regime.As magnetic field strengths and electron energies increase, the Klein-Nishina cutoff on SSC cross-section becomes relevant and serves to suppress SSC losses.SSC in the Klein-Nishina regime has been considered by many authors, but their work assumed an injected electron distribution rather than calculating it self-consistently (Nakar et al. 2009;Wang et al. 2010;Barniol Duran et al. 2012).We are unaware of an attempt to simultaneously treat electron acceleration and SSC cooling in the Klein-Nishina regime, and so we also ignore these complications to focus instead on the much simpler synchrotron loss process.
We include another means of limiting particle energy in the code: a free escape boundary (FEB) upstream of the shock.Any particles reaching the boundary are considered to decouple from the shock, and are removed from the acceleration process. 2 Particle escape from a shock can only be identified after the fact; particles scattering ahead of a shock have a non-zero probability of returning back to the shock regardless of how far from the shock they get, and this probability is increased with an expanding shock (Drury 2011).An additional complication arises since particle diffusion lengths increase in the weaker magnetic fields far upstream from the shock.Rather than track all particles indefinitely, or introduce a random chance of escape between each scattering event, we approximate particle escape with the use of a sharp boundary.The FEB acts as an additional limitation on particle energy since diffusion lengths are proportional to energy (Equation 2).The location of the FEB for all MC simulations is where the negative sign denotes a position upstream from the shock, and That is, the physical location of the FEB depends on both the shock speed and the ambient magnetic field B 0 , which taken here to be 10 −5 G. Beyond setting the physical scale of the FEB, the ambient magnetic field is largely irrelevant because it is dominated everywhere in the MC simulations by the amplified turbulence.
It is easier for particles to scatter out to the FEB in mildly relativistic shocks for several reasons.First, Equation 5shows that the physical location is closer to the shock when Γ 0 is lower.Second, according to Equation 3 the magnetic field is weaker around slower shocks.This means longer mean free paths per Equation 2, so the same physical distance represents fewer diffusion lengths when the magnetic field is reduced.The weaker field also means that electrons experience smaller radiative losses; they therefore keep more of their energy as they scatter towards the FEB and can make the trip in fewer scattering steps.Taken as a whole, if the FEB acts to limit the maximum electron energy we expect that maximum to be an increasing function of Γ 0 β 0 ; once electrons can no longer reach the FEB, the relationship between shock speed and maximum energy will change.
All shocks simulated in this work assume that acceleration occurs in the test-particle limit; that is, the presence of accelerated particles upstream from the shock does not modify the velocity profile.In actuality, pressure from these accelerated particles acts to slow down inflowing plasma, generating a shock precursor via a nonlinear feedback loop.This nonlinear shock modification occurs at all shock speeds, from nonrelativistic to ultrarelativistic (Ellison & Eichler 1984;Reville et al. 2009;Ellison et al. 2013).However, the precursor in a relativistic shock is weak compared to those of nonrelativistic shocks (Ellison et al. 2013); so while nonlinear effects may be present, their treatment is deferred to future work.

ESTIMATING p max FROM THE MONTE CARLO RESULTS
The MC code described above was used to simulate particle acceleration at numerous points in a three-dimensional parameter space.The proper speed of the shocks, Γβ, was allowed to vary between 0.83 (corresponding to Γ 0 = 1.3) to 340.The plateau value of B in the vicinity of the shock ranged between 10 −4 and 10 −1 .The acceleration efficiency η mfp varied between 1 and 10.The strongest magnetic fields in the simulations occur at the shock.Their strength is, by the definition of B , a product of that parameter and the local energy density: This strength depends on both Γβ, by way of Γ 2 0 , and on the product B n 0 (where n 0 is the ambient density of the material being swept up by the shock).Each individual MC run can thus be characterized by a parameter triplet φ MC = (Γβ, B n 0 , η mfp ).
Over this wide, physically-relevant range of φ MC parameters, the MC code was used to simulate 20 iterations (with a different initial random number seed) for each parameter triplet.From the resulting 20 electron spectra, we estimated a posterior likelihood distribution (PLD) of the parameters π pl = (log 10 p max , σ, χ, log 10 K) via Markov chain Monte Carlo (MCMC), as detailed in Appendix A.
An example of this process is given in Figure 1, showing both the MC realizations and the MCMC-estimated spectrum.The MC electron distribution has two components: a "thermal peak" (composed of electrons that crossed the shock but did not enter the acceleration process) and the non-thermal tail (made up of electrons that did shock-accelerate).The sharpness of the thermal peak is an artifact of how the MC code handles energy gain upon shock crossing3 ; it does not impact the acceleration process that produces the non-thermal tail.The non-thermal distribution is The curves correspond to the p max values predicted by the empirical expression (Equations 8-13) in its full (15-parameter, dashed line) and reduced (10-parameter, solid line) forms; the values of the parameters are provided in Table 1 below.Gray shading marks regions of the parameter space where electrons were able to escape from the shock in the upstream direction.The rightmost column compares MC runs with the same B n 0 but different η mfp .The bottom row compares MC runs with the same η mfp but different B n 0 .
well-captured by Equation 1, including the extended power-law tail and the exponential turnover at p max .
The momentum value of the thermal peak in Figure 1 is a free parameter in this work.Its position can be fixed using either a PIC simulation or by assuming a specific fraction of energy placed in the electron distribution ( e in the literature), but we are not interested here in the low-energy portion of the spectrum.We instead place the thermal peak at a sufficiently small momentum that the power law generated by shock acceleration can cover several decades, allowing us to adequately estimate the σ parameter in Equation 1.
The values of p max for each φ MC are shown in Figure 2.While all four of Equation 1's π pl parameters are estimated for each φ MC , hereafter we focus only on p max .Each dot in Figure 2 represents the highest-likelihood value of p max from a single φ MC (cf. Figure 1).Within each subpanel, corresponding to a specific ( B n 0 , η mfp ) doublet, p max versus Γβ appears to follow a broken power law.The shaded area in each subpanel marks the region of Γβ values for which electrons were able to scatter upstream to the FEB discussed in Section 2. It appears from Figure 2 that electron escape is, at the very least, strongly correlated with the break in the p max curves.That is, below the break the maximum electron energy is limited by escape at the FEB, while above the break the limiting process is radiative losses.

AN EMPIRICAL EXPRESSION FOR p max
As seen in Figure 2, p max appears to follow a broken power law in Γβ: In this form, the parameters A, B, C, D, and E do not have well-defined physical associations.For example, while setting y = B results in p max = A, it is not the case that B corresponds to the location of the peak and A to the height at the peak.Instead, both the location and the height of the peak also depend on C, D, and E. It is therefore useful to re-arrange Equation 7 so that the five parameters describe different, mostly independent features of the broken power law.This yields whose five features are p pk , the maximum value of the broken power law; Γ pk β pk , the value of Γβ at which the function peaks; L and H, the asymptotic power-law indices far below and above (respectively) the peak; and finally W , which controls the width of the peak.
We seek an empirical formula for p max over the three-dimensional φ MC = (Γβ, B n 0 , η mfp ) parameter space.Since Equation 8 already explicitly depends on Γβ, this can be accomplished by assuming that all five features of Equation 8 depend on the remaining two parameters, B n 0 and η mfp , such that The subscript C here denotes the coefficient parameter for each feature, and the subscripts η and identify the power parameters that quantify the dependence of each feature on η mfp and B n 0 respectively.Note that p pk is the only feature with units: specifically, it is scaled to m p c.
The PLDs for the 15 parameters in Equations 9-13 were estimated via MCMC as described in Appendix A. The results are presented in Table 1.The estimates suggest that p pk depends only weakly on B n 0 , since p is close to 0. Similarly, L and H have a very weak dependence on both B n 0 and η mfp .Therefore the parameter estimation procedure was repeated with a reduced 10-parameter expression, making the following modifications: Parameter estimates for the reduced 10-parameter expression also appear in Table 1.
Parameter 15-parameter 10-parameter p C 7.88×10 The parameter estimates in Table 1 suggest that when electron energy is limited by the location of the FEB, the value of p max scales as roughly (Γβ) 2 ; once radiation losses are the limiting factor, p max scales as approximately (Γβ) −1 .Neither the value of L nor that of H depends strongly on the acceleration efficiency (η mfp ) or the peak magnetic field strength ( B n 0 ), motivating the reduction from a 15-parameter set to just 10.The location of the break in the power law, Γ pk β pk , decreases as the magnetic field strength goes up, and increases as acceleration efficiency drops.The height of the break, p pk , decreases as acceleration efficiency drops, and shows almost no dependence on the magnetic field strength.As for the width parameter W , it increases with larger values of B n 0 and with smaller values of η mfp ; since larger W means a sharper peak, however, the trend is for broader breaks in the power law when either the magnetic field is weak or acceleration is inefficient.
With the values presented in Table 1, it is possible to compute the value of p max for any desired φ MC .First η mfp and B n 0 are used in Equations 9-13 to compute the five features of Equation 8.The computed features are, in turn, used to obtain p max by evaluating Equation 8 at the desired value of Γβ.

DISCUSSION
We now consider two applications of Equation 8to GRB afterglows.The first is the upper limit on synchrotron photon energy produced by the external shock of a GRB, and the second is a brief discussion of the physical process that limits electron energy in our simulations.

Maximum energy of synchrotron photons
GRB afterglows are routinely associated with emission above 100 MeV, and analytical treatments of particle acceleration suggest that electron synchrotron photons can exceed ∼GeV energies while the forward shock of the GRB is still moving ultrarelativistically (Piran & Nakar 2010;Kumar et al. 2012).The MC simulations presented in Section 3 serve as a numerical check on previous analytical calculations, as they take place under largely the same assumed conditions (c.f.Section 2).We now calculate hν syn,max = 3γ 2 max hqB/(4πm e c) for electron distributions whose maximum momentum p max is given by Equation 8.This is the characteristic synchrotron frequency of such electrons; the resultant photon spectrum would peak at a frequency approximately 0.3ν syn,max (Rybicki & Lightman 1979).
Let us now restrict our focus only to the portion of Equation 8 above p pk , where radiative losses rather than escape act to limit the maximum electron energy.At shock Lorentz factors far above Γ pk β pk , Equation 8 simplifies to where we have reduced the independent variable to Γ rather than Γβ.Using the 10-parameter values from Table 1, the above equation can be expressed having recast the equation in terms of Lorentz factor rather than momentum.(The feature W cannot be evaluated without knowing the value of B n 0 ; the dependence on η mfp is weak enough to be ignorable.For the range of conditions simulated in this work W varies between 1.73 and 4.08.We have taken 2.90 as a fiducial value to eliminate the dependence of γ max on W .The eventual result for hν syn,max differs by a factor close to unity as a consequence of this simplification.) We now assume that the large-scale hydrodynamics are governed by the Blandford & McKee (1976) solution for a relativistic impulsive explosion.Further assuming an adiabatic shock, with neither energy injection nor significant radiative losses, one may derive the following expression for the shock Lorentz factor Γ: In this equation, k controls the density dependence of the circumburst medium (0 for constant density, and 2 for a wind-like medium); E iso is the isotropic-equivalent explosion energy of the GRB; N is the density parameter (such that the number density encountered by the shock is n 0 = N R −k at some distance R); and t obs is the time since the start of the GRB in the observer frame.(Cosmological redshift may be significant, but we omit the factors of (1 + z) for the sake of simplicity.)A similar application of the Blandford-McKee solution results in an expression for the ambient density n 0 : The characteristic frequency of synchrotron emission depends on the electron Lorentz factor, γ max as shown in Equation 18, and on the magnetic field strength B as defined by Equation 3. We can specify k = 0 or k = 2 in these two equations, replacing N with n ism or A as appropriate: where we have used the notation Q = Q x ×10 x , e.g.A = A ,34 ×10 34 cm −1 , and all quantities use cgs units.We remind the reader that the observed peak of the photon distribution would occur at energies roughly 30% of those given in Equation 24.
In the k = 0 case, the exponents of energy and time suggest that they may be due to uncertainty in the parameter H C , and that the value of hν syn,max may not depend on either quantity: if H C were precisely 1 instead of 1.01, neither factor would appear in Equation 24.That is, the results presented in Figure 2 suggest that the maximum observed synchrotron energy is essentially independent of both explosion energy and time.For the k = 2 case, the maximum synchrotron frequency grows with time; though both Doppler factor and peak magnetic field decay with time, they are more than compensated for by the growth in maximum electron Lorentz factor.
The above limits should hold as long as (1) the motion of the blast wave is given by the adiabatic Blandford-McKee solution, (2) the shock Lorentz factor is above Γ pk β pk given by Equation 10, and (3) there is amplified microturbulence that decays on either side of the shock as prescribed by Equation 4. Condition (1) means that Equation 24may not apply very early on in the afterglow, if the shock is still in the acceleration or coasting phases of its expansion.Condition (2) is satisfied for Γβ > ∼ 10.Using this limitation in Equation 19 leads to the following restriction on t obs : As shown above, both ISM and wind media allow for electron synchrotron photons that reach GeV energies; but the extremely high energy photons detected in events like GRB 090902B (33.4 GeV; Abdo et al. 2009), GRB 130427 (95 GeV;Ackermann et al. 2014), and GRB 190114C (> 300 GeV;Mirzoyan 2019) suggest that an additional channel produces the highest-energy photons detected in GRB afterglows.(Alternately, electron acceleration could proceed in a more exotic scenario than modeled in this work, e.g.Khangulyan et al. 2020, who consider acceleration in the presence of clumpy magnetic field condensations.)Though thermal electrons were ignored in this work, the synchrotron self-Compton mechanism in this population is a straightforward way to produce emission that peaks in the GeV-TeV range well into the afterglow (Warren et al. 2017).

What process limits electron energy?
The idea of using the synchrotron cooling break to interpret afterglow observations is nearly as old as afterglow observations themselves (Sari et al. 1998); for a recent example of a work using this formalism, see Fraija et al. (2020), or the numerous citations within.There are multiple approaches to determining the location of the cooling break.One method is to assume that the electron distribution immediately behind the shock extends to infinity (Granot & Sari 2002); synchrotron losses then quickly cool the highest-energy electrons and limit the extent of the electron distribution downstream from the shock.Another, more realistic, way is to equate the electron acceleration time to the synchrotron cooling time to determine the energy at which acceleration stops (Lemoine 2013;Khangulyan et al. 2020).
None of the previous analytical methods would seem to apply to the acceleration scenario presented in this work.For example, consider the simplest case: electrons scattering in a homogeneous magnetic field.In this setting, the synchrotron loss time is given by for electrons with energy γ elec moving in a homogeneous field of strength B. The maximal electron energy may be calculated by comparing the above loss time to the typical residence time of an electron in the downstream region of the shock.In keeping with the Bohm approximation above, we assume that the residence time is proportional both to the gyroperiod of the electron and to the Bohm factor: For electrons with the maximum achievable energy, t syn and t res|d should be roughly equal, yielding The magnetic field is determined by Equation 6. Substituting that equation into Equation 28 leads to the proportionality Equation 29 has a different dependence on Γ 0 (−1/2) than displayed in Figure 2 and Table 1 (≈ −1).This suggests that the above scenario does not reflect the energy-limiting process occurring in the MC simulations.This is unsurprising, as Equations 26 and 27 assumed a homogenous magnetic field of strength B, while the MC simulations took place in a decaying microturbulent magnetic field.Another straightforward limit on the maximum attainable energy of a shock-accelerated particle is to equate the acceleration time to the dynamical time of the shock.The former quantity is conservatively given by the downstream residence time in Equation 27.The latter is proportional to t shock /Γ for a shock of age t shock .Setting the two to be roughly equal leads to the proportionality where we have substituted Equation 6 into Equation 27.This approach also fails to reproduce the dependence on Γ seen in Figure 2 and Table 1, suggesting either that the assumptions made are not correct or that the maximum electron energy is not limited by time.
The broad structure visible in Figure 2 is that of a broken power law, which suggests that two processes control p max in our model of electron acceleration.There are hints, however, of more structure-particularly at the top right corner (high η mfp and low B n 0 ), where the empirical equation seems to slightly overestimate p max , and at the bottom left corner (low η mfp and high B n 0 ), where there may be an additional power law segment at higher Γβ values.These possible additional features could be physical in origin, something we intend to investigate in future work.

CONCLUSIONS
We have presented here the results of numerous Monte Carlo (MC) simulations of diffusive shock acceleration of electrons.Acceleration took place in a highly turbulent magnetic field that decayed upstream and downstream from the shock, as observed in particle-in-cell simulations of relativistic shocks.The accelerated electron spectra are well-described by power laws with an exponential turnover (Equation 1, Figure 1), and the likelihood distributions of their parameters π pl = (p max , σ, χ, log 10 K) were estimated using a Markov chain Monte Carlo (MCMC) approach.We focused here only on the parameter p max , i.e. the location of the exponential turnover and the energy at which acceleration gains balance loss processes.When the estimates of p max are plotted across a 3-dimensional parameter space of astrophysical interest, they appear to form broken power laws as a function of the shock speed Γβ (Figure 2).The peak in each curve lies close to the shock speed at which electrons cease to escape upstream from the shock, which suggests that upstream escape and radiation losses are the two limiting processes that control the maximum electron energy.
Based on the broken power law observed for p max as a function Γβ, we proposed an empirical expression in terms of physically-motivated features (Equation 8).Another round of MCMC was used to estimate its parameters (Table 1).The resulting expression allows for prediction of the peak electron energy-and correspondingly the synchrotron cooling energy-in scenarios where relativistic shocks accelerate electrons, without the need for time-consuming numerical simulations.Extensions to this work will involve increasing the size and dimensionality of the MC-simulated parameter space φ MC = (η mfp , B n 0 , Γβ) from what was considered here, as well as quantitative treatments of the σ and χ parameters in Equation 1.

Figure 1 .
Figure1.Comparison between the electron spectra obtained from the MC code, and the spectrum described by Equation 1 for the MCMC-identified highest-likelihood π pl (see Section 4).Histograms of electron spectra from 20 different MC iterations (using the same π pl but different random number seeds) are shown with light grey lines.Their geometric average is plotted in red.The black curve is Equation 1 using the highest-likelihood parameters identified via MCMC.The input parameters for the MC simulation were φ MC = (Γβ = 100, η mfp = 1, B n 0 = 10 −4 ).The corresponding highest-likelihood parameters estimated by MCMC were π pl = (log 10 p max = 5.13, σ = 2.18, χ = 1.86, log 10 K = 60.7).

Figure 2 .
Figure 2. Dependence of Equation 1's p max on Γβ for the various combinations of B n 0 and η mfp simulated.The curves correspond to the p max values predicted by the empirical expression (Equations 8-13) in its full (15-parameter, dashed line) and reduced (10-parameter, solid line) forms; the values of the parameters are provided in Table1below.Gray shading marks regions of the parameter space where electrons were able to escape from the shock in the upstream direction.The rightmost column compares MC runs with the same B n 0 but different η mfp .The bottom row compares MC runs with the same η mfp but different B n 0 .

Table 1 .
The mode of each parameter's PLD; and the boundaries of the 68% credible region of the PLD for each parameter, marginalized over all others.For the five fixed parameters of the 10-parameter set, no credible region is given.