*Article* **Particle Acceleration in Mildly Relativistic Outflows of Fast Energetic Transient Sources**

**Andrei Bykov \*, Vadim Romansky and Sergei Osipov**

Ioffe Institute, 194021 St. Petersburg, Russia; romanskyvadim@gmail.com (V.R.); osm2004@mail.ru (S.O.) **\*** Correspondence: byk@astro.ioffe.ru

**Abstract:** Recent discovery of fast blue optical transients (FBOTs)—a new class of energetic transient sources—can shed light on the long-standing problem of supernova—long gamma-ray burst connections. A distinctive feature of such objects is the presence of modestly relativistic outflows which place them in between the non-relativistic and relativistic supernovae-related events. Here we present the results of kinetic particle-in-cell and Monte Carlo simulations of particle acceleration and magnetic field amplification by shocks with the velocities in the interval between 0.1 and 0.7 c. These simulations are needed for the interpretation of the observed broad band radiation of FBOTs. Their fast, mildly to moderately relativistic outflows may efficiently accelerate relativistic particles. With particle-in-cell simulations we demonstrate that synchrotron radiation of accelerated relativistic electrons in the shock downstream may fit the observed radio fluxes. At longer timescales, well beyond those reachable within a particle-in-cell approach, our nonlinear Monte Carlo model predicts that protons and nuclei can be accelerated to petaelectronvolt (PeV) energies. Therefore, such fast and energetic transient sources can contribute to galactic populations of high energy cosmic rays.

**Keywords:** fast blue optical transients; non-thermal particle acceleration; particle-in-cell plasma modeling; high energy cosmic rays

#### **1. Introduction**

The time domain astronomy operating now at all wavebands from radio to gamma rays has provided unique information on highly energetic processes in transient astrophysical objects such as supernovae (SNe), gamma-ray bursts (GRBs), fast radio bursts and others. The pioneer program of a dedicated search for supernovae in the optical band started by Fritz Zwicky in 1936, which has later allowed, in particular, obtaining fundamental results on the accelerated expansion of the Universe via spectroscopy of SN type Ia, is now ongoing with great perspectives [1]. The capabilities of fast and sensitive wide field imaging suggested for the forthcoming Large Synoptic Survey Telescope (LSST) would allow detecting many thousands of luminous SNe and tidal disruption events (TDEs) per year as well as studying other types of transient sources [2]. Together with the LSST, the currently operating wide field [3] and survey [4] X-ray observatories along with the future high energy missions [5–7], gravitational wave and neutrino observatories will allow revealing the physical nature of various types of energetic space transients.

The fast blue optical transients (FBOTs) [8–11] are among the most interesting recent discoveries. Their appearance is somewhat different from most of the core-collapse SNe [8,12]. Together with the low-luminosity GRBs they possibly belong to the intermediate class of phenomena filling the gap between non-relativistic SNe and "standard" long duration gammaray bursts, and which could have volumetric rates well above that of the GRBs [13,14]. The duration of both the energy-momentum release from the central engine and the interaction of the anisotropic ejecta with the outer layers of the progenitor star and its circumstellar matter determine the transient appearance (see, e.g., [15]).

Three recently studied powerful FBOT sources AT2018cow [9], CSS161010 [10] and ZTF18abvkwla [11] were characterized by a low ejected mass and fast outflows. Indeed,

**Citation:** Bykov, A.; Romansky, V.; Osipov, S. Particle Acceleration in Mildly Relativistic Outflows of Fast Energetic Transient Sources. *Universe* **2022** , *8* , 32. https://doi.org/10.3390/ universe8010032

Academic Editors: Galina L. Klimchitskaya, Vladimir M. Mostepanenko and Nazar R. Ikhsanov

Received: 19 November 2021 Accepted: 3 January 2022 Published: 5 January 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

97

Margutti et al. [9] found that to explain the fast rise of the optical and radio emission together with the persistent photosphere appearance in AT2018cow a wide range of velocities in the range from about 0.02 c to 0.2 c is needed, where *c* is the speed of light. The aspherical ejecta with the range of velocities had the estimated mass ∼0.1–1 *M* [9]. The long low frequency *uGMRT* radio observations [16] allowed estimating the shock radius to be *<sup>R</sup>* = (6.1 − 14.4) × 1016 cm and the speed of the fast ejecta to be above 0.2 c at 257 days after the shock breakout. The mass loss rate of the progenitor star was found to be a few times 10−<sup>6</sup> *<sup>M</sup>* yr−<sup>1</sup> during the period of 20–50 years assuming the wind velocity ∼1000 km s<sup>−</sup>1, while it was possibly 100 times higher in a period of a few years just before the event [16].

Another interesting FBOT source is CSS161010 located in a dwarf galaxy at a distance about 150 Mpc [10]. On the basis of the synchrotron interpretation of its radio emission, which is peaked at about 100 days after the FBOT event, the authors suggested a presence of a mildly relativistic outflow of four-velocity 0.55 c driving a blast wave. They estimated the ejected mass to be in the range of 0.01–0.1 *M*. The outflow is faster than that estimated for AT2018cow and is similar to that in ZTF18abvkwla. The origin of the bright X-ray luminosity is attributed to an emission component, which is likely different from the primary one, which produced the synchrotron radio emission. Recent radio, millimeter wave and X-ray observations [17,18] of a short-duration luminous FBOT transient ZTF20acigmel (AT2020xnd) located at z = 0.2433 indicated a presence of a fast ejecta with a ∼0.2 c speed shock with estimated energy above 1049 ergs. AT2020xnd has shown high radio luminosity of *<sup>L</sup><sup>ν</sup>* ∼ 1030 ergs−<sup>1</sup> Hz−<sup>1</sup> at 20 GHz almost 75 days after the event [18]. The observational data suggested a shock driven by a fast outflow of velocity 0.1–0.2 c interacting with the dense circumstellar matter shaped by an intense wind of *<sup>M</sup>*˙ <sup>≈</sup> <sup>10</sup>−<sup>3</sup> *<sup>M</sup>* yr−<sup>1</sup> with velocity of *vw* = 1000 km s−<sup>1</sup> from the progenitor star [18]; the presence of a steep density profile of *ρ*(*r*) ∝ *r*−<sup>3</sup> in the wind was suggested by [17]. Similar to AT2018cow, the detected X-ray emission is in excess compared to the extrapolated synchrotron spectrum and constitutes a different emission component, possibly powered by accretion onto a newly formed black hole or neutron star.

The distinctive features of the four FBOT transients discovered so far are their high peak bolometric luminosity *<sup>L</sup>* <sup>&</sup>gt;<sup>∼</sup> <sup>1043</sup> erg s−<sup>1</sup> and a rapid timescale of a few days duration. Multi-wavelength observations of these objects uncovered the presence of powerful subor mildly relativistic outflows which are likely originated from either rare SN-type or TDE-type events with intermediate or stellar mass black holes (see, e.g., [19]). The kilonova type sources with the neutron stars mergers typically eject the masses in the range of <sup>10</sup><sup>−</sup>4–10−<sup>2</sup> *<sup>M</sup>* (and even <sup>∼</sup>0.1 *<sup>M</sup>* for the black hole—neutron star mergers) with velocities ∼0.1–0.3 c [20]. It is interesting that the recent model of a supernova from a primordial population III star of a 55,500 *M* mass with general relativistic instability [21] predicts the ejecta velocities of about 0.3 c. Earlier mildly and moderately relativistic ejecta outflows were found in a few broad line type Ic SNe (e.g., [22,23]). The geometry and structure of the outflows producing FBOTs depend on the source of their power and it is a subject of modeling [9]. Relativistic mass ejection in spherically symmetric shock outflows of core-collapse supernovae was studied in detail in [24–26]. Their high-velocity solutions demonstrated rather a steep dependence of the deposited kinetic energy *Ek* as a function of the ejecta four-speed *Ek* ∝ (*β*Γ)−5.2. Much flatter energy—ejecta velocity distribution *Ek* ∝ (*β*Γ)−2.4 can be obtained for engine-driven asymmetric supernovae with a powerful activity of compact stellar remnants [27–30]. Recently, numerical models of supernova explosions where the supernova ejecta interacts with the relativistic wind from the central engine were constructed in [31]. The formation of relativistic flows in the interaction of the powerful non-thermal radiation produced by the central machine with the supernova ejecta was considered in the papers [32,33].

Very luminous optical FBOT events with light curves extending to a couple of weeks can be expected in the case of shock breakout into a dense circumstellar shell produced by the dense progenitor wind a few years before the SN event [34–36]. The presence of an

hour timescale central engine activity with a luminosity of about 10<sup>47</sup> erg s−<sup>1</sup> producing mildly relativistic jet was proposed in [15] to model SN 2006aj associated with a lowluminosity GRB. The bright high frequency radio emission of FBOTs and, possibly, the hard X-ray component detected in some events of this type, is produced by relativistic electrons accelerated by shocks driven by fast moderately relativistic outflows from the central engine. Moreover, the synchrotron self-absorption effects [37] are apparent in some FBOT spectra.

The SN shock breakout is usually accompanied by a bright ultraviolet (UV) flash, and it is likely that afterwards the shock enters a collisionless regime [38] with the Xray dominated spectrum. Some models of relativistic shock breakout which consider a multifluid structure of a relativistic shock mediated by radiation in a cold electron-proton plasma are currently under discussion [39,40]. A MeV gamma-ray flash of the total energy ∼10<sup>48</sup> erg lasting from a few seconds to a few hours was predicted in [41,42] for some SN events. At the later SN stages (after a few days from the event) the radiative shock is transforming into a collisionless plasma shock regulated by kinetic plasma instabilities. A specific feature of the collisionless shocks is their ability to create a powerful non-thermal particle population and to accelerate relativistic particles [43]. The collisionless shock structure and the efficiency of particle acceleration depend on the shock speed, on magnetic field inclination to the shock normal and on the plasma magnetization parameter [44,45]. Modeling of particle acceleration by non-relativistic shocks of velocities below 0.1 c in supernova remnants was discussed in [46] while studies of relativistic shocks was presented in [47].

Kinetic simulations of the efficiencies of the shock ram pressure conversion to magnetic fluctuations and relativistic particles are needed to provide an adequate interpretation of the observed non-thermal radiation. The mildly relativistic magneto-hydrodynamic (MHD) flows were shown (see, e.g., [48,49]) to be the most efficient environment providing the maximum energies of the accelerated nuclei for a given magnetic/kinetic luminosity of the power engine. Recent models suggested a possibility of cosmic ray acceleration to ultra high energies in the low-luminosity GRBs associated with SNe (see, e.g., [50,51]) and in relativistic SNe [52]. Fast outflows from SNe with dense circumstellar shells could accelerate cosmic rays up to the high energy regime on a few weeks timescale (e.g., [53,54]).

The structure and particle acceleration in the fast collisionless shocks can be successfully modeled with the kinetic particle-in-cell (PIC) technique [43,55] at many thousands of the particle gyro-scales. On the other hand, as the observed multi-wavelength spectra of supernova remnants and GRBs compellingly demonstrated that the spectra of accelerated particles extend to many decades in the particle momentum, a combination of both microscopic (see Section 2) and macroscopic kinetic models (see Section 3) is necessary to construct realistic models of such sources. Cosmic ray acceleration by supernova remnants is a subject of extensive modeling [46,52,56–59]. One of the most uncertain points is at what stage PeV regime cosmic rays can be accelerated. We discuss in Section 3 high energy cosmic ray acceleration in FBOT-type sources as potential pevatrons.

Here we will present 2D PIC simulation of sub and mildly relativistic shocks within the shock speed interval 0.1–0.7 c in proton-electron plasmas which could be applied to FBOT sources modeling. The fiducial case in our work is 0.3 c.

Earlier, Park et al. [60] presented a set of 1D PIC models of shocks in a range of velocities up to 0.1 c. They performed in the 1D case long simulation runs up to 5 × 105 *<sup>ω</sup>*−<sup>1</sup> *pe* . For higher shock speed of 0.75 c Crumley et al. [61,62] published 2D PIC simulation results. In the quasi-parallel shock of a velocity 0.75 c simulation of about 4100 *ω*−<sup>1</sup> *pi* duration allowed them to model Bell-mediated shock. In particular, they noted that acceleration of electrons is likely initially associated with the shock drift acceleration and then switches to diffusive shock acceleration regime as it was seen in [60].

To model radio emission observed from CSS161010, where the shock velocity of about 0.3 c was suggested by observations, we make 2D PIC simulations with *mp*/*me* ratio of 100, which are limited to the timescale ∼<sup>7</sup> × 104 *<sup>ω</sup>*−<sup>1</sup> *pe* . We discuss possible extrapolations of electron spectra to the energy range needed to model the radio emission.

#### **2. Particle-in-Cell Simulations of Fast Mildly Relativistic Shocks**

In this section, we present the results of particle-in-cell simulations of mildly relativistic shocks at gyro-scales with application to the observed non-thermal emission from FBOTs.

To simulate the structure and non-thermal particle acceleration by a mildly relativistic collisionless shock wave we have employed the publicly available code Smilei developed by Derouillat et al. [63]. This code is based on explicit Finite Differences Time Domain approach for solving Maxwell equations and on a relativistic solver for particle movement with charge-conserving algorithm proposed by [64].

For model shock initialization we employ a common approach: a homogeneous plasma flow collides with an ideal reflecting wall. The simulation box is two-dimensional, with the number of cells *Nx* = 204,800 in the direction along the flow velocity and *Ny* = 400 in the transverse direction, while particle velocities and the electromagnetic field are represented by full 3D vectors. Homogeneous plasma flows in through the right boundary and the reflecting wall is placed at the left boundary. Boundary conditions in the transverse direction are periodic. The initial magnetic field - *B* lies in the plane of the simulation inclined by the angle *θ* to the velocity of the flow. The electric field is initialized to compensate the Lorentz force in the laboratory frame - *E* = −*<sup>v</sup>* <sup>×</sup> - *B*/*c*. The velocity of plasma flow *v* is 0.3 c and its Lorentz-factor *γ* is 1.05. The magnetization *σ* = *B*2/4*πnγmpc*2, where *n* is upstream concentration and *mp* is mass of proton, in all the modeled configurations is about 10<sup>−</sup>4. The electron mass *me* is artificially increased up to *mp*/*me* = 100 in order to save computational resources. All the quantities obtained from the simulation can be scaled with respect to the plasma concentration, which can be chosen arbitrary during the data analysis. The spatial grid step is *dx* = 0.2*c*/*ω<sup>e</sup>* and the time step is *dt* = 0.09 *ω*−<sup>1</sup> *<sup>e</sup>* , where *ω<sup>e</sup>* is electron plasma frequency *ω<sup>e</sup>* = 4*πne*2/*meγ*, *e* is the absolute value of the electron charge. Within such a setting the simulation box size along the *x*-axis corresponds to 500 gyroradii of protons in plasma flow *rg* = *mpvcγ*/*eB* and along the *y*-axis—to 1 gyroradius. The maximum simulation time is 7 <sup>×</sup> <sup>10</sup><sup>4</sup> *<sup>ω</sup>*−<sup>1</sup> *<sup>e</sup>* or about 100 inverse proton gyrofrequencies. For plasma concentration *<sup>n</sup>* ≈ 1 cm−3, this corresponds to timescale about 1 s, which is much smaller than the typical activity period of FBOTs.

The efficiency of particle acceleration depends on the inclination angle *θ*, especially in the case of relativistic flows [44,61,62,65]. To participate in the diffusive acceleration process, a particle needs to escape from the shock front and if it moves along the magnetic field, the maximum velocity along the *x*-axis is *c* cos(*θ* ), and it should be larger than the shock velocity *v sh* (all quantities here are measured in the upstream frame). As one can see in Figure 1, the high energy tail of the electron distribution is much higher for a quasiparallel shock. Such an angle dependence may lead to the presence of different electron distributions within one object, and this may possibly explain the difference between the spectral indices of synchrotron radio and inverse Compton X-ray radiation observed in FBOT AT2018cow [9].

For further modeling of synchrotron radiation from FBOT CSS161010 we have used a setup with parameters *v* = 0.3 c, *σ* = 0.0002 and *θ* = 30◦ (a quasiparallel shock). Time evolution of the concentration profile in such a shock averaged in the transverse direction is shown in Figure <sup>2</sup> and the corresponding magnetic field at the moment 7 <sup>×</sup> <sup>10</sup><sup>4</sup> *<sup>ω</sup>*−<sup>1</sup> *<sup>e</sup>* is shown in Figure 3.

With the computed concentration profile we can determine the coordinate of the shock front: we consider that *xsh* is the first point from the right where the concentration is two times larger than that in the far upstream. The average shock velocity is the ratio of the shock coordinate and the simulation time *vsh* = *xsh*/*t*. At later times the shock wave propagation is close to a stationary regime. The shock velocity measured in the downstream frame is close to the constant value *vsh* ≈ 0.085 c. In the upstream (observer) frame it corresponds to *v sh* = (*vsh* + *v*)/ 1 + *vshv*/*c*<sup>2</sup> ≈ 0.38 c.

**Figure 1.** Electron distribution function in the shock downstream with initial parameters *v* = 0.3 c, *σ* = 0.0002 and different inclination angles.

**Figure 2.** Time evolution of concentration normalized to the far upstream concentration.

**Figure 3.** Magnetic field normalized to the far upstream magnetic field.

The simulations have allowed us to compute model distributions of electrons. In the downstream region such distributions have a complex shape, where one can see two main components: a thermal peak and a power law tail at high energies. At low energies (*E* < 5 *mec* 2) we approximate the distribution with a Maxwell-Juttner function: an iterative process is employed to minimize the functional *<sup>f</sup>*(*T*) = <sup>5</sup>*mec*<sup>2</sup> *mec*<sup>2</sup> (*F*(*E*) <sup>−</sup> *Fmj*(*E*, *<sup>T</sup>*))2*dE*, where *<sup>F</sup>*(*E*) is simulated electron distribution function and *Fmj*(*E*, *T*) is Maxwell-Juttner distribution function, and find the effective temperature. For high energies (20 *mec* <sup>2</sup> < *E* < 50 *mec* 2) we follow a least squares approach for linear regression in double logarithmic coordinates and obtain the power law spectral index. The electron spectrum in a close downstream of the shock (5000 grid cells behind the shock) at time *t* = 70,000 *ω*−<sup>1</sup> *<sup>e</sup>* for the setup with *σ* = 0.0002 and *<sup>θ</sup>* <sup>=</sup> <sup>30</sup>◦ and its approximation with temperature *Te* <sup>=</sup> <sup>5</sup> <sup>×</sup> 1010 K and spectral index *<sup>s</sup>* <sup>=</sup> 3.59 are shown in Figure 4.

Time dependence of electron distribution function at later stages of simulation is shown in Figure 5. Distribution is not stationary and one can see the oscillation of the electron spectral distribution due to the influence of magneto-hydrodynamic instabilities, as described in [61,66]. So we have chosen the electron distribution at time *t* = 70,000 *ω*−<sup>1</sup> *<sup>e</sup>* for the further modeling.

Recent observations of FBOTs [8–11] discussed above have shown that their synchrotron spectra are strongly influenced by synchrotron self-absorption. At the given time moment, observable radio fluxes rise as *F<sup>ν</sup>* ∝ *ν*5/2 for low frequencies, then have a peak and fade with a power law tail, which depends on the particular electron distribution function. Additionally, the time dependence of such spectrum is very specific and its details could be used to imply a number of source parameters: the magnetic field and the shock radius can be determined from the measured maximum of the light-curve *Fmax* at given frequency

*ν*, once the fractions of energy in accelerated electrons and magnetic field *<sup>e</sup>* and *<sup>B</sup>* are established, with the relations derived by Chevalier [37].

**Figure 4.** The electron distribution function in the shock downstream with initial parameters *v* = 0.3 c, *<sup>σ</sup>* <sup>=</sup> <sup>2</sup> <sup>×</sup> <sup>10</sup>−<sup>4</sup> and *<sup>θ</sup>* <sup>=</sup> <sup>30</sup>◦.

**Figure 5.** The electron distribution function in the shock downstream at different time moments.

**]**

**Log10 E [mec**

$$R = \left(\frac{6\epsilon\_B c\_6 \, ^{s+5}F\_{max} \, ^{s+6}D^{2s+12}}{\epsilon\_\ell f(s-2)\pi^{s+5}c\_5 \, ^{s+6}E\_1 \, ^{s-2}}\right)^{\frac{1}{2s+13}} \frac{2c\_1}{\nu} \, ^{\prime} \tag{1}$$

$$B = \left(\frac{\varepsilon\_B 2 \text{36} \pi^3 c\_5}{\varepsilon\_\varepsilon^2 f^2 (s-2)^2 c\_6 ^3 E\_1 ^{2s-4} F\_{\text{max}} D^2}\right)^{\frac{2}{2s+15}} \frac{\nu}{2c\_1} \tag{2}$$

where *s* is the electron spectral index, derived from the power law tail of emissivity spectrum and its time dependence, *c*<sup>1</sup> , *c*5, *c*<sup>6</sup> are the Pacholczyk constants [67], depending on *s*, *E*<sup>1</sup> is the minimum energy of the electron power law distribution, which usually equals to the electron rest energy, *f* is the filling factor—the emitting fraction of the source volume, and *D* is the distance to the source. A problem of this method is that it needs a lot of assumptions about the shock structure and electron distribution. In this model the emitting region is considered to be a homogeneous flat disk with the radius *r* and depth *f* × *r*. There is also a model modification for a spherical source with an account for source inhomogenuity [17]. The electron distribution is considered to be a power law. Additionally, fractions of energy *<sup>e</sup>* and *<sup>B</sup>* are unknown and are often chosen according to the equipartition rule *<sup>e</sup>* = *<sup>B</sup>* = 1/3. However the simulation results show, that these values are unlikely to be that large. Following Chevalier [68] we define the energy fractions in terms of the upstream kinetic energy density (*mp* + *me*)*nuv sh* 2 , where *nu* is the upstream concentration. In this notation *<sup>e</sup>* = *Ee*/*ρuv sh* <sup>2</sup> and *<sup>B</sup>* <sup>=</sup> *<sup>B</sup>*2/ 8*π*(*mp* + *me*)*nuv sh* 2 , *Ee* is the accelerated electron energy evaluated by subtracting from the total electron energy, that corresponds to the Maxwell–Juttner fit of the distribution function. The values obtained by particle-in-cell simulations in our fiducial setup are *<sup>e</sup>* = 0.014 and *<sup>B</sup>* = 0.03. They depend on the initial conditions, but for a wide variety of parameters we see that the fraction of energy in magnetic field is lower than 10%. Particle-in-cell simulations have rather small scales and cannot describe the influence of long wave upstream instabilities caused by high energy particles, but Monte Carlo simulations show similar values, as described below.

Using the electron distribution function obtained from the PIC simulation we can evaluate the spectral density of the energy flux of synchrotron radiation from the source, taking into account the effect of synchrotron self absorption. Standard formulae for that effect are described in detail in [69]. Emitted power per unit frequency per unit volume is

$$I(\nu) = \int\_{E\_{\min}}^{E\_{\max}} dE \frac{\sqrt{3}e^3 nF(E)B\sin(\phi)}{m\_\varepsilon c^2} \frac{\nu}{\nu\_\varepsilon} \int\_{\frac{\nu}{\nu\_\varepsilon}}^{\infty} \mathcal{K}\_{5/3}(\mathbf{x})d\mathbf{x},\tag{3}$$

where *φ* is the angle between the magnetic field and the line of sight, *ν<sup>c</sup>* is the critical frequency *ν<sup>c</sup>* = 3*e* <sup>2</sup>*B* sin(*φ*)*E*2/4*πme* <sup>3</sup>*c*5, and *K*5/3 is the modified Bessel function. The absorption coefficient is

$$k(\nu) = \int\_{E\_{\rm min}}^{E\_{\rm max}} dE \frac{\sqrt{3}e^3}{8\pi m\_\varepsilon \nu^2} \frac{nB\sin(\phi)}{E^2} \frac{d}{dE} E^2 F(E) \frac{\nu}{\nu\_\varepsilon} \int\_{\frac{\nu}{\nu\_\varepsilon}}^{\infty} \mathcal{K}\_{5/3}(\mathbf{x}) d\mathbf{x}.\tag{4}$$

The obtained spectral density is further integrated over the volume of the source, which is described as a spherical shell, whose volume is determined by a filling factor *f* = 0.5. Hence the total emissivity and observable flux at distance *D* can be derived. The magnetic field *B* is considered constant and perpendicular to the line of sight. The concentration in the stellar wind depends on radius as *n* ∝ *r*<sup>−</sup>2, but in the downstream of the shock we assume it constant. The electron distribution function is also constant in the volume of the source. This modeling is further applied to explain the spectrum of the FBOT transient CSS161010. As described in [10] at time *t* = 357 days after explosion its shock velocity was *vsh* = 0.36*c*, which corresponds to the modeled upstream plasma flow velocity *vs*. = 0.3*c*, measured in the downstream frame. At this moment, parameters derived with Equations (1) and (2), assuming the equipartition regime *<sup>e</sup>* = *<sup>B</sup>* = 1/3, are as follows: the magnetic field *<sup>B</sup>* <sup>=</sup> 0.052 G, the outer radius *<sup>R</sup>* <sup>=</sup> 3.3 <sup>×</sup> 1017 cm, the concentration *n* = 1.9 cm−3, the electron spectral index *s* = 3.5, which is similar to that obtained from a PIC simulation, and the minimum energy *E*<sup>1</sup> = 4 *mec*2.

The thermal electrons may contribute substantially to the observed synchrotron spectrum and it was shown in [70] that models with the power law electron distribution might be oversimplified. Additionally, equipartition regime is not very reasonable assumption. Thus, the use of the realistic kinetic simulations is important to model broad band radiation spectra of FBOTs. We have compared the observed fluxes to three models of source emission, evaluated with power law electron distribution and two distributions obtained from PIC simulation: with extrapolation of power law tail to higher energies (*Emax* = 500 *mec*2) and without it. For the first case, we have used the parameters described above. In order to evaluate parameters *B* and *R* of the simulated particle distribution, we minimized the functional *<sup>g</sup>*(*B*, *<sup>R</sup>*) = ∑(*F*(*νi*, *<sup>B</sup>*, *<sup>R</sup>*) <sup>−</sup> *Fobs*(*νi*))2, where *<sup>ν</sup><sup>i</sup>* are the observed frequencies, *Fobs*(*νi*) are the observed fluxes and *F*(*νi*, *B*, *R*)—the modeled fluxes, using the gradient descent algorithm. We employed six measurements from [10]: four made with VLA at day 357 at frequencies 1.5, 3.0, 6.05 and 10 GHz and two made by GMRT at day 350 at frequencies 0.33 and 0.61 GHz. Concentration in these equations is determined by the magnetic field *B*, obtained from the PIC simulation. One can see that distribution, obtained directly from particle-in-cell simulation, cannot correctly fit the power law tail in observational data, so we had to extrapolate the simulated electron distribution to higher energies because the PIC approach requires a lot of computational resources and thus is not suitable to simulate the considered system up to timescales long enough to form a long tail of accelerated particles. The modeled parameters for extrapolated distribution are *<sup>B</sup>* <sup>=</sup> 0.069 G, *<sup>R</sup>* <sup>=</sup> 3.0 <sup>×</sup> 1017 cm and concentration *n* = 210 cm−3. The values of the magnetic field and radius are rather close to the values from [10], while the concentration is much higher. These results are illustrated in Figure 6.

**Figure 6.** Observed (green circles) and modeled spectral energy distribution of CSS161010 at day 357 for various model electron distributions.

One can see that the parameters of the shock (especially the concentration), obtained from the radiation model strongly depend on the electron distribution function. This should be kept in mind during interpretation of the observational data.

#### **3. A Monte Carlo Model of Cosmic Ray Acceleration in Fast Transient Sources**

We have developed a Monte Carlo model of particle acceleration by collisionless shock waves. Acceleration occurs according to the first-order Fermi mechanism when particles are scattered by magnetic fluctuations and cross the shock front many times—the so-called diffusive shock acceleration (DSA). The distribution function of the accelerated particles has a significant anisotropy in the upstream of the shock, and this anisotropy leads to development of plasma instabilities in the upstream. The development of plasma instabilities also leads to a magnetic field amplification (MFA) in the upstream. The pressure of the accelerated particles can be on the order of the total momentum flux flowing onto the shock. In this case, the pressure gradient of the accelerated particles leads to modification of the plasma flow in the upstream.

The Monte Carlo code employed to describe the DSA and MFA is one-dimensional, stationary, nonlinear and plane-parallel. We consider acceleration of protons. The nonrelativistic model described in [71] has been modified to be applicable for the case of relativistic shocks. Within this model all the particles are divided into background and accelerated ones: the accelerated particles are those that have crossed the shock front from the downstream to the upstream at least once. The accelerated particles are treated individually. Background particles are described via macroscopic parameters under the assumption of their local Maxwell distribution function up to the injection point near the shock front in the upstream. After this point, the background particles are described as particles, which allows us not to add any additional parameters to describe the injection of particles into the acceleration process. During the propagation of particles, they are scattered elastically and isotropically in the rest frame of the scattering centers according to the pitch-angle scattering approach [72–76]. The reference frame of the scattering centers for background particles moves with a speed *u*(*x*) relative to the rest frame of the shock. The rest frame of the scattering centers for accelerated particles moves with a speed *u*(*x*) + *vscat*(*x*) relative to the rest frame of the shock. *u*(*x*) is the speed of the background plasma flow. The presence of *vscat*(*x*) is due to the fact that with the development of resonant instability, the modes propagating only in a certain direction relative to the background plasma are amplified. *vscat*(*x*) will be determined below.

Between the scatterings the particles uniformly move straightforward. The distance between particle scatterings is proportional to its mean free path, which we define as:

$$\lambda(\mathbf{x}, p) = \frac{1}{\frac{1}{\lambda\_{B, \text{st}}(\mathbf{x}, p)} + \frac{1}{\lambda\_{\text{ss}}(\mathbf{x}, p)}},\tag{5}$$

where *k* is the absolute value of the wavenumber of modes that make up magnetic fluctuations, *p* is the particle momentum, *x* is the particle coordinate that is counted from the front of the shock wave. Negative values of *x* correspond to the upstream, positive values correspond to the downstream.

$$
\lambda\_{B,st}(\mathbf{x}, p) = \frac{pc}{eB\_{ls,st}(\mathbf{x}\_\prime k\_{res})} \,\prime \tag{6}
$$

$$B\_{ls,st}(\mathbf{x},k) = \sqrt{4\pi \int\_0^k \mathcal{W}(\mathbf{x},k)dk},\tag{7}$$

where *W*(*x*, *k*) is the spectral density of the turbulence energy.

$$
\lambda\_{\rm ss}(\mathbf{x}, p) = \left(\frac{pc}{\pi e}\right)^2 \frac{1}{\int\_{k\_{\rm res}}^{\infty} \frac{W(\mathbf{x}, k)}{k} dk} \,\tag{8}
$$

$$\frac{k\_{\rm res}pc}{eB\_{\rm ls}(\chi\_{\prime}k\_{\rm res})} = 1,\tag{9}$$

*Bls* is the large-scale magnetic field:

$$B\_{ls}(\mathbf{x},k) = \sqrt{4\pi \int\_0^k \mathcal{W}(\mathbf{x},k)dk + B\_{0'}^2} \tag{10}$$

where *B*<sup>0</sup> is the constant longitudinal magnetic field. The particle propagates until it gets outside of the model box to the far downstream or crosses the free escape boundary (FEB) located at *x* = *xFEB* in the far upstream.

The iterative scheme of the employed numerical model allows us to keep the conservation laws of momentum and energy fluxes near the shock front. In the stationary relativistic case these are are formulated as follows. The law of conservation of particle flux (automatically kept here) has the form:

$$
\gamma(\mathbf{x})\beta(\mathbf{x})\pi(\mathbf{x}) = F\_{\mathbf{n}0} \tag{11}
$$

where *γ*(*x*) = 1/ 1 − *β*2(*x*) is the Lorentz factor of the background plasma flow, *β*(*x*) = *u*(*x*)/*c*, *n*(*x*) is the background plasma number density, *Fn*<sup>0</sup> is the particle flux in the far upstream.

The conservation law of momentum flux takes the form:

$$\gamma^2(\mathbf{x})\beta^2(\mathbf{x})\left[m\_P c^2 n(\mathbf{x}) + \Phi\_{\rm l\dot{m}}(\mathbf{x}) + \Phi\_{\rm w}(\mathbf{x})\right] + P\_{\rm l\dot{n}}(\mathbf{x}) + P\_{\rm w}(\mathbf{x}) + F\_{\rm px}^{cr}(\mathbf{x}) = F\_{\rm px0} + Q\_{\rm px}^{\rm ESC} \tag{12}$$

where *Pth*(*x*) is the background plasma pressure, Φ*th*(*x*) = Γ*th*(*x*)*Pth*(*x*)/(Γ*th*(*x*) − 1), Γ*th* is the background plasma adiabatic index, *Pw*(*x*) is the turbulence pressure, *Pw*(*x*) = 0.5 (*k*) *W*(*x*, *k*)*dk*, Φ*w*(*x*) = 3*Pw*(*x*), *Fpx*<sup>0</sup> is the momentum flux in the far unperturbed upstream, *QESC px* is the momentum flux carried away by the accelerated particles through the FEB (*QESC px* = *Fcr px*(*xFEB*)), *Fcr px*(*x*) is the momentum flux of the accelerated particles.

The conservation law of energy flux takes the form:

$$\gamma^2(\mathbf{x})\beta(\mathbf{x})\left[m\_p c^2 n(\mathbf{x}) + \Phi\_{th}(\mathbf{x}) + \Phi\_w(\mathbf{x})\right] + F\_{en}^{cr}(\mathbf{x}) = F\_{en0} + Q\_{en}^{ESC} \tag{13}$$

where *Fen*<sup>0</sup> is the energy flux in the far unperturbed upstream, *QESC en* is the energy flux carried away by accelerated particles through the FEB (*QESC en* = *Fcr en*(*xFEB*)), *Fcr en*(*x*) is the energy flux of the accelerated particles.

To keep the conservation laws (12) and (13), it is necessary to determine the profile of the background plasma flow *u*(*x*) in the upstream and the full compression by the shock *Rtot* by means of an iterative process. At the initial iteration, approximate profiles *u*(*x*), *vscat*(*x*) in the upstream, *W*(*x*, *k*) and the full shock compression *Rtot* are set. Then the particles are propagated and their distribution function is calculated accordingly. In the far downstream, the momentum distribution function of all the particles in the rest frame of the flow is isotropic. Thus, it is possible to determine the adiabatic index of the entire plasma in the downstream Γ*p*<sup>2</sup> based on the obtained particle distribution function. Here and below, the subscript 0(2) denotes the values in the far upstream (downstream). In these designations *Rtot* = *β*0/*β*2.

To find the full compression, we write down the flux conservation laws (11), (12) and (13) for the far upstream and for the downstream.

$$
\gamma\_2 \beta\_2 \eta\_2 = \gamma\_0 \beta\_0 \eta\_0 \tag{14}
$$

$$\begin{aligned} \gamma\_2^2 \beta\_2^2 \left[ m\_p c^2 n\_2 + \Phi\_{p2} + \Phi\_{w2} \right] + P\_{p2} + P\_{w2} &= \gamma\_0^2 \beta\_0^2 \left[ m\_p c^2 n\_0 + \Phi\_{th0} + \Phi\_{w0} \right] + \\ &+ P\_{th0} + P\_{w0} + Q\_{p\chi}^{\mathrm{ESC}}, \end{aligned} \tag{15}$$

$$
\gamma\_2^2 \beta\_2 \left[ m\_p c^2 n\_2 + \Phi\_{p2} + \Phi\_{w2} \right] = \gamma\_0^2 \beta\_0 \left[ m\_p c^2 n\_0 + \Phi\_{th0} + \Phi\_{w0} \right] + Q\_{\varepsilon n}^{\text{ESC}} \tag{16}
$$

where *Pp*<sup>2</sup> is the pressure of the entire plasma in the downstream, Φ*p*<sup>2</sup> = Γ*p*<sup>2</sup>*Pp*2/ Γ*p*<sup>2</sup> − 1 . We determine the current values *QESC en* , *QESC px* and Γ*p*<sup>2</sup> at this iteration of the quantities from the distribution function obtained after particle propagation. The value *Pw*<sup>2</sup> is also calculated based on the equation given below. Thus, three unknowns *Rtot*, *Pp*<sup>2</sup> and *n*<sup>2</sup> remain in the Equations (14)–(16). Solving these equations, we find a new value of *Rtot*. The value *Rtot* for the next iteration is found by averaging the new value and the old one.

A new profile of the flow velocity in the upstream is determined according to the formula based on (12):

$$\gamma\_{new}(\mathbf{x})\beta\_{new}(\mathbf{x}) = \gamma\_{old}(\mathbf{x})\beta\_{old}(\mathbf{x}) + \frac{F\_{px0} + Q\_{px}^{ESC} - F\_{px}^{cr}(\mathbf{x}) - F\_{px}^{th}(\mathbf{x}) - F\_{px}^{w}(\mathbf{x})}{\gamma\_0 \beta\_0 m\_p c^2 n\_0},\tag{17}$$

where the values obtained after propagation of particles are in the right part of the expression. The background plasma momentum flux is *Fth px*(*x*) = *γ*2(*x*)*β*2(*x*)Φ*th*(*x*) + *Pth*(*x*). The turbulence momentum flux is *F<sup>w</sup> px*(*x*) = *γ*2(*x*)*β*2(*x*)Φ*w*(*x*) + *Pw*(*x*). *γold*(*x*)*βold*(*x*) is determined by the flow profile at the previous iteration. Selection of the flow profile based on the expression (17) in the area where the background flow is described in the form of particles, works well in the case of non-relativistic motion, when calculating in this area *Fth px*(*x*) based on the distribution function of background particles. In the relativistic case, as shown in [77], the momentum flow in the iterative process converges to a greater value than in the far upstream, when using (17) near the shock wave front. Hence, following [77] we smooth out the flow profile *u*(*x*) near the shock front. The new speed profile is then averaged with the old one. Below we describe the equations that are used to calculate the values included in *Fth px*(*x*) and *F<sup>w</sup> px*(*x*).

The turbulence energy spectrum defining *F<sup>w</sup> px*(*x*) is found based on the solution of the following equation:

$$\begin{split} \gamma(\mathbf{x})u(\mathbf{x})\frac{\partial \mathcal{W}(\mathbf{x},k)}{\partial \mathbf{x}} + \frac{3}{2} \frac{\partial (\gamma(\mathbf{x})u(\mathbf{x}))}{\partial \mathbf{x}} \mathcal{W}(\mathbf{x},k) + \frac{\partial \Pi(\mathbf{x},k)}{\partial k} &= \\ &= G(\mathbf{x},k)\mathcal{W}(\mathbf{x},k) - \mathcal{L}(\mathbf{x},k), \end{split} \tag{18}$$

where L(*x*, *k*) is the turbulent energy dissipation. The spectral flux of the turbulent energy (turbulent cascade) is:

$$\Pi(\mathbf{x},k) = -\frac{\mathbb{C}^\*}{\sqrt{\rho(\mathbf{x})}} k^{\frac{11}{2}} W(\mathbf{x},k)^{\frac{1}{2}} \frac{\partial}{\partial k} \left(\frac{W(\mathbf{x},k)}{k^2}\right) \tag{19}$$

where *ρ*(*x*) is the background plasma density,

$$\mathcal{C}^\* = \frac{\mathfrak{Z}}{11} \mathcal{C}^{-\frac{3}{2}}\_{\text{Kolm}^\prime} \tag{20}$$

where *C*Kolm is the Kolmogorov's constant, which is here taken equal to *C*Kolm = 1.6. The expression (19) is derived from [78]. *G*(*x*, *k*) is the growth rate of the turbulent energy in the background plasma rest frame due to plasma instabilities. Similar to [71], here we considered the growth rate of current instabilities—the Bell's non-resonant [79] and resonant instability. The accelerated particle current, which determines the growth rates, is calculated in the rest frame of the scattering centers, after propagation of particles.

Differentiating the Equations (12) and (13) by *x*, using (11), and excluding from the equations the term proportional to *n*(*x*) we get the following relation:

$$\begin{split} \beta(\mathbf{x}) \Big( \frac{d\Phi\_{th}(\mathbf{x})}{d\mathbf{x}} + \frac{d\Phi\_{\mathbf{w}}(\mathbf{x})}{d\mathbf{x}} \Big) + \frac{dF\_{en}^{cr}(\mathbf{x})}{d\mathbf{x}} + \frac{1}{\gamma(\mathbf{x})} \frac{d(\gamma(\mathbf{x})\beta(\mathbf{x}))}{d\mathbf{x}} (\Phi\_{th}(\mathbf{x}) + \Phi\_{\mathbf{w}}(\mathbf{x})) = \\ = \beta(\mathbf{x}) \Big( \frac{dP\_{th}(\mathbf{x})}{d\mathbf{x}} + \frac{dP\_{\mathbf{w}}(\mathbf{x})}{d\mathbf{x}} + \frac{dF\_{\mathbf{p}\mathbf{x}}^{cr}(\mathbf{x})}{d\mathbf{x}} \Big). \end{split} \tag{21}$$

If we assume that each of the components of the system is affected only by the change of (*γ*(*x*)*β*(*x*)), that is, the change in energy and momentum flux occurs adiabatically due to the change of the flow velocity, we can divide the Equation (21) into separate adiabatic equations for the components:

$$
\beta(\mathbf{x}) \frac{d\Phi\_{th}(\mathbf{x})}{d\mathbf{x}} + \frac{1}{\gamma(\mathbf{x})} \frac{d(\gamma(\mathbf{x})\beta(\mathbf{x}))}{d\mathbf{x}} \Phi\_{th}(\mathbf{x}) = \beta(\mathbf{x}) \frac{d P\_{th}(\mathbf{x})}{d\mathbf{x}},\tag{22}
$$

$$
\beta(\mathbf{x}) \frac{d\Phi\_w(\mathbf{x})}{d\mathbf{x}} + \frac{1}{\gamma(\mathbf{x})} \frac{d(\gamma(\mathbf{x})\beta(\mathbf{x}))}{d\mathbf{x}} \Phi\_w(\mathbf{x}) = \beta(\mathbf{x}) \frac{dP\_w(\mathbf{x})}{d\mathbf{x}},\tag{23}
$$

$$\frac{dF\_{\rm en}^{cr}(\mathbf{x})}{d\mathbf{x}} = \beta(\mathbf{x}) \frac{dF\_{\rm px}^{cr}(\mathbf{x})}{d\mathbf{x}}.\tag{24}$$

The Equation (23) in this case is equivalent to the Equation (18) integrated by *k* in the absence of MFA and dissipation. In the presence of MFA and dissipation after integration by *k*, the Equation (18) will take the form (similar to (23)):

$$\begin{split} \beta(\mathbf{x}) \frac{d\Phi\_{\mathbf{w}}(\mathbf{x})}{d\mathbf{x}} + \frac{1}{\gamma(\mathbf{x})} \frac{d(\gamma(\mathbf{x})\beta(\mathbf{x}))}{d\mathbf{x}} \Phi\_{\mathbf{w}}(\mathbf{x}) &= \beta(\mathbf{x}) \frac{dP\_{\mathbf{w}}(\mathbf{x})}{d\mathbf{x}} + \\ &+ \frac{1}{c\gamma(\mathbf{x})} \int\_{(k)} G(\mathbf{x}, k) \mathcal{W}(\mathbf{x}, k) d\mathbf{k} - \frac{1}{c\gamma(\mathbf{x})} L(\mathbf{x}), \end{split} \tag{25}$$

where the dissipation term is *L*(*x*) = (*k*) L(*x*, *k*)*dk*. Here we use the following expression for the turbulent energy dissipation:

$$\mathcal{L}(\mathbf{x},k) = v\_{\Gamma}(\mathbf{x}) \frac{k^2}{k\_{th}} \mathcal{W}(\mathbf{x},k)\_{\prime} \tag{26}$$

$$v\_{\Gamma}(\mathbf{x}) = \frac{B\_{ls}(\mathbf{x}, k\_{th})}{\sqrt{4\pi\rho(\mathbf{x})}},\tag{27}$$

$$\frac{k\_{th}c\sqrt{k\_bT\_{th}(\mathbf{x})}}{eB\_{ls}(\mathbf{x}\_sk\_{th})} = 1,\tag{28}$$

where *kb* is the Boltzmann's constant, *Tth*(*x*) is the background plasma temperature (*Pth*(*x*) = *n*(*x*)*kbTth*(*x*)). In this model, it is assumed that in the downstream Π(*x*, *k*) = 0 and L(*x*, *k*) = 0.

From the comparison of Equations (23) and (25) one can see that there are two additional terms in the right side of (25). To fulfill the total energy conservation Equation (21), these terms must be compensated by introducing additional terms into the equations for the remaining components of the system. We assume that the dissipation of the turbulent energy flow leads to heating of the background plasma. An increase in the energy flow turbulence due to MFA can be compensated by scattering accelerated particles in the frame

of scattering centers moving with the speed *u*(*x*) + *vscat*(*x*). Then the equation for the energy flux of the accelerated particles will take the form:

$$\mathcal{L}\frac{dF\_{\text{eff}}^{\text{cr}}(\mathbf{x})}{d\mathbf{x}} = [u(\mathbf{x}) + v\_{\text{scat}}(\mathbf{x})] \frac{dF\_{\text{px}}^{\text{cr}}(\mathbf{x})}{d\mathbf{x}}.\tag{29}$$

Within the considered geometry we introduce *vscat*(*x*) as follows:

$$v\_{scat}(\mathbf{x}) = \min\left(v\_{ampl}(\mathbf{x}), v\_{A,eff}(\mathbf{x})\right),\tag{30}$$

$$v\_{ampl}(\mathbf{x}) = -\frac{\int\_{(k)} G(\mathbf{x}, k) W(\mathbf{x}, k) dk}{\gamma(\mathbf{x}) \frac{d F\_{px}^{cr}(\mathbf{x})}{d \mathbf{x}}} \,\tag{31}$$

$$v\_{A,eff}(\mathbf{x}) = -\frac{B\_{eff}}{\sqrt{4\pi\rho'}}\tag{32}$$

$$B\_{eff}(\mathbf{x}) = \sqrt{4\pi \int\_{(k)} W(\mathbf{x}, k) dk + B\_{0\prime}^2} \tag{33}$$

up to the injection point of the background particles in the upstream. After the injection point and before the shock front *vscat*(*x*) = *vampl*(*x*). In the downstream *vscat*(*x*) = 0. In this model, we assume that part of the energy flux taken from the accelerated particle flux in the region where *vampl*(*x*) < *vA*,*eff*(*x*) , goes to heat the background plasma. Thus, the equation determining the change in the energy flux of the background plasma has the form:

$$\begin{split} \beta(\mathbf{x}) \frac{d\Phi\_{th}(\mathbf{x})}{d\mathbf{x}} + \frac{1}{\gamma(\mathbf{x})} \frac{d(\gamma(\mathbf{x})\beta(\mathbf{x}))}{d\mathbf{x}} \Phi\_{th}(\mathbf{x}) = \beta(\mathbf{x}) \frac{d P\_{th}(\mathbf{x})}{d\mathbf{x}} + \frac{v\_{diss}(\mathbf{x})}{c} \frac{d F^{\sigma}\_{\mathbf{px}}(\mathbf{x})}{d\mathbf{x}} + \\ &+ \frac{1}{c\gamma(\mathbf{x})} L(\mathbf{x}), \end{split} \tag{34}$$

where *vdiss*(*x*) = *vA*,*eff*(*x*) − *vampl*(*x*) if *vampl*(*x*) < *vA*,*eff*(*x*) , in the opposite limit *vdiss*(*x*) = 0. After the injection point of background particles and in the downstream *vdiss*(*x*) = 0.

Substituting the expression for Φ*th*(*x*) into the Equation (34), we obtain the following equation for the background plasma pressure:

$$\begin{split} \gamma(\mathbf{x})\beta(\mathbf{x})\frac{d}{d\mathbf{x}}\frac{P\_{th}(\mathbf{x})}{\Gamma\_{th}(\mathbf{x})-1} + \frac{d(\gamma(\mathbf{x})\beta(\mathbf{x}))}{d\mathbf{x}}\frac{\Gamma\_{th}(\mathbf{x})P\_{th}(\mathbf{x})}{\Gamma\_{th}(\mathbf{x})-1} &= \\ &= \gamma(\mathbf{x})\frac{v\_{dis}(\mathbf{x})}{c}\frac{dF\_{px}^{\mathbf{r}}(\mathbf{x})}{d\mathbf{x}} + \frac{1}{c}L(\mathbf{x}).\end{split} \tag{35}$$

The solution of the Equation (35) defines *Fth px*(*x*) in the expression (17).

To solve the Equation (18), one needs to set the profile *W*(*x*, *k*) on the FEB *x* = *xFEB*. We define *<sup>W</sup>*(*xFEB*, *<sup>k</sup>*) <sup>∼</sup> *<sup>k</sup>*<sup>−</sup> <sup>5</sup> <sup>3</sup> —the Kolmogorov's spectrum at the FEB with the energy-carrying scale *Len*. *W*(*xFEB*, *k*) is normalized as follows:

$$\int\_{(k)} \mathcal{W}(\mathfrak{x}\_{FEB}, k) dk = \frac{B\_0^2}{4\pi}. \tag{36}$$

Based on the developed Monte Carlo model, we have performed calculations of proton acceleration by mildly relativistic shocks, which are thought to come out during explosions of some SN kinds. In this case, the shocks often propagate into the wind of the pre-supernova star. The number density of particles (protons) in the stellar wind can be estimated as:

$$m\_w(r) = \frac{\dot{M}}{4\pi v\_w m\_p r^2},\tag{37}$$

where *M*˙ is the mass-loss rate of the pre-supernova, *vw* is the speed of the stellar wind, *r* is the distance from the center of the star. For *<sup>M</sup>*˙ <sup>=</sup> <sup>10</sup>−4*M*yr−1, *vw* = 1000 km s−<sup>1</sup> and *<sup>r</sup>* <sup>=</sup> <sup>10</sup><sup>15</sup> cm: *nw* <sup>≈</sup> <sup>3</sup> <sup>×</sup> <sup>10</sup><sup>6</sup> cm−3. This estimate is used to estimated the number density *<sup>n</sup>*<sup>0</sup> in the far upstream. *Len* <sup>=</sup> <sup>3</sup> <sup>×</sup> 1016 cm in all the calculations. We assume that the free escape boundary is at 0.2 of the current radius of the shock. The calculation parameters and their results are presented in Table 1.

$$
\varepsilon\_B' = \frac{B\_{eff,2}^2}{8\pi\gamma\_0^2 u\_0^2 m\_p n\_0}.\tag{38}
$$

The spatial coordinate *x* in the figures is measured in *rg*<sup>0</sup> = *mpcu*0/*eB*0. For the model A1: *xFEB* <sup>=</sup> <sup>−</sup>4.782 <sup>×</sup> <sup>10</sup><sup>6</sup>*rg*0.

As can be seen from Figure 7, with an increase of the shock speed, the maximum momentum of accelerated particles increases. The value *<sup>B</sup>* also increases with an increase in the velocity of the shock wave in the range of simulated shock velocities 0.1–0.5 *c* and it decreases slightly for the shock velocity of 0.7*c* (see Table 1) where the transition to relativistic shock regime occurs. Particle-in-cell simulations of electron and proton spectra in the trans-relativistic regime were discussed in [61,62].


**Table 1.** The grid of Monte Carlo calculation parameters.

Figure 8 shows particle distribution functions used for evaluations of the models A1, A2 and A3, with varying values of *xFEB* and the magnetic field at the FEB. Note, that with an increase of the FEB distance, the strength of the magnetic field on the FEB decreases proportionally. It can be seen that in these configurations the maximum momenta of the accelerated particles are almost the same. The current of the highest energy particles near the FEB amplify turbulent fluctuations due to the small-scale Bell's instability on scales much smaller than their own gyroradius. In Figure 9, the spectral energy density of the turbulence is shown at various upstream points for the model A1. Furthermore, it takes time to significantly amplify the magnetic field, and thus, the amplitude of the amplified turbulent field in a significant part of the upstream differs slightly from its value at the FEB (see Figure 10). The gyroradius of the highest energy particles near the FEB turns out to be significantly larger than the scale of the amplified fluctuations, which leads to a small contribution of the expression (8) to the free path (5). Thus, the free path (5) of the highest energy particles in a significant part of the upstream is determined by the contribution (6) with the initial turbulent field. That is, the free path (5) of the highest energy particles in a significant part of the upstream is *λ*(*x*, *p*) ∼ *p*/*B*<sup>0</sup> according to the normalization (36). The non-relativistic theory of DSA gives a simple estimate of the maximum momentum *pmax* in the case of an upstream path independent of the coordinate: *λ*(*pmax*)*c*/3*u*<sup>0</sup> ≈ *xFEB*. Accordingly, our model has a good estimate: *pmax* ∼ *B*0*xFEB*. In Figure 11, the independence of the maximum particle momentum from the number density *n*<sup>0</sup> in the far upstream is illustrated. In Figure 12, it is shown, that with an increase of the magnetic field *B*<sup>0</sup> with the other parameters kept constant, the maximum particle momentum increases.

**Figure 7.** Particle distribution function in the shock rest frame, thick curves correspond to the point *x* = 0 (the front shock), thin curves correspond to a point *x* = *xFEB*. The correspondence of the certain model Table 1 is reflected in the legend.

We have made the following estimate of the acceleration time *τ<sup>a</sup>* of the particles to their maximum momentum. In the estimation, we assume that the mean free path of the highest energy particles with momentum *pmax* in most of the upstream can be estimated by the gyroradius in the magnetic field *B*0. Thus, taking into account that the magnetic field in the downstream is much stronger and, accordingly, the diffusion coefficient of particles is much smaller than in the upstream *<sup>τ</sup><sup>a</sup>* <sup>≈</sup> <sup>3</sup>*D*(*pmax*)/*u*<sup>2</sup> 0, where *<sup>D</sup>*(*pmax*) <sup>=</sup> *<sup>λ</sup>*(*pmax*)*c*/3 <sup>≈</sup> *pmaxc*2/3*eB*<sup>0</sup> is the diffusion coefficient of the highest energy particles in the upstream. Accordingly, the estimate of the particle acceleration time to the momentum *pmax* has the form:

$$
\tau\_a \approx 3.5 \times 10^5 \left(\frac{u\_0}{0.1c}\right)^{-2} \left(\frac{B\_0}{3 \times 10^{-3} \text{G}}\right)^{-1} \left(\frac{p\_{\text{max}}}{10^5 m\_p c}\right) \text{s.} \tag{39}
$$

**Figure 8.** The particle distribution function in the shock rest frame, the thick curves correspond to the point *x* = 0 (the front shock), the thin curves correspond to the point *x* = *xFEB*. The correspondence to Table 1 is reflected in the legend.

In the considered models we assume that the current radius of the shock is *Rf* = 5|*xFEB*|. The expansion time of the supernova remnant can be estimated as the ratio of the current shock radius to the shock velocity. Accordingly, the estimate for the expansion time has the form:

$$
\pi\_{\rm exp} \approx 8.3 \cdot 10^5 \left( \frac{|\chi\_{\rm FB}|}{5 \cdot 10^{14} \text{cm}} \right) \left( \frac{u\_0}{0.1c} \right)^{-1} \text{s.} \tag{40}
$$

It can be seen from the expressions (39), (40) that the expansion time is longer than the acceleration time for the model (see Table 1), the maximum momentum can be estimated from Figure 7, which confirms the self-consistency of the stationary model we have employed.

**Figure 9.** Spectral turbulent energy density for the model A1. The green solid curve corresponds to the point *x* = *xFEB*. The blue dashed curve corresponds to the point *x* ≈ 0.75 · *xFEB*. The red dash-point curve corresponds to the point *x* ≈ 0.5 · *xFEB*. The thin purple dashed curve corresponds to the point *x* ≈ 0.1 · *xFEB*. The thin magenta solid curve corresponds to the point *x* = 0 (the shock front).

**Figure 10.** Background plasma velocity profile (upper panel) and magnetic field profile (bottom panel) for the model B1.

**Figure 11.** Particle distribution function in the shock rest frame, the thick curves correspond to the point *x* = 0 (the front shock), the thin curves correspond to the point *x* = *xFEB*. The correspondence to Table 1 is reflected in the legend.

**Figure 12.** Particle distribution function in the shock rest frame, the thick curves correspond to the point *x* = 0 (the front shock), the thin curves correspond to the point *x* = *xFEB*. The correspondence to Table 1 is reflected in the legend.

#### **4. Discussion and Conclusions**

In order to understand the processes of non-thermal particle acceleration at work in fast and energetic transients we have carried out and presented above a microscopic gyro-scale modeling of the collisionless shock structure and non-thermal particle spectra, which show a strong dependence of both electron and proton acceleration efficiency on the shock obliquity.

In this context one should keep in mind that the particle-in-cell simulation assumed a fixed homogeneous magnetic field and cold flow of the incoming particles in the shock upstream at the boundary of the simulation box. On the other hand, fast non-thermal particles can penetrate into the upstream plasma flow. It is clear from the results of the Monte Carlo modeling presented in Section 3 (see, in particular, Figure 10) that the efficient cosmic ray acceleration at the shock may provide a strong magnetic field amplification and modification of the upstream plasma flow. The cosmic ray pressure gradient in the shock upstream decelerates the incoming plasma flow and the cosmic ray driven instabilities may highly amplify the fluctuating magnetic fields well outside the particle-in-cell simulation box which is limited by a few hundred of the proton gyroradii around the plasma sub-shock. Such simulations show the potential importance of the feedback effects which cannot be modeled with the microscopic particle-in-cell simulations so far. Therefore we have used here the particle-in-cell model to simulate the electron injection and acceleration in sub and mildly relativistic shocks of fast energetic transients together with Monte Carlo modeling. While the Monte Carlo technique describes the structure of upstream flow modified by the accelerated particles at scales well above the proton gyro-scales, it cannot be used to simulate the electron injection where particle-in-cell approach is required. Thus, we used the combination of the two techniques to model radio emission.

In our study, of electron spectrum convergence in quasi-parallel shock of 0.3*c* speed we observed some non-monotonic temporal behavior of the electron spectra as it is shown in Figure 5. This could be due to the development of Bell's instability as it was found earlier in [61] and which mediate the quasi-parallel shock structure at proton gyro-scales. The maximum electron energy achievable in our simulation was about 50 *mec*<sup>2</sup> and we extrapolate the spectral slope to larger energy. We extended the spectra from 50 to 500 *mec*<sup>2</sup> using the same spectral slope, it is enough to model radio spectra. This is an assumption in our case. However, the extrapolation seems to be justified by the presence of a powerlaw such as a spectral component with a similar slope right after the thermal peak in Crumley et al. simulation (see the electron spectra presented in Figure 7 [61] in the electron momentum range between 0.1 and 30 *mpc* ). A similar spectral component is also apparent right after the peak of the spectrum in our steady Monte Carlo simulations (see, e.g., our Figure 8).

We applied 2D kinetic PIC simulations for sub-relativistic flows with velocity of 0.3 c to model the spectra of electrons and protons. The dynamical range of the full kinetic PIC simulations is limited but the electron spectrum extrapolated with account for the results of Monte Carlo simulations is enough to model the observed radio emission in the fast transients where such fast outflows were found. In the Monte Carlo model the particle-in-cell simulated domain where the electron injection occurs correspond to the sub-shock structure which is apparent at *x* = 0 in Figure 10. The plasma compression ratio at the sub-shock of the cosmic ray modified shock is about 3 and the spectra of particles accelerated at the sub-shock in the low energy regime (c.f. the proton spectral slope right after the peak in Figure 8) are consistent with the ones obtained within the particle-in-cell simulation dynamical range. The gyroradius of electrons with Lorentz-factor of 1000 in our simulations is about the gyroradius of protons with Lorentz-factor of 10. The spectral index of proton distribution simulated with Monte Carlo model at Lorentz-factor of 10 is consistent with the PIC electron spectrum extrapolated to ∼<sup>500</sup> *mec*2.

This allows to model the non-thermal emission and to understand the origin of the synchrotron radio emission of FBOTs. Additionally, one may assume an explanation of their hard X-ray spectrum using the dependence of the electron power law distribution

index on the shock obliquity. As can be seen in Figure 1 it leads to the presence of relativistic electron populations of comparable intensities within a range of spectral indices along the curved shock surface expected in a wide angle outflow. In this case the synchrotron radio emission could have rather steep spectral indices, while the X-ray component produced by the inverse Compton scattering of the intense optical radiation by the relativistic electrons with the harder spectral index may appear to be flatter.

The simplified macroscopic Monte Carlo kinetic model was used in Section 3 to estimate the maximum energies of the accelerated cosmic ray nuclei. The model accounted for both the nonlinear feedback effects of the flow modification by the cosmic ray pressure gradient and magnetic turbulence amplification by cosmic ray driven instabilities in the upstream of a plain mildly relativistic shock. The model demonstrated a possibility of PeV regime proton acceleration in the shocks driven by mildly relativistic outflows of fast energetic transients on a few weeks timescale.

Transient objects of different types such as gamma-ray bursts and tidal disruption events [80–82] are widely discussed as possible sources of the observed high energy neutrinos. The search of PeV gamma-ray sources have recently revealed a number of potential galactic pevatron sources (see, e.g., [83–86]). The analysis presented in [59] allowed the authors to conclude that the extended galactic supernova remnants are not likely to be pevatrons, while some other types of galactic sources associated with the SN events in the compact clusters of young massive stars [87] and gamma-ray binaries with compact relativistic stars [88,89] can produce PeV photons and neutrinos. While PeV photons undergo strong attenuation and can be detected mostly from the galactic sources, high energy neutrinos from extragalactic energetic transients can be detected with currently operating and future neutrino observatories simultaneously with LSST optical transients detection.

**Author Contributions:** All of the authors have equally contributed into the paper. A.B. constructed the model concept. PiC simulations were performed by V.R. Monte Carlo simulations were performed by S.O. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received funding from Russian Science Foundation grant 21-72-20020.

**Institutional Review Board Statement:** The study did not involve humans or animals.

**Informed Consent Statement:** The study did not involve humans.

**Acknowledgments:** We are honored to dedicate this paper to the memory of professor Yury Nikolaevich Gnedin (1935–2018). His great scientific enthusiasm, competence and the unfailing goodwill for many years were extremely important for us. We are grateful to the referees for constructive comments and suggestions. The authors were supported by the RSF grant 21-72-20020. Some of the modeling was performed at the Joint Supercomputer Center JSCC RAS and some—at the "Tornado" subsystem of the Peter the Great Saint-Petersburg Polytechnic University Supercomputing Center.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

