# Numerical modeling and physical interplay of stochastic turbulent
acceleration for non-thermal emission processes.

Sayan Kundu  Discipline of Astronomy, Astrophysics and Space Engineering
Indian Institute of Technology, Indore
Madhya Pradesh, India - 452020  Bhargav Vaidya  Discipline of Astronomy,
Astrophysics and Space Engineering
Indian Institute of Technology, Indore
Madhya Pradesh, India - 452020  Andrea Mignone  Dipartimento di Fisica
Generale, Universita degli Studi di Torino
Via Pietro Giuria 1, 10125 Torino, Italy<EMAIL_ADDRESS>
###### Abstract

Particle acceleration is an ubiquitous phenomenon in astrophysical and space
plasma. Diffusive shock acceleration (DSA) and stochastic turbulent
acceleration are known to be the possible mechanisms for producing very high
energetic particles, particularly in weakly magnetized regions. An interplay
of different acceleration processes along with various radiation losses is
typically observed in astrophysical sources. While DSA is a systematic
acceleration process that energizes particles in the vicinity of shocks,
stochastic turbulent acceleration (STA) is a random energizing process, where
the interaction between cosmic ray particles and electromagnetic fluctuations
results in particle acceleration. This process is usually interpreted as a
biased random walk in energy space, modelled through a Fokker-Planck equation.
In the present work, we describe a novel Eulerian algorithm, adopted to
incorporate turbulent acceleration in the presence of DSA and radiative
processes like synchrotron and Inverse-Compton emission. The developed
framework extends the hybrid Eulerian-Lagrangian module in a full-fledged
relativistic Magneto-hydrodynamic (RMHD) code PLUTO. From our validation tests
and case studies, we showcase the competing and complementary nature of both
acceleration processes. Axisymmetric simulations of an RMHD jet with this
extended hybrid framework clearly demonstrate that emission due to shocks is
localized while that due to turbulent acceleration originates in the backflow
and is more diffuse, particularly in the high energy X-ray band.

Acceleration of particles – Radiation mechanisms: Non-thermal – Plasma –
Turbulence – Methods: Numerical

##  1 Introduction

From giving a universal power-law trend to the cosmic ray spectrum to
explaining the observed emission features of various astrophysical sources,
particle acceleration process plays a crucial role in shaping our
understanding of the nature of various space and astrophysical phenomena.
Several observations require particles to be accelerated to very high energies
in order to explain the energetics in different astrophysical sources. Due to
high electrical conductivity, astrophysical plasma is incapable of sustaining
a global electric field, making it challenging to energize particles in this
scenario. Particle acceleration processes provide an alternative way to
accelerate particles in the absence of a global electric field. The existing
literature (Blandford, 1994; Kirk et al., 1994; Melrose, 1996) suggests three
main approaches to accelerate charged particles in an astrophysical plasma
environment: shock acceleration (DSA), coherent electric field acceleration,
and stochastic acceleration (STA).

In Fermi (1949), Fermi first gave a proper mechanism for accelerating charged
particles to explain the cosmic ray spectrum and the possible origin of high-
energy cosmic ray particles. The mechanism considers relativistic particles
getting scattered by moving inhomogeneities, mainly various plasma waves (MHD
waves for highly relativistic cosmic ray particles (Parker, 1955; Sturrock,
1966; Kulsrud & Ferrari, 1971)), and gaining energy (accelerate) in a
randomized manner. This process is known as stochastic turbulent acceleration
(STA) process. The randomness in the acceleration makes this process
inefficient to energize particles, as suggested by the emission timescales
observed in various astrophysical sources. Nevertheless, STA is considered to
be an important source of turbulence damping in plasma and because of the
omnipresence of turbulence in various astrophysical sources, STA has been
invoked in order to explain the particle acceleration process in solar flares
(Petrosian, 2012), corona above accretion disk of compact object (Dermer et
al., 1996; Liu et al., 2004; Belmont et al., 2008; Vurm & Poutanen, 2009),
supernova remnant (Bykov & Fleishman, 1992; Kirk et al., 1996; Marcowith &
Casse, 2010; Ferrand & Marcowith, 2010), gamma-ray burst (Schlickeiser &
Dermer, 2000), emission from blazars(see Asano & Hayashida (2018) and
references therein), radio lobes of AGN Jets (O’Sullivan et al., 2009), the
diffuse X-ray emission from AGN jets (Fan et al., 2008) along with fermi
bubbles of galaxies (Mertsch & Petrosian, 2019), galaxy clusters (Brunetti &
Lazarian, 2007; Donnert & Brunetti, 2014). Recently STA has also been
suggested as a candidate for the spectral gradient observed in galaxy clusters
(Rajpurohit et al., 2020).

On the other hand, DSA gives a proper framework where particles can interact
with the magnetic inhomogeneities in a way that could only increase the
particles’ energy (Bell, 1978; Drury, 1983; Blandford & Eichler, 1987; Malkov
& Drury, 2001). Due to it’s efficiency, DSA has been used to describe the
particle acceleration process in various astrophysical systems, for example
interplanetary helio-spheric shocks (Jokipii et al., 2007; Perri & Zimbardo,
2015), shock wave of supernova remnant (Bell, 2014), stellar bow shock
(Rangelov et al., 2019), oblique shock in AGN jets (Meli, A. & Biermann, P.
L., 2013), radio relics of galaxy clusters (Kang et al., 2017; van Weeren et
al., 2017; Zimbardo & Perri, 2017). Though DSA is more efficient compared to
STA mechanism, it is believed to only give rise to localized emission where
STA is thought to produce large scale diffusive emission (Fan et al., 2008).

To study these particle acceleration processes in various astrophysical
systems, a numerical approach is imperative because of the multi-scale nature
of the astrophysical plasma. Numerical study for plasma systems can broadly be
categorized into different classes. Direct computation, mainly known as
Particle in Cell (PIC) method, where Newton-Lorenz force law is solved along
with Maxwell’s equation describing the dynamical evolution of the electric and
magnetic field (Giacalone & Ellison, 2000; Nishikawa et al., 2007; Spitkovsky,
2008; Sironi & Spitkovsky, 2011). This first principle approach has been taken
by various researchers to study the particle acceleration processes (Comisso &
Sironi, 2018; Wong et al., 2019; Marcowith et al., 2020). The next numerical
scheme studies the plasma by solving the Vlasov equation for particle
distribution evolution along with Maxwell’s equations (Palmroth et al., 2018).
This scheme provides the advantage to study various plasma behaviour
distinctively. This approach also enables us to study particle acceleration
processes in different physical settings. Similar to this approach, another
approach is often taken to study particle acceleration process in the quasi-
linear approximation where a Fokker-Plank equation is solved in order to
evolve the cosmic ray spectrum due to interaction with MHD waves (Miniati,
2001; Donnert & Brunetti, 2014; Winner et al., 2019; Vazza et al., 2021).

Another numerical procedure studies the plasma in the fluid regime, also known
as magneto-hydrodynamic (MHD) regime. This numerical procedure assumes plasma
to be sufficiently collisional. That is why this procedure is incapable of
capturing the physics of particle acceleration because collisions would make
them to follow a Maxwellian which is in contrast to the observed power-law
trend for the distribution of the accelerated particles. Though fluid approach
fails to capture the particle acceleration process, it provides the background
for the particles to interact with various MHD waves and accelerate. Recently
some research has been devoted to combine the fluid and the PIC approaches
(Bai et al., 2015) to study the DSA (Mignone et al., 2018). The final
numerical method uses Monte-Carlo technique to study particle acceleration by
shock wave (Achterberg & Krulls, 1992; Baring et al., 1994; Marcowith & Kirk,
1999; Wolff & Tautz, 2015) and turbulence (Giacalone & Jokipii, 1999; Teraki &
Asano, 2019). Among all the numerical techniques available the Particle in
Cell method has an advantage (Ostrowski, 1988; Ellison et al., 1990; Ellison &
Double, 2002; Lemoine & Pelletier, 2003; Baring, 2004; Niemiec & Ostrowski,
2006) over all other techniques because PIC not only can model the particle
acceleration process, it also determine the self-generated magnetic
turbulence, and treat them self-consistently with the cosmic ray particles.
But the disadvantage of the PIC technique is, it is computationally very
expensive (Ellison et al., 2013). And in order to bypass this problem other
numerical techniques are used. Among them the kinetic test particle approach
is one of the most efficient one because it could easily be incorporated with
multi-scale simulations.

As most of the sources of particle acceleration act simultaneously in
different regions of astrophysical sources, it is imperative to develop a
framework that can study such region to understand role of individual
acceleration process. In this work, we use the kinetic test particle approach
to study the competing and complimentary actions of DSA and STA. Other
complimentary approaches have focused on studying the role either of the
acceleration processes individually, for example, Miniati et al. (2001);
Miniati (2003); Donnert & Brunetti (2014) have demonstrated the role of STA in
large scale galaxy clusters.

Recently, the existing Lagrangian particle module developed by Vaidya et al.
(2018) in the PLUTO Code (Mignone et al., 2007) has been applied to AGN jets
at kpc scales to study the impact of instabilities and subsequent shocks on
particle acceleration and non-thermal emission (Borse, Nikhil et al., 2021;
Mukherjee et al., 2021). In the present work, we extend this Lagrangian
framework by incorporating the STA process, to study the effect of both DSA
and STA along with their roles in shaping the emission structure in
astrophysical sources. In this context, a macro-particle is a Lagrangian
entity that moves along with the fluid and collects an ensemble of real
particles (e.g. leptons) that are distributed in 1D momentum space.

The paper is organised as follows; in section 2, we discuss the fundamental
theory and necessary equations to describe the STA process. In section 3, we
propose and describe a numerical algorithm to solve the cosmic ray transport
equation. We validate our algorithm and discuss it’s accuracy in section 4. We
analyze STA process in presence and absence of shocks in section 5 and also
discuss the role of several STA parameters through applications to test
situations. Section 6 discusses our findings and summarizes this work.

##  2 Turbulent Particle Acceleration : Theory

This paper aims to study the effect of MHD turbulence and shocks on cosmic ray
transport and their effect on the spectral signature of various astrophysical
systems. The process of interaction between cosmic ray particles and turbulent
plasma is stochastic in nature. Due to the random nature of the interaction,
the energy of a cosmic ray particle follows a biased random walk, which leads
the particle distribution to follow a diffusion equation (Tverskoǐ, 1967):

$\frac{{\partial}f_{0}}{{\partial}t}=\frac{1}{p^{2}}\frac{{\partial}}{{\partial}p}\left(p^{2}D_{pp}\frac{{\partial}f_{0}}{{\partial}p}\right)=\frac{\partial}{\partial
p}\Bigg{(}D_{pp}\frac{\partial f_{0}}{\partial
p}\Bigg{)}+\frac{2D_{pp}}{p}\frac{\partial f_{0}}{\partial p},$ (1)

where, $f_{0}$ is the particle distribution function that depends on time $t$
and momentum $p$. $D_{pp}$ is the diffusion coefficient in momentum space. The
above equation resembles a Fokker-Planck equation (Blandford & Eichler, 1987).
In a magnetized medium charged cosmic rays are also prone to loose their
energy via various radiative and adiabatic losses. Inclusion of these loss
effects along with the random interactions with turbulent magnetic fields
results in the evolution of the distribution of relativistic cosmic ray
particles as follows (Webb, 1989),

${}\nabla_{\mu}(u^{\mu}f_{0}+q^{\mu})+\frac{1}{p^{2}}\frac{\partial}{\partial
p}\Big{[}-\frac{p^{3}}{3}f_{0}\nabla_{\mu}u^{\mu}+\langle\dot{p}\rangle_{L}f_{0}-\Gamma_{visc}p^{4}\tau\frac{\partial
f_{0}}{\partial p}-p^{2}D_{pp}\frac{\partial f_{0}}{\partial
p}-p(p^{0})^{2}\dot{u}_{\mu}q^{\mu}\Big{]}=0.$ (2)

The various terms of the equation are described below:

  1. 1.

$\nabla_{\mu}(u^{\mu}f_{0}+q^{\mu})$ represents the change in $f_{0}$, due to
the spatial transport. $q^{\mu}$ is the spatial diffusion flux, $u^{\mu}$ is
the bulk four-velocity;

  2. 2.

$\frac{p^{3}}{3}f_{0}\nabla_{\mu}u^{\mu}$ defines the energy loss due to
adiabatic expansion;

  3. 3.

$\langle\dot{p}\rangle_{L}f_{0}$ describes the radiative losses, such as
synchrotron and various Inverse Compton (IC) processes;

  4. 4.

$\Gamma_{\rm visc}p^{4}\tau\frac{\partial f_{0}}{\partial p}$ is the particle
acceleration term due to fluid shear (Rieger & Duffy, 2019);

  5. 5.

$p^{2}D_{pp}\frac{\partial f_{0}}{\partial p}$ represents the Fermi $\rm II$
order particle acceleration or STA process (see Eq. (1));

  6. 6.

$p(p^{0})^{2}\dot{u}_{\mu}q^{\mu}$ originates because of the frame
transformation.

Following Vaidya et al. (2018), we neglect the spatial diffusion flux
$q^{\mu}$ as well as the acceleration due to frame transformation (i.e., terms
1 and 6). Also, acceleration due to shear flow ($\Gamma_{\rm visc}=0)$ is not
considered in the present study. Furthermore, the omission of the spatial
diffusion term is compromised by an inclusion of a momentum independent escape
term in Eq. (2) (Achterberg & Krulls, 1992), so that Eq. (2) takes the form,

$\nabla_{\mu}(u^{\mu}f_{0})+\frac{1}{p^{2}}\frac{\partial}{\partial
p}\Big{[}-\frac{p^{3}}{3}f_{0}\nabla_{\mu}u^{\mu}+\langle\dot{p}\rangle_{L}f_{0}-p^{2}D_{pp}\frac{\partial
f_{0}}{\partial p}\Big{]}=-\frac{f_{0}}{T_{\rm esc}},$ (3)

where $T_{\rm esc}$ is the escape timescale. The above equation is same one
used in Vaidya et al. (2018) to update the spectral distribution of a single
macro-particle with the additional contributions related to Fermi $\rm II$
order acceleration and the escape term.

Note that, for relativistic flows, the convective derivative can be expressed
as,

$u^{\mu}\nabla_{\mu}\equiv\gamma\left[\frac{\partial}{\partial
t}+v^{i}\frac{\partial}{\partial x^{i}}\right]=\frac{d}{d\tau},$ (4)

where $\tau$ is the proper time. Assuming pitch angle isotropy in momentum
space ($p$), the distribution function can be written in terms of the number
density of the relativistic particles as $N(p,\tau)dp=4\pi p^{2}f_{0}dp$ with
$N(p,\tau)$ being the number density of non-thermal particles with momentum
between $p$ and $p+dp$. Accordingly Eq. (3) can be written as,

$\displaystyle\frac{dN}{d\tau}+\frac{\partial}{\partial
p}\Bigg{[}-N\nabla^{\mu}u_{\mu}\frac{p}{3}+\frac{\left<\dot{p}\right>_{l}}{p^{2}}N-D_{pp}\frac{\partial
N}{\partial p}$
$\displaystyle+\frac{2ND_{pp}}{p}\Bigg{]}=-N\nabla^{\mu}u_{\mu}-\frac{N}{T_{\rm
esc}}$ (5)

Transforming the independent variable from momentum ($p$) to Lorentz factor
($\gamma$) following $p\approx\gamma m_{0}c$, with $c$ being the speed of
light in vacuum and $m_{0}$ being the mass of the ultra relativistic cosmic
ray particles, Eq. (2) can be expressed as (see Eq. 11 of Tramacere et al.,
2011):

$\begin{array}[]{lcl}\displaystyle\frac{\partial\chi_{p}}{\partial\tau}+\frac{\partial}{\partial\gamma}\left[(S+D_{A})\chi_{p}\right]&=&\displaystyle\frac{\partial}{\partial\gamma}\left(D\frac{\partial\chi_{p}}{\partial\gamma}\right)-\frac{\chi_{p}}{T_{\rm
esc}}\\\ \vskip 6.0pt plus 2.0pt minus
2.0pt\cr&&\displaystyle+Q(\gamma,\tau)\,,\end{array}$ (6)

where $\chi_{p}=N/n$, with $n$ being the number density of the fluid at the
position of macro-particle, $S$ corresponds to radiative and adiabatic losses
and $D_{A}=2D/\gamma^{2}$ corresponds to the acceleration due to Fermi $\rm
II$ order with $D=D_{pp}/m_{0}^{2}c^{2}$. We also include $Q(\gamma,\tau)$ as
a source term in Eq. (6), which accounts for particle injection process from
external sources.

A numerical approach to solve Eq. (6) without the terms on the right hand side
and $D_{A}$ has been discussed in an earlier work (Vaidya et al., 2018), along
with the particle energization through $1^{\rm st}$-order Fermi acceleration
at shocks. The numerical method for DSA has then recently been improved to
account for the history of particle spectra by Mukherjee et al. (2021) and
will be repeated here for completeness.

The improved version of the DSA routine includes a convolution of the upstream
spectra to the downstream region of the shock in an instantaneous steady state
manner. In particular, as the macro-particle crosses the shock, its downstream
spectra is updated as follows:

$\displaystyle\chi_{p}^{\rm
down}(\gamma)\propto\int_{\gamma_{\min}}^{\gamma}\chi_{p}^{\rm
up}(\gamma^{\prime})G(\gamma,\gamma^{\prime})\frac{d\gamma}{\gamma}$ (7)

where, $\chi_{p}^{\rm up}(\gamma)$ is the distribution function far upstream
and $\chi_{p}^{\rm down}(\gamma)$ is the steady state downstream distribution
function, $G(\gamma,\gamma^{\prime})=(\gamma/\gamma^{\prime})^{-m+2}$, with
$m=3r/(r-1)$ and $r$ is the compression ratio. Here, $\gamma_{\min}$ is the
minimum value of Lorentz factor obtained from the upstream spectrum. The value
of $\gamma_{\max}$, the upper-limit of the convolution, is evaluated by
equating timescales due to radiative losses and various acceleration processes
(i.e., DSA and STA) (Böttcher & Dermer, 2010; Mimica & Aloy, 2012; Vaidya et
al., 2018). Further, it is also ensured that the Larmor radius of the highest
energetic lepton within a macro-particle has a radius equal to or less than
one grid cell width. Further details are explicitly mentioned in (Vaidya et
al., 2018; Mukherjee et al., 2021).

###  2.1 Momentum diffusion coefficient ($D$)

The micro-physical processes of the turbulent interaction are encapsulated in
the transport coefficients of Eq. (6). The mathematical form of these
transport coefficients due to different interactions of cosmic ray and
turbulent magnetized medium have been derived for Alfvènic turbulence (see,
for instance, Schlickeiser, 2002; Brunetti & Lazarian, 2007; O’Sullivan et
al., 2009).

In this work, we will consider STA following a 1D energy spectrum expressed as
a power-law in terms of wave vector norm $|\boldsymbol{k}|=k$ with exponent
$-q$,

$\displaystyle W(k)\sim k^{-q},$ (8)

where, $W(k)$ is the turbulent energy spectrum in Fourier space. The momentum
diffusion coefficient can therefore be expressed as (Schlickeiser, 1989;
O’Sullivan et al., 2009),

$\displaystyle D_{pp}\approx\beta_{A}^{2}\frac{\delta
B^{2}}{B^{2}}\Big{(}\frac{r_{g}}{\lambda_{\max}}\Big{)}^{q-1}\frac{p^{2}c^{2}}{r_{g}c}\propto
p^{q},$ (9)

where $p$ is the momentum of the cosmic ray particles, $D_{pp}$ is the
momentum diffusion coefficient, $\beta_{A}$ is the Alfv́en velocity normalized
to the speed of light, $B$ is the mean magnetic field, $\delta B$ its
fluctuations, $r_{g}$ is the particle gyroradius and $\lambda_{\max}$ is the
maximum correlation length of the turbulent medium.

With the definitions above, the systematic acceleration timescale ($t_{A}$)
for STA can be written as

$\displaystyle t_{A}\approx\beta_{A}^{-2}\frac{l}{c}.$ (10)

where $l$ (the mean free path of the cosmic ray particle) can be expressed as

$\displaystyle l\approx\frac{B^{2}}{\delta
B^{2}}\Big{(}\frac{r_{g}}{\lambda_{\max}}\Big{)}^{1-q}r_{g}.$ (11)

Therefore, the acceleration timescale (Eq. (10)) in terms of $\gamma$ could be
expressed as,

$\displaystyle t_{A}\approx\frac{A^{2}}{2}\rho c(m_{0}\gamma
c^{2})^{2-q}B^{q-4}\lambda_{\max}^{q-1},$ (12)

where, $A=B/\delta B$ defines the turbulence level whose value is set to unity
for the present study (O’Sullivan et al., 2009).

###  2.2 Timescales

The processes described in Eq. (6) involve separate timescales due to
different radiative losses and STA process. These timescales can be expressed
in terms of the particle Lorentz factor $\gamma$ as follows:

  1. 1.

Radiative losses time due to Inverse Compton (IC) in Thompson limit and
synchrotron radiation, $t_{L}\propto 1/\gamma$;

  2. 2.

Diffusion time due to Fermi $\rm II$ order momentum diffusion
$t_{D}\propto\big{(}\frac{\gamma}{\gamma_{s}})^{2-q}$, for the chosen
diffusion coefficient $D\propto\Big{(}\frac{\gamma}{\gamma_{s}}\Big{)}^{q}$.
The value of $t_{D}$ therefore becomes a constant, $t_{D}=1/D_{0}$ with a
choice of $q=2$, where $D_{0}$ is the proportionality constant. Here,
$\gamma_{s}$ defines scale Lorentz factor which we have taken it to be unity
for all the cases considered in this work;

  3. 3.

The acceleration timescale $t_{A}=t_{D}/2$, estimated from Eq. (6) with the
acceleration coefficient $D_{A}=2D/\gamma$.

These considerations are of crucial importance in devising a numerical scheme
for the solution of Eq. (6), since an explicit method would demand $\Delta
t<{\min}\\{t_{L},\,t_{D},\,t_{A}\\}$ for stability reason.

##  3 Turbulent Particle Acceleration : Algorithm

###  3.1 Numerical Method

Eq. (6) is a non-homogeneous, convection-diffusion like partial differential
equation (PDE) with variable coefficients. This equation combines both
hyperbolic and parabolic terms. The non-homogeneous character of the equation
is attributed to the presence of the source and sink terms.

While various numerical methods for the numerical solution of Eq. (6) have
been proposed (see, for instance Chang & Cooper, 1970; Winner et al., 2019),
here we take a more up-to-date and refined approach based on the employment of
Runge-Kutta IMplicit-EXplicit (RK-IMEX) schemes whereby the hyperbolic term of
the PDE are treated using an upwind Godunov-type explicit formalism while the
parabolic (diffusion) term is handled implicitly.

Also, in order to account for the large range of values taken by the particle
Lorentz factor $\gamma$, we employ a logarithmically spaced grid to provide
equal resolution per decade.

To this end, we first introduce a coordinate transformation for the
independent coordinate $\gamma\in[\gamma_{\min},\,\gamma_{\max}]$ in the
following way,

$\xi(\gamma)=\frac{\log(\gamma/\gamma_{\min})}{\log(\gamma_{\max}/\gamma_{\min})},$
(13)

where, $\xi\in[0,1]$ is the transformed (logical) coordinate. Eq. (6) is then
rewritten as,

$\frac{\partial\chi}{\partial\tau}+\xi^{\prime}\frac{\partial}{\partial\xi}(H\chi)=\xi^{\prime}\frac{\partial}{\partial\xi}\Big{[}D\xi^{\prime}\frac{\partial\chi}{\partial\xi}\Big{]}-\frac{\chi}{T_{esc}}+Q\,$
(14)

where we have dropped the subscript $p$ for ease of notation, while
$\xi^{\prime}$ is the Jacobian of this transformation given by Eq. (13),

$\xi^{\prime}=\frac{d\xi}{d\gamma}=\frac{1}{\gamma\log(\gamma_{\max}/\gamma_{\min})}\,,$
(15)

while $H=S+D_{A}$, from Eq. (6).

In order to apply the RK-IMEX scheme, we discretize Eq. (14) on a one-
dimensional mesh of $N$ points using the method of lines,

$\frac{d\chi_{i}}{dt}={\cal A}_{i}+{\cal D}_{i}+{\cal S}_{i},$ (16)

so that the original PDE becomes a system of ordinary differential equations
at the nodal points $i=i_{b},\,...,\,i_{e}$, with $N=i_{e}-i_{b}+1$. In Eq.
(16), ${\cal A}_{i}$ is the advection term, ${\cal D}_{i}$ is the diffusion
term and ${\cal S}_{i}$ accounts for accounts for source and sink terms.

The advection term ${\cal A}_{i}$ is discretized in conservative fashion using
the nonlinear Van Leer flux limiter scheme (Van Leer, 1977),

${\cal A}_{i}=-\xi^{\prime}_{i}\frac{{\cal F}^{\rm adv}_{i+\frac{1}{2}}-{\cal
F}^{\rm adv}_{i-\frac{1}{2}}}{\Delta\xi},$ (17)

where the advection flux follows an upwind selection rule,

${\cal F}_{i+\frac{1}{2}}^{\rm
adv}=\begin{cases}H(\gamma_{i+\frac{1}{2}})\chi_{i+\frac{1}{2}}^{L}&\quad{}\quad
H(\gamma_{i+\frac{1}{2}})>0\\\
H(\gamma_{i+\frac{1}{2}})\chi_{i+\frac{1}{2}}^{R}&\quad{}\quad
H(\gamma_{i+\frac{1}{2}})<0\,.\end{cases}$ (18)

The left and right states $\chi^{L}_{i+\frac{1}{2}}$ and
$\chi^{R}_{i+\frac{1}{2}}$ are constructed up to $2^{\rm nd}$-order accuracy
in space using a slope limiter to prevent oscillations around extrema,

$\begin{array}[]{l}\displaystyle\chi^{L}_{i+\frac{1}{2}}=\chi_{i}+\frac{\delta\chi_{i}}{2},\\\
\vskip 6.0pt plus 2.0pt minus
2.0pt\cr\displaystyle\chi^{R}_{i+\frac{1}{2}}=\chi_{i+1}-\frac{\delta\chi_{i+1}}{2},\end{array}$
(19)

with the $\Delta\chi_{i}$ is the harmonic mean slope limiter (Van Leer, 1977),

$\delta\chi_{i}=\begin{cases}\displaystyle\frac{2\Delta\chi_{i+\frac{1}{2}}\Delta\chi_{i-\frac{1}{2}}}{\Delta\chi_{i+\frac{1}{2}}+\Delta\chi_{i-\frac{1}{2}}}&\quad{\rm
if}\quad\Delta\chi_{i+\frac{1}{2}}\Delta\chi_{i-\frac{1}{2}}>0\\\
0&\quad{}\quad{\rm otherwise}\end{cases}$ (20)

where, $\Delta\chi_{i\pm\frac{1}{2}}=\pm(\chi_{i\pm 1}-\chi_{i})$. Note that
this scheme is $2^{\rm nd}$-order accurate away from discontinuities and that
the reconstruction step demands for 2 ghost zones beyond the active domain
cells.

For the diffusion term ${\cal D}_{i}$, we also adopt a conservative formalim
and choose a central differencing approach yielding $2^{\rm nd}$-order
accuracy in the uniform $\xi$ grid:

${\cal D}_{i}=\xi^{\prime}_{i}\frac{{\cal F}^{\rm diff}_{i+\frac{1}{2}}-{\cal
F}^{\rm diff}_{i-\frac{1}{2}}}{\Delta\xi},$ (21)

where,

${\cal F}_{i+\frac{1}{2}}^{\rm
diff}=\left(\xi^{\prime}D(\gamma,t)\right)_{i+\frac{1}{2}}\Bigg{(}\frac{\chi_{i+1}-\chi_{i}}{\Delta\xi}\Bigg{)},$
(22)

is the diffusion flux constructed following a central difference approach.

In the RK-IMEX approach, the advection is carried out explicitly while the
diffusion operator and the source terms are handled implicitly. This allows to
overcome the restrictive time step limitation $\Delta
t\lesssim\Delta\xi^{2}/(\xi^{\prime}D)$ imposed by a typical explicit
discretization.

We have implemented two similar approaches for the temporal integration of Eq.
(16) in the PLUTO code. The first one is the Strong Stability Preserving (SSP)
scheme (2,2,2) of Pareschi & Russo (2005).

Omitting the subscript $i$ for simplicity,

$\begin{array}[]{lcl}\chi^{(1)}&=&\displaystyle\chi^{(n)}+\Delta t\alpha{\cal
D}^{(1)}\\\ \vskip 6.0pt plus 2.0pt minus
2.0pt\cr\chi^{(2)}&=&\displaystyle\chi^{(n)}+\Delta t\Big{[}{\cal
A}^{(1)}+(1-2\alpha){\cal D}^{(1)}+\alpha{\cal D}^{(2)}\Big{]}\\\ \vskip 6.0pt
plus 2.0pt minus 2.0pt\cr\chi^{(n+1)}&=&\displaystyle\chi^{(n)}+\frac{\Delta
t}{2}\Big{[}{\cal A}^{(1)}+{\cal A}^{(2)}+{\cal D}^{(1)}+{\cal
D}^{(2)}\Big{]},\end{array}$ (23)

where $\Delta t$ is the time-step, $\alpha=1-1/\sqrt{2}$.

For the second approach we choose ARS(2,2,2) scheme due to Ascher et al.
(1997):

$\begin{array}[]{lcl}\chi^{(1)}&=&\displaystyle\chi^{(n)}+\Delta
t\Big{[}\alpha{\cal A}^{(n)}+\alpha{\cal D}^{(1)}\Big{]}\\\ \vskip 6.0pt plus
2.0pt minus 2.0pt\cr\chi^{(n+1)}&=&\displaystyle\chi^{(n)}+\frac{\Delta
t}{2}\Big{[}\delta{\cal A}^{(n)}+(1-\delta){\cal A}^{(1)}\Big{]}\\\ \vskip
6.0pt plus 2.0pt minus 2.0pt\cr&&\displaystyle+\frac{\Delta
t}{2}\Big{[}(1-\alpha){\cal D}^{(1)}+\alpha{\cal
D}^{(n+1)}\Big{]},\end{array}$ (24)

where, $\alpha=1-1/\sqrt{2}$, $\delta=1-\frac{1}{2\alpha}$.

Both time-stepping methods require the inversion of two tri-diagonal matrices
per step, which we perform following the Thomas algorithm (Press et al.,
1992). In the present work, we will only show results from the SSP(2,2,2)
scheme since results obtained with the ARS(2,2,2) are similar. Furthermore,
for the sake of comparison, we have also implemented the standard Chang-Cooper
algorithm (Chang & Cooper, 1970; Park & Petrosian, 1996) for solving the
Fokker-Planck Equation.

####  3.1.1 Boundary conditions

In order for our numerical method to operate correctly, boundary conditions
(b.c.) must be specified in the guard (or ghost) zones for $i=i_{b}-1,i_{b}-2$
and likewise for $i=i_{e}+1,i_{e}+2$. Two common b.c. have been routinely
employed (Marcowith et al., 2020). The first one (zero-particle) is a
Dirichlet b.c. requiring the value of the distribution function $\chi$ to
vanish in the ghost zones. This kind of boundary condition in solving the
cosmic ray transport problem is used, for instance, by Winner et al. (2019).
Another boundary condition is a Neumann-like condition requiring zero-flux
across the boundary interface. This condition has been used, for instance, by
Chang & Cooper (1970) to solve the Fokker-Planck equation. The zero-flux b.c.
conserves the integral of $\int\chi d\gamma$ (the analogous of particle number
conservation). For more discussion on the boundary conditions for cosmic ray
transport see Park & Petrosian (1995). Unless otherwise states, we will employ
the zero-flux b.c. to ensure that without the presence of source and sink
terms in Eq. (6), the total number of particles remain conserved. At the
implementation level, we enforce the zero-flux b.c. separately according to
the implicit/explicit stage level in our RK-IMEX update:

  * •

during the implicit diffusion step we impose zero-gradient b.c.:

$\begin{dcases}\begin{array}[]{ll}\chi_{i}^{\rm diff}=\chi_{i_{b}}^{\rm
diff}&\quad{\rm for}\quad i<i_{b}\\\ \vskip 6.0pt plus 2.0pt minus
2.0pt\cr\chi_{i}^{\rm diff}=\chi_{i_{e}}^{\rm diff}&\quad{\rm for}\quad
i>i_{e}\end{array}\end{dcases}.$ (25)

where $\chi^{\rm diff}$ is the solution array immediately before the implicit
step.

  * •

during the explicit hyperbolic update we impose reflective condition

$\begin{dcases}\begin{array}[]{ll}\chi_{i}^{\rm adv}=-\chi_{2i_{b}-i-1}^{\rm
adv}&\quad{\rm for}\quad i<i_{b}\\\ \vskip 6.0pt plus 2.0pt minus
2.0pt\cr\chi_{i}^{\rm adv}=-\chi_{2i_{e}-i+1}^{\rm adv}&\quad{\rm for}\quad
i>i_{e}\end{array}\end{dcases}$ (26)

together with

${\cal F}_{i_{b}-\frac{1}{2}}^{\rm adv}={\cal F}_{i_{e}+\frac{1}{2}}^{\rm
adv}=0\,.$ (27)

In Eq. (26) $\chi^{\rm adv}$ represents the solution array immediately before
the explicit advection step.

A third b.c. is used to assess the accuracy of our algorithm against a
reference or analytical solution. In this case, the value of $\chi$ in the
ghost zones is set to the corresponding analytical value in those zones,
unless otherwise stated.

##  4 Results : Code Validation Tests

In this section we proceed to assess the accuracy of our newly proposed
algorithm. For accuracy calculation, errors will be computed using the $L_{1}$
norm, defined as (Winner et al., 2019):

$L_{1}(N)=\frac{\displaystyle\sum_{i=1}^{N}\left|\chi^{\rm ref}_{i}-\chi^{\rm
num}_{i}\right|\Delta\gamma_{i}}{\displaystyle\sum_{i=1}^{N}\chi^{\rm
ref}_{i}\Delta\gamma_{i}},$ (28)

where, $N$ is the number of energy bins. To further ensure that the scheme
accuracy is not get dominated by the spatial discretization, the increment in
$N$ is compensated by the decrement in $\Delta t$ such that the ratio
$N/\Delta t$ stays constant (Vaidya et al., 2017). In section 5 all the tests
are performed following the zero-flux boundary prescription. Furthermore all
the simulations in this work are performed using the SSP(2,2,2) scheme with
Courant number 0.4, unless otherwise specified.

###  4.1 Simple Advection

Figure 1: Evolution of the particle distribution function and their
corresponding $L_{1}$ error for the simple advection following $S=\gamma^{2}$
(_Top panel_) and $S=-\gamma^{2}$ (_Bottom panel_) case with IMEX-SSP
algorithm. _Left panel:_ shows the numerical (solid lines) and analytical
(black dotted lines) solutions at different times. _Right panel:_ L1 norm
errors at different resolutions (blue dots) and $2^{\rm nd}$-order reference
slope (dashed lines).

We start by considering a simple advection benchmark by setting
$S=k\gamma^{2}$, $D_{A}=D=0$ in Eq. (6). Here we consider two cases, owing to
two diffrent values of $k=\pm 1$. The analytical solution for the case of
$k=-1$ is given by (Kardashev, 1962; Sarazin, 1999):

$\chi_{p}=\begin{dcases}N_{0}\gamma^{-s}(1-\gamma/\gamma_{\rm
cut})^{s-2},&\,\gamma\geq\gamma_{\rm cut}\\\ 0,&\,\gamma\leq\gamma_{\rm
cut}\end{dcases}$ (29)

where, $\gamma_{\rm cut}=1/\tau$, while for $k=1$ we do not encounter such
discontinuity in the result,

$\chi_{p}=N_{0}\gamma^{-s}(1+\gamma/\gamma_{\rm cut})^{s-2}.$ (30)

The initial condition consists of a power-law spectrum,
$\chi_{p}(\gamma,0)=N_{0}\gamma^{-s}$ with $s=3.3$. For the numerical
calculations, we consider the range of $\gamma\in[10,10^{3}]$ as our
computational domain. We show the evolution of $\chi_{p}$ and the
corresponding error for both values of $k$ in Fig. 1, using $128$ bins and
fixed time step $\Delta\tau=0.00375$. The top left panel of Fig. 1 shows the
evolution of $\chi_{p}$ for $k=1$, while the bottom left panel depicts the
same for $k=-1$. The solid curves represent the numerical solutions while the
black dotted curves depict the analytical solution at the corresponding time.
For $k=1$, the distribution function follows the analytical results closely,
while, for $k=-1$ some deviations are observed at a later stage ($\tau=0.03$)
between the analytic and numerical solution, owing to the steepening of the
solution (Eq. 29). A convergence test is shown for both cases in the right
panel of Fig. 1 where we plot the $L_{1}$ error as a function of the number of
bins. Blue dots and the black dashed curve represent, respectively, the
computed $L_{1}$ error and a reference for the $1/N^{2}$ slope. For $k=1$ (top
right) results converge with $2^{\rm nd}$-order accuracy for all resolutions,
while for $k=-1$ (bottom right) a slight deviation from the $2^{\rm nd}$-order
convergence can be observed. This discrepancy is attributed to the
discontinuous nature of analytic solution presented in Eq. (29).

###  4.2 Simple Diffusion

Figure 2: _Left:_ Simple diffusion case for different times where solid lines
show the numerically computed particle distribution function and black dotted
curve depicts analytical solutions. _Right:_ $L_{1}$ error convergence plot
for the Simple diffusion case with IMEX-SSP algorithm.

Next, we solve Eq. (6) in the case of simple diffusion where, $S=D_{A}=0$ and
$D=\gamma^{2}$. The analytical solution for this case can be written as (Park
& Petrosian, 1995),

$\chi_{p}=\frac{1}{\gamma\sqrt{4\pi\tau}}\exp\Bigg{\\{}-\frac{[\log(\gamma_{0}/\gamma)+\tau]^{2}}{4\tau}\Bigg{\\}}$
(31)

We define the computational domain as $\gamma\in[1,10^{6}]$ and employ $128$
logarithmically spaced bins with a fixed time-step $\Delta\tau=0.0375$. The
initial condition is given by the analytical solution (Eq. 31) at $\tau=1.0$
and $\gamma_{0}=100.0$. The results are shown in Fig. 2. The left panel shows
the evolution of the distribution function at different times with solid
(black dotted) curve representing the numerical (analytical) solution. In the
right panel of Fig. 2 the corresponding $L_{1}$ error is shown by varying the
grid size from $32$ to $4096$ bins. Here $2^{\rm nd}$-order convergence is
observed uniformly at all resolutions.

###  4.3 Hard-sphere Equations

The next numerical benchmark is intended to verify the correctness of our
implementation when source and sink terms are present in the Fokker-Planck
equation. Additionally, we also compare our code with the standard Chang-
Cooper algorithm (Chang & Cooper, 1970). For this purpose, we solve the
following Fokker-Planck equation

Figure 3: _Left:_ Evolution of the particle distribution following Eq. (32)
with $\theta=1$. Dashed curves plot results obtained with the Chang-Cooper
scheme, red curves correspond to the SSP(2,2,2) scheme. Different shades
correspond to different times. Black dotted curve depicts the analytical
solutions at the corresponding times. _Right:_ $L_{1}$-norm error convergence
for both Chang-Cooper (blue dots) and SSP(2,2,2) (red dots) schemes. Black
curves shows the reference slopes for the corresponding schemes. Figure 4:
Time evolution of the integral $\int\chi_{p}(\gamma,\tau)d\gamma$ is shown for
the proposed boundary condition (zero flux boundary) along with the boundary
condition where the value of the distribution functions in the ghost zones are
computed from the analytic expression (analytic boundary).

$\frac{\partial\chi_{p}}{\partial\tau}=\frac{\partial}{\partial\gamma}\Big{(}\gamma^{2}\frac{\partial\chi_{p}}{\partial\gamma}-\gamma\chi_{p}(\gamma,\tau)\Big{)}-\theta\chi_{p}\,.$
(32)

The analytical solution of the previous equation can be written as (Park &
Petrosian, 1995),

$\chi_{p}=\frac{e^{-\theta\tau}}{\gamma\sqrt{4\pi\tau}}\exp\Bigg{\\{}-\frac{[\log(\gamma_{0}/\gamma)+2\tau]^{2}}{4\tau}\Bigg{\\}}\,.$
(33)

For the present purpose, we take the inverse escape timescale $\theta=1$ and
the initial particle distribution is obtained by setting $\tau=1.0$,
$\gamma=\gamma_{0}=100.0$ in Eq. (33). The computational domain is taken as
$\gamma\in[1,10^{6}]$ using $128$ (log-spaced) energy bins and a fixed time
step $\Delta\tau=0.0375$.

Numerical solutions obtained via the Chang-Cooper algorithm (dashed curves)
and the SSP(2,2,2) algorithm (solid lines) are shown in the left panel of Fig.
3 at different time (colors). The analytical solution (dotted lines) is also
superposed. The corresponding resolution study is reported in the right panel
of the same figure using $L{{}_{1}}$ error. From the plots it clearly appears
that the Chang-Cooper algorithm converges at $1^{\rm st}$-order rate while the
SSP(2,2,2) scheme gives full $2^{\rm nd}$-order convergence, so that even at
low resolutions the latter yields an error which is already one order of
magnitude smaller than the former. At the resolution of $N=4096$ the SSP
method outperforms the Chang-Cooper scheme by more than $3$ orders of
magnitude.

Notice that, although we employ a conservative discretization, particle number
is not strictly conserved for this test, owing to the chosen boundary
condition which allows a non-zero net flux through the endpoints of the
computational domain. In order to check particle conservation, we have
therefore repeated the same test in absence of sink ($\theta=0$) and by
prescribing the zero-flux b.c. (see section 3.1.1). Results for the previous
and current b.c. are shown in Fig. 4. It can be observed from the figure that
while the integral due to the previous b.c (depicted by green dots),
decreasing with time, the integral due to the zero-flux b.c. (depicted by
black dots) remains constant. This validates the particle number conserving
nature of the proposed boundary condition.

###  4.4 Log-Parabolic Nature of Particle Spectra

Figure 5: _Top left_ : evolution of the particle distribution function with
turbulent acceleration and synchrotron losses with two magnetic field values.
_Top right_ : evolution of the curvature of the distribution function fitted
with a log-normal density profile (Eq. 36). Analytic solution is shown in
solid orange line. _Bottom panel_ : $\chi_{p}(\gamma,\tau)/\gamma^{2}$ as a
function of $\gamma$ at steady state ($\tau=30\,t_{s}$), in agreement with Eq.
A2. The plot shows the increase as $\gamma^{2}$ (black dashed lines) followed
by an exponential cut-off.

It has been shown (Massaro, E. et al., 2004; Massaro et al., 2006) that the
hump structure in the spectral energy distribution (SED) of blazars could be
described with a log-parabolic curve and this log-parabolicity is speculated
to have originated from STA (Tramacere et al., 2011). Here we validate the
log-parabolic nature of the particle distribution due to STA which
consequently translates to log-parabolic nature of observed SED. In
particular, we numerically solve the transport equation (6), in its
conservative form (without source and sink terms) using the zero-flux boundary
prescription, for STA including synchrotron losses. We choose our grid as
$1.0\leq\gamma\leq 10^{9}$ with $5000$ computational bins and
$\Delta\tau=0.003$ with the following transport coefficients,

$S=-C_{0}\gamma^{2}B^{2}\,,\quad D=D_{0}\gamma^{2}\,,\quad
D_{A}=\frac{2D}{\gamma}\,,$ (34)

where $C_{0}=1.28\times 10^{-9}$, $D_{0}=10^{-4}$ sec-1 is the diffusion
constant. We employ $1/D_{0}$ as our unit time ($t_{s}$).

Here, we consider the one-zone model for the blazar emission (Tramacere et
al., 2011) where the geometry of the acceleration region is taken as spherical
with radius $R=5\times 10^{13}$ cm threaded by a magnetic field $B_{\rm mag}$.
In this region, the acceleration is accompanied by the radiative losses.
Moreover, in order to solve Eq. (6) we consider a mono-energetic initial
distribution $\chi_{p}$ corresponding to a total power $L_{\rm inj}=10^{39}$
erg/sec, where

$L_{\rm inj}=N_{\rm part}\frac{4}{3}\pi R^{3}\int\gamma
m_{e}c^{2}\delta(\gamma-\gamma_{\rm inj})d\gamma,$ (35)

where, $N_{part}$ is the total number of particles injected per unit volume
and $\gamma_{\rm inj}=10.0$. The Dirac delta is approximated with a Gaussian
distribution with $\sigma=0.5$ and $\mu=10$ and it is shown by the purple
solid line in left panel of Fig 5. Furthermore, Eq. (6) is solved by adopting
two different magnetic field values $B_{\rm mag}=1$G, $0.1$G and the
corresponding distribution of $\chi_{p}$ for time $\tau=30\,t_{s}$ is shown in
the top left panel of Fig. 5.

The numerical solution is shown in the top left panel of Fig. 5 for different
magnetic field strengths. We point out that the steady-state distribution is
expected to have an ultra-relativistic Maxwellian form as described in Eq.
(A2) in Appendix A. This is confirmed in the bottom panel of Fig. 5 where we
plot $\chi_{p}/\gamma^{2}$ as a function of $\gamma$, showing that our results
correctly reproduce the $\gamma^{2}$-dependence of the spectrum.

Also, in order to quantify the effects of acceleration and radiative losses on
the spectral evolution, we estimate the curvature of the distribution
function. The curvature is measured by finding the peak value of the
distribution function at each time-step which is also the point at which
$t_{L}=t_{A}$ (Katarzyński et al., 2006, see also Sec. 2.2) and subsequently
fitting a log-normal curve through $10$ points centered around $\gamma_{c}$
(the energy at which the maximum occurs). The curvature is then taken as the
inverse of the variance of the best fit. In particular, we adopt the fitting
curve (Kardashev, 1962) as follows:

${}\chi_{\rm
fit}=\frac{A}{\gamma\sigma}\exp{\Big{\\{}-\frac{(\log(\gamma)-\mu-\sigma^{2})^{2}}{4\sigma^{2}}\Big{\\}}}\,,$
(36)

with curvature parameter defined as $r=1/(4\sigma^{2})$. The fitting curve is
a solution to the Fermi $\rm II$ order transport equation (Eq. 6 with
$S=0,D=\gamma^{2}$ and $D_{A}=2D/\gamma$ without sources and sinks) when
$\sigma^{2}=\tau$, therefore the evolution of the curvature $r$ goes as $\sim
1/(4\tau)$. In the top right panel of Fig. 5 we compare $r$ in the
acceleration region (yellow solid line) with $r$ numerically calculated by
fitting Eq. (36) with the particle distribution, at each time, for different
$B_{\rm mag}$ values (red and black dotted lines) .

Our results show that the fitted curvature initially decays with time as
$r\propto t_{s}/4\tau$, following a trend of curvature in the acceleration
region, and then a sudden jump of the curvature to the steady value of
$r=0.25$ can be observed. The results therefore confirm that, during the
earlier stages, STA dominates the evolution of the particle distribution
function and, later, that steady state is reached much faster for stronger
magnetic fields, as confirmed by the curvature evolution (black dots in the
top right plot of Fig. 5).

Summarizing, the numerical benchmarks proposed in this section validate our
implementation and demonstrate that the proposed SSP(2,2,2) scheme is fully
conservative and it provides full $2^{\rm nd}$-order accuracy, in contrast to
its predecessors (i.e. Chang & Cooper, 1970; Winner et al., 2019) with typical
$1^{\rm st}$-order accuracy.

##  5 Effect of Turbulent acceleration in presence of Shocks

In this section, we describe the effect of STA on particle spectra in presence
of shock. In particular, we consider several test situations where the
equations of classical or relativistic MHD are solved using the PLUTO code
(Mignone et al., 2007) along with Lagrangian particles to model the non-
thermal emission (Vaidya et al., 2018; Mukherjee et al., 2021) in presence of
DSA and radiative losses. To study the effects of STA, the newly developed
algorithm (see section 3) has been incorporated into the Lagrangian framework.
The effects of DSA and STA on particle spectra and subsequent non-thermal
emission signatures are compared for various test situations and discussed in
the following.

###  5.1 Non-relativistic MHD Planar shock

Here we perform a simulation of a non-relativistic MHD planar shock
interacting with a single macro-particle in a turbulent medium. We solve the
2D ideal MHD equations with adiabatic equation of state on a Cartesian grid
$x\in[0,40]$ and $y\in[0,2]$ using $1024\times 128$ grid zones. Initially, we
place a shock wave at $x=1$ which moves towards the increasing $x$ direction.
The upstream density and pressure, $\rho_{u}$ and $P_{u}$, are taken as $1$
and $10^{-4}$, respectively, in dimensionless units. A random density
perturbation is added to simulate a non-homogeneous upstream medium. The
magnetic field is defined as
$\boldsymbol{B}=B_{0}(\cos{\theta},\sin{\theta})$, where $\theta$ (the
obliquity) is the angle between $\boldsymbol{B}$ and the direction of shock
normal. For our purpose, we have considered $\theta=30^{\circ}$ while $B_{0}$
is computed from the plasma beta, $\beta=10^{2}=2P_{u}/B_{0}^{2}$.

The physical units adopted for this test are: length $\hat{L}_{0}=100\,{\rm
pc}$, density $\hat{\rho}_{0}=10^{-2}\,{\rm amu}$ while the unit velocity is
taken to be the speed of light $c$. With this choice, pressure will be given
in units of $\hat{P}_{0}=1.5\times 10^{-5}\,{\rm dyne/cm}^{2}$, magnetic field
in units of $\hat{B}_{0}=1.4\times 10^{-2}\,{\rm G}$ and time in units of
$\hat{\tau}_{0}=326.4\,{\rm yrs}$.

The particle is initially located at $(x,y)\equiv(1.5,1.0)$ with an energy
distribution following a steep decreasing power-law profile with index $9$.
The grid ranges in $10\leq\gamma\leq 10^{10}$ using $128$ (log-spaced) bins.
The particle spectrum (Eq. 6) is evolved accounting for synchrotron, inverse-
Compton and adiabatic losses along with the diffusion effect, modelled
following the STA timescale (Eq. 12). Additionally, the effect of shock is
captured via the steady state update convolution, Eq. (7). We also vary the
index $q$ for various turbulent spectra $W(k)\propto k^{-q}$ in three
different scenarios: a) with only STA and no shock, b) both shock and STA and
c) both shock and STA with the latter active only in the downstream region.
The value of $\lambda_{\max}$ is taken to be $\hat{L}_{0}/10^{5}$ for all the
simulations.

Figure 6: _Top section_ : Density map of a fluid with a lagrangian particle
(shown in white dot). The upstream region is shown in blue, and the downstream
region is shown in green. _Bottom section_ : Particle spectra in various
scenarios with $q=2$ turbulence spectrum. Particle spectra _Middle left_ : For
the case of only DSA with a compression ratio of $3.89$ and various losses.
_Middle right_ : In a turbulent medium with various losses but no shock.
_Bottom left_ : With the both shock of same compression ratio, turbulence and
various losses. _Bottom right_ : For turbulence present only at the downstream
region. The black dashed curve shows the particle energy spectrum for the time
when the density map snapshot is taken. Figure 7: Steady-state particle
distribution with shock and turbulence acceleration for various turbulence
spectra. _Left_ : For $q=2$, _Middle_ : for $q=5/3$ and _Right_ : for $q=3/2$.
The solid blue line depicts the case of turbulent acceleration without shock;
the orange line describes the case of shock and turbulence acceleration
considering both regions ahead and behind of shock are turbulent, and the
green line also describes the shock and turbulence acceleration scenario where
only the post-shock region is turbulent. Figure 8: Dependence of $\gamma_{c}$
on various parameters for turbulent acceleration. _Left_ : Dependence of
$\gamma_{c}$ on various $B$ field, _Middle_ : Dependence of $\gamma_{c}$ on
various $\rho$ values and _Right_ : Dependence of $\gamma_{c}$ on various
values of $\lambda_{\max}$. Data point from corresponding simulations are
shown as dots and the result from analytic calculations (see Eq. (37)) is
shown with a dashed line for reference.

The result in the case of weak turbulence ($q=2$) is shown in Fig. 6 where
$t_{A}$ (see section 2.2) is independent of $\gamma$. The top panel shows the
Lagrangian particle position on top of the background gas density distribution
at $t=56.13$. The evolution of the particle energy spectra with various
radiative losses and different acceleration scenarios are shown in the bottom
four panels using different colors (as indicated by the colorbar). The upper
plot depicts the evolution of the particle spectra for the situation when only
DSA is effective. As the shock hits the particle, the spectra becomes flatter
and radiative and adiabatic losses give rise to a cut-off that gradually
shifts from larger values of $\gamma$ to lower values.

The evolution of the particle spectra due to STA alone is shown in the
corresponding right panel. The spectra is now considerably different when
compared to the previous case since, owing to turbulence and losses, particle
energization occurs continuously rather than just when crossing the shock. The
spectra evolves towards the typical steady state of the ultra-relativistic
Maxwellian, as observed in §4.4, with a peak value $\gamma_{c}\sim 10^{8}$
when $t_{A}=t_{L}$. We also notice that the high energy cut-off does not ever
decreases to lower values of $\gamma$ (as for the pure DSA) but, rather, it
settles into a steady state as the result of mutual compensation between
losses and STA.

In the bottom left plot, we show the evolution of the energy spectrum in the
presence of both shock and STA. Both the upstream and the downstream are
turbulent. In this scenario, the distribution function becomes harder than the
initial one owing to the presence of upstream turbulence. The height of the
spectrum now considerably increases if compared to the previous two cases.
Such an increase is primarily due to the sub-grid modeling adopted at the
shock front: the particle enters the shock with a pre-accelerated spectrum and
eventually ends up in the downstream region with a different steady state
(when compared to the STA alone case).

Finally, the particle energy evolution for the case in which STA is active
only in the downstream region is shown in bottom right panel. As expected, the
particle distribution does not significantly change until the particle crosses
the shock and then enters in the downstream region where turbulence is active.
Here steady state is attained due to STA. In this sense, the evolution
resembles the previous case.

Further notice that, for all the cases but the pure DSA one, the particle
distribution functions eventually seem to achieve steady states of similar
kind. This is expected as the predicted steady state spectrum depends on the
functional form of the transport coefficients which are not affected by the
presence of the shock.

####  5.1.1 Effect of turbulence on evolution of particle spectra

Additionally, in Fig. 7 we compare the particle steady-state distribution for
turbulent spectra with $q=5/3$ (middle), and with $q=3/2$ (right) with that
obtained for $q=2$ (left).

The main difference between the acceleration scenario for turbulent spectrum
with $q=2$, on one side, and $q=5/3$ or $q=3/2$, on the other, is that the
latter achieve steady state more rapidly because of the dependence of $t_{A}$
on $\gamma$.

Furthermore, the steady-state spectra for $q=5/3,3/2$ in the case of shock and
STA are not significantly different from the ones computed with STA alone (see
blue and orange solid line in the middle and right plot of Fig. 7). Owing to
the smaller acceleration timescale, in fact, the spectra for $q=5/3,3/2$
approach the steady state only when the particle arrives in the upstream
region making the shock injection less effective (see section 6) compared to
the $q=2$ case. However, for the case where turbulence is present only in the
downstream region, shock injection can clearly be observed (solid green line
in Fig. 7) as no significant turbulent energization took place in the upstream
region.

Additionally, we analyze the behaviour of $\gamma_{c}$, with various values of
$B_{0}$, $\rho_{u}$ and $\lambda_{\max}$. Analytically the value of
$\gamma_{c}$ can be calculated by equating $t_{A}$ to $t_{L}$ and yielding

$\gamma_{c}=\left\\{2\times
10^{3}\times\frac{\Big{(}\frac{eB\lambda_{\max}}{m_{e}c^{2}}\Big{)}^{2-q}}{\rho\lambda_{\max}}\right\\}^{\frac{1}{3-q}}$
(37)

Plots of $\gamma_{c}$ computed from simulation data with different values of
$B$ , $\rho$ and $\lambda_{\max}$ are compared in Fig. 8 toghether with the
analytic form (Eq. 37). We observe a good correspondence between the results.

####  5.1.2 Interplay of DSA and STA

Figure 9: Dependence of shock injection on the upstream spectrum for various
shock compression ratio with $\beta=100.0$. The obliquity is made fixed at
$30^{\circ}$. In the inset the downstream distribution function is shown for
two different values of $t_{A}/t_{L}$.

In the previous section we found that the shock acceleration depends on the
upstream spectrum. With this motivation here we try to analyze the impact of
STA on particle shock energization by modulating the acceleration timescale
$t_{A}$ and display its effect on the shock injection with different
compression ratios. Moreover, we define the value of $t_{A}$ in terms of
$t_{L}$ at $\gamma=1.0$ and for each choice of $t_{A}$, we perform the
simulation up to time $\tau=100\,\hat{\tau}_{0}$. Owing to the conserving
nature of the boundary condition, the number of micro-particles in a macro-
particle remains same once the shock takes place, thus by calculating the
number of micro-particles after shock we estimate the effect of shock
injection when STA is in process. The variation of total number of particles
after shock is shown with ratio $t_{A}/t_{L}$ at $\gamma=1.0$ for different
shock compression ratio in Fig. 9 with a fixed magnetic field calculated using
$\beta\,=\,100.0$. Further, the corresponding particle spectra at
$\tau=100\,\hat{\tau}_{0}$ is plotted for two values of the ratio and is shown
in the inset of Fig. 9.

When $t_{A}$ is much less than $t_{L}$ at $\gamma=1.0$ (or the ratio
$t_{A}/t_{L}$ is small) the particle spectrum reaches the log-parabolic
steady-state (see section. 5.1), before shock hits the particle. making the
shock injection less effective. On the other hand when the ratio $t_{A}/t_{L}$
is comparatively high, one observe very minute effect of STA on the particle
distribution in the upstream making the shock injection very effective for
this case. Furthermore, notice that for any value of $t_{A}/t_{L}$ shock with
higher compression ratio injects more number of particles than the lower ones.
Also from the distribution functions shown in the inset, for two different
values of $t_{A}/t_{L}$, it can be observed that the spectra that were hit by
strong shock (high compression ratio) reach to the steady state much faster
compared with the spectra hit by moderate shock (moderate compression ratio).
Moreover, the decrement of the $\gamma_{c}$ (see section 5.1.1) with
increasing $t_{A}/t_{L}$ could also be seen. Additionally, the number could be
seen to achieve a steady state, around $N\sim 10^{-6}$, at the higher values
of $t_{A}/t_{L}$ implies an upper bound of the particle injection at the shock
for different compression ratios.

In summary, we observe that the effect of shock injection on the particle
distribution function depends on the nature of the upstream particle
distribution spectra. If the timescale of the STA in the upstream region is
such that the particle distribution converges to steady-state spectra before
the DSA could take place, the effect of shock injection becomes minimal.
However, if in the upstream region the particle spectra do not reach the
steady-state before the shock hits the particle, then a considerable effect of
shock injection on particle spectra could be seen. This analysis spanning a
wide parameter base, therefore showcases the interplay of these two particle
acceleration processes.

###  5.2 Relativistic Blast Wave

Figure 10: Temporal evolution of particle distribution of a Lagrangian
particle in a turbulent medium for relativistic blast wave with different $B$
fields. The turbulent spectrum is taken as $\propto k^{-2}$, so the value of
$q$ is $2$ and the value of $\lambda_{\max}=\hat{L}_{0}/10$. _Left_ :
Corresponds to $B_{0}=5\times 10^{-2}\hat{B}_{0}$, _Middle_ : Depicts the
evolution of the particle distribution for $B_{0}=5\times 10^{-3}\hat{B}_{0}$
and _Right_ : Corresponds to the evolution for $B_{0}=5\times
10^{-4}\hat{B}_{0}$. Dashed blue line corresponds to the initial distribution
function which is $\propto\gamma^{-9}$. Figure 11: Spectral slope distribution
of particles initially placed at different angle ($\phi$) at the final time
($\tau=6$) with $B_{0}=5\times 10^{-4}\hat{B}_{0}$ for the relativistic blast
wave test.

Here we focus on the impact of a relativistic blast wave on the evolution of
the spectral distribution in the presence of both shock and turbulence. Due to
the underlying symmetry of the problem we choose a single quadrant with
$512^{2}$ Cartesian computational zones with $x,y\in[0,6]$. The initial
condition consists of an over-pressurized central region of circular radius
$0.8\hat{L}_{0}$ filled with pressure and density
$\\{P_{c},\rho_{c}\\}=\\{1,1\\}$ surrounded by a uniform medium with
$\\{P_{e},\rho_{e}\\}=\\{3\times 10^{-5},10^{-2}\\}$. The magnetic field is
taken perpendicular to the $\\{x,y\\}$ plane, $\boldsymbol{B}=B_{0}\hat{z}$ as
in Vaidya et al. (2018). The boundary condition is set to be reflecting at
$x=y=0$ and outflow elsewhere. We initially place $360$ Lagrangian macro-
particles uniformly over $0<\phi<\pi/2$ at the radius of
$\sqrt{x^{2}+y^{2}}=2$. Physical units are chosen such that
$\hat{L}_{0}=10\,{\rm pc}$, $\hat{\rho}_{0}=0.01\,{\rm amu}$,
$\hat{P}_{0}=1.5\times 10^{-5}\,{\rm dyne/cm}^{2}$, $\hat{v}_{0}=c$,
$\hat{B}_{0}=1.37\times 10^{-2}\,{\rm G}$ and $\hat{\tau}_{0}=32.64\,{\rm
yrs}$. The initial distribution function for each macro-particle is taken to
be a steep decreasing power-law profile with index $9$ covering a range in
Lorentz factor $\gamma\in\\{1,10^{8}\\}$ discretized using $128$ bins. Similar
to the MHD planar shock test (section 5.1), the diffusion coefficient is
modelled following the acceleration timescale and the losses are modelled
following the synchrotron, Inverse-Compton and adiabatic loss processes.

The evolution of the particle distribution for a macro-particle initially
placed at $65^{\circ}$, for $q=2$, is shown in Fig. 10, where the particle
evolution is shown for $3$ different magnetic fields: $B_{0}=5\times 10^{-2}$
(left panel), $B_{0}=5\times 10^{-3}$ (middle panel) and $B_{0}=5\times
10^{-4}$ (right panel). Furthermore, in all three cases the value of
$\lambda_{\max}=\hat{L}_{0}/10$.

For the case with strongest magnetic field, the particle distribution
initially evolves due to STA and, after crossing the shock, a steady-state
ultra-relativistic Maxwellian-like spectral distribution can be seen to emerge
eventually with a sharp cut-off beyond $\gamma_{c}\sim 10^{8}$. On the
contrary, for the weakest magnetic field case, the spectral evolution shows
distinct signatures of DSA only. Indeed, STA signature can hardly be observed
as the timescale obeys $t_{A}\propto B^{-2}$ (see Eq. 12), thus very large for
the simulation time. In this case, the initial steep spectra is accelerated
and the spectral slope is flattened and cooling due to synchrotron and IC
emission is evident from the cut-off. Moreover, it should be noted that the
particle can be energized beyond $\gamma>10^{9}$. For the intermediate case,
we observe effects of both shock and STA in shaping the particle spectra.

Additionally, we quantified grid orientation effects by estimating the slope
of the distribution functions for each macro-particle as a function of their
initial angular positions. This is shown, at time $\tau=6$ for $B_{0}=5\times
10^{-4}\hat{B}_{0}$, in Fig. 11. The final slope for all the macro-particles
approximately fall in the same range ($\approx-4$) with additional variations
due to discretization error ($\sim 2\%$). Therefore all macro-particles will
have similar spectral distribution as shown for the typical macro-particle in
Fig. 10, apart from the minor variations due to discretization error.

###  5.3 Relativistic Magneto-hydrodynamic Jet

Figure 12: Temporal evolution of the spectrum of a Lagrangian particle which
has gone through shock atleast once, in the RMHD Jet. _Top:_ For the case of
only DSA _Bottom:_ For the case with STA along with DSA. Figure 13: Comparison
between the emission from turbulence and DSA and only DSA for radio frequency,
$1.4\,GHz$ at time $\tau=200$. Notice that the radial coordinate has been
mirrored in the left plot. Figure 14: Same as Fig. 13 but for optical blue
light frequency $6.59\times 10^{5}\,GHz$ at time $\tau=200$. Figure 15: Same
as Fig. 13 but for $0.4\,KeV$ X-Ray at time $\tau=200$.

In this section, we describe a toy model of a relativistic magneto-
hydrodynamic jet and analyze its emission signatures due to the DSA and STA of
cosmic rays. In particular, we employ a $2D$ cylindrical grid
$\\{R,Z\\}\in\\{0,0\\}$ to $\\{20,50\\}$ using $160\times 400$ grid cells. The
ambient medium is initially static ($\boldsymbol{V}_{m}=0$) with constant
density $\rho_{m}=10^{3}\hat{\rho}_{0}$, where, $\hat{\rho}_{0}=1.67\times
10^{-24}\,{\rm gr}\,{\rm cm}^{-3}$. An under-dense beam with
$\rho_{j}=\hat{\rho}_{0}$ is injected into the ambient medium with velocity
$v_{z}$ along the vertical direction through a circular nozzle of unit radius,
$R_{j}=\hat{L}_{0}$ from the lower $Z$ boundary. The value of $v_{z}$ is
prescribed using the Lorentz factor $\gamma_{j}=10$ and $\hat{L}_{0}=100\,{\rm
pc}$ implying an unit timescale of $\hat{\tau}_{0}=326.4\,{\rm yrs}$. The
magnetic field is purely poloidal,
$\boldsymbol{B}=B_{z}\hat{\boldsymbol{e}}_{z}$ and is initially prescribed in
jet nozzle and also in the ambient medium,

$B_{z}=\sqrt{2\sigma_{z}P_{j}}.$ (38)

where, $P_{j}$ is the jet pressure at $R=R_{j}$ estimated from the Mach number
$M=v_{j}\sqrt{\rho_{j}/(\Gamma P_{j})+1/(\Gamma-1)}=6$ and adiabatic index
$\Gamma=5/3$. The values for $\sigma_{z}$ is taken to be $10^{-4}$ for the
present simulation.

We further inject $25$ Lagrangian macro-particles every two time steps with an
initial power-law spectral distribution with index $-9$ on a initial $\gamma$
grid with $\\{\gamma_{\min},\gamma_{\max}\\}\equiv\\{1,10^{5}\\}$ discretized
with $128$ bins.

The energy spectrum of the macro-particles are calculated for two different
scenarios: i) considering only DSA and different losses and ii) considering,
in addition, also stochastic processes. For scenario (i) we follow the
numerical algorithm developed in Vaidya et al. (2018); Mukherjee et al. (2021)
to estimate the particle spectral distribution, while for scenario (ii) we
solve Eq. (6) without the source and sink terms, along with the diffusion
coefficient $D\propto\gamma^{2}$, where the proportionality constant is
computed from the value of $t_{A}$ following Eq. (12) and with the value of
$\lambda_{\max}=\hat{L}_{0}/100$. The loss terms account for synchrotron,
Inverse Compton and adiabatic losses. Also, compared to the previous test
problems here we take Courant number $0.8$ when solving Eq. (6). Moreover, for
both scenarios we compute the emissivity for each macro-particle based on
their local spectral distribution and interpolated it on the underlying grid
(Vaidya et al., 2018).

In Fig. 12, we show the spectral evolution of representative particles, that
have been shocked at least once, for each of the scenarios. The top panel
shows spectral evolution of a representative particle for the case where
acceleration is due to shocks alone. The effect of DSA and radiative losses
are clearly visible, respectively, from the spectral flattening and from high
energy cut-offs. Here the cut-off can be observed clearly, as during DSA, the
maximum energy get shifted according to the prescription described in Sec. 2.
When the maximum $\gamma$ exceeds its initial value, cooling processes become
effective so that the macro-particle quickly cools accounting for sharp
spectral cut-off.

The bottom panel shows the spectral evolution of similar particle for the case
where STA is also included (besides DSA). the distribution reveals a hump-like
structure in the low-energy end of the spectrum that slowly shifts towards
higher $\gamma$ values. With time, this eventually leads the distribution
function to reach a steady state, as described by Eq. (A2). Notice that our
choice of parameters (Eq. 12) is such that the acceleration timescale $t_{A}$
is larger or comparable to the dynamical time, leading to feeble acceleration.
We also point out that, during the initial stages, the particle spectrum
exhibits a pile-up effect at low $\gamma$, because of the finite grid
constraint, as discussed in section 4.3. This spurious effect dims with time
as lower $\gamma$ particles starts to accelerate toward higher $\gamma$. The
impact of DSA (in addition to STA) can be distinguished from the flattening of
the spectral distribution. The more pronounced low-energy cutoff is attributed
to the lower energy particles being accelerated by STA, eventually creating a
deficiency in the number of particles at low $\gamma$.

From the instantaneous spectral distribution of Lagrangian macro-particles
spread across the computational domain, we estimate the synchrotron emissivity
by convolving the macro-particle spectra with single electron synchrotron
spectra and interpolated it on the computational grid (see Eq. 36-37 in Vaidya
et al., 2018). In Figs. 13, 14 and 15, the emissivity $J_{\nu}$ computed from
the Lagrangian macro-particles is shown for different frequencies at time
$\tau=200\hat{\tau}_{0}$ for the two different scenarios (left and right
halves, respectively).

In Fig. 13, with $1.4\,{\rm GHz}$ radio frequency, the emission due to
turbulence and shock (right half) is very similar to the case with DSA only
(left half). For the case with optical frequency ($\nu=6.59\times 10^{5}$ GHz)
(Fig. 14), the emission becomes less than the radio frequency (Fig. 13) for
both cases with and without STA. This is expected because of the faster
cooling time with higher energy. However, a significant larger emission can be
seen in case ii) in the region $Z\lesssim 10$. The material in this region
originates from the back-flow dynamics of the jet (Cielo et al., 2014;
Matthews, James H. et al., 2019). If only shock energization is accounted for,
the particle spectra become very steep in this region owing to radiative
losses and the absence of strong shocks. However, if STA is also taken into
account, the spectra remain hard because of the competing effects of STA and
radiative losses. Similar high emission features are observed in X-ray
($\nu=10^{8}$ GHz) as well (right panel of Fig. 15). On the contrary, in the
presence of DSA only, a significant reduction in the X-ray emission can be
seen (left half). Here most of the emission originates from the regions near
jet head as well as isolated spots in the cocoon. In addition, smaller
emission centers can be observed in the region around the re-collimation
shocks along the beam. This differs from the case with DSA + STA, where the
emission pattern was wider and more uniformly distributed throughout the jet
and the backflow region.

##  6 Discussion and Summary

In this paper we have focused on the numerical modeling of stochastic
turbulent acceleration (STA) and its physical contribution to the spectral
evolution of highly energetic particles. The numerical formulation is based on
the fluid-particle hybrid framework of Vaidya et al. (2018); Mukherjee et al.
(2021) developed for the PLUTO code, where the non-thermal plasma component is
modeled by means of Lagrangian macro-particles embedded in a classical or
relativistic magnetized thermal flow.

The particle distribution function is evolved by solving numerically a Fokker-
Planck equation in which STA is modelled by two components: a hyperbolic term
describing the systematic acceleration (Fermi $\rm II$) and a parabolic
contribution accounting for random resonant interaction between particles and
plasma turbulent waves. While Vaidya et al. (2018) presented a Lagrangian
method for the solution of the Fokker-Planck equation in the presence of
hyperbolic terms only, here we have introduced a novel Eulerian algorithm to
account also for an energy-dependent diffusion coefficient $D\sim\gamma^{2}$
which can become stiff in the high-energy limit. To overcome the explicit time
step restriction, the new method takes advantage of $2^{\rm nd}$-order Runge
Kutta Implicit-Explicit (IMEX) methods, so that hyperbolic terms (e.g.
adiabatic expansion / radiative losses / Fermi $\rm II$) are treated
explicitly while parabolic terms (modelling turbulent diffusion) are handled
implicitly.

Selected numerical benchmarks validated against analytical solutions and grid
resolution studies demonstrate that our implementation has improved stability
and accuracy properties when compared to previous solvers (see for example
Chang & Cooper, 1970; Winner et al., 2019). In addition, due to the presence
of boundary condition our algorithm respects physical constraints (for
example, $\gamma\geq 1$) which are not always satisfied in the Lagrangian
method (Vaidya et al., 2018; Mukherjee et al., 2021) with an evolving grid.
STA modeling has also been validated against radiative synchrotron loss
process by studying the evolution of curvature of particle spectrum (Tramacere
et al., 2011).

With these motivations, we have studied the effect of STA as well as other
energization processes, on the particle spectrum in the presence of shocks,
using toy-model applications. Such an interplay is commonly believed to
operate in supernova remnants, AGN radio lobes, galaxy clusters and radio
relics.

As a first application example, we considered a simple planar shock in four
different acceleration scenarios. We found that when STA and DSA both are
considered, the former seems to affect the shock injection by changing the
macro-particle distribution function. Further tests with different forms of
the diffusion coefficient reveal a similar behavior. Additionally, we have
also quantified the effect of STA time scale on the radiative losses and its
influence on the interplay with DSA. In particular, we observe that the effect
of shocks on particle distribution weakens with decreasing STA time scales.
Similar interplay of DSA and STA was also evident in case of spherical shock
formed in the test case of RMHD blast wave.

Finally, we have extended our algorithm to explore the emission properties of
the axisymmetric RMHD jet using a toy model. We find a significant difference
both in the evolution of the spectral distribution and the ensuing emission
signatures due to the presence or absence of the STA process. In particular,
inclusion of STA results in diffuse emission within the jet back-flow,
particularly in the high-energy X-ray band. Consequences of such an important
finding will be further explored in forthcoming works focusing on
astrophysical systems along with comparison with observed signatures.

## ACKNOWLEDGEMENTS

We would like to thank the anonymous referee for the helpful comments, and
constructive remarks on this manuscript. All simulations were performed at the
computing facility at Indian Institute of Technology, Indore. We would like to
thank the financial support from the Max Planck partner group award at Indian
Institute of Technology, Indore.

## References

  * Achterberg & Krulls (1992) Achterberg, A., & Krulls, W. M. 1992, A&A, 265, L13
  * Asano & Hayashida (2018) Asano, K., & Hayashida, M. 2018, ApJ, 861, 31, doi: 10.3847/1538-4357/aac82a
  * Ascher et al. (1997) Ascher, U. M., Ruuth, S. J., & Spiteri, R. J. 1997, Applied Numerical Mathematics, 25, 151, doi: https://doi.org/10.1016/S0168-9274(97)00056-1
  * Bai et al. (2015) Bai, X.-N., Caprioli, D., Sironi, L., & Spitkovsky, A. 2015, ApJ, 809, 55, doi: 10.1088/0004-637X/809/1/55
  * Baring (2004) Baring, M. G. 2004, Nuclear Physics B - Proceedings Supplements, 136, 198, doi: https://doi.org/10.1016/j.nuclphysbps.2004.10.008
  * Baring et al. (1994) Baring, M. G., Ellison, D. C., & Jones, F. C. 1994, International Astronomical Union Colloquium, 142, 547–552, doi: 10.1017/S0252921100077794
  * Bell (1978) Bell, A. R. 1978, MNRAS, 182, 147, doi: 10.1093/mnras/182.2.147
  * Bell (2014) Bell, A. R. 2014, Brazilian Journal of Physics, 44, 415.  https://doi.org/10.1007/s13538-014-0219-5
  * Belmont et al. (2008) Belmont, R., Malzac, J., & Marcowith, A. 2008, A&A, 491, 617, doi: 10.1051/0004-6361:200809982
  * Blandford & Eichler (1987) Blandford, R., & Eichler, D. 1987, Physics Reports, 154, 1 , doi: https://doi.org/10.1016/0370-1573(87)90134-7
  * Blandford (1994) Blandford, R. D. 1994, ApJS, 90, 515, doi: 10.1086/191869
  * Borse, Nikhil et al. (2021) Borse, Nikhil, Acharya, Sriyasriti, Vaidya, Bhargav, et al. 2021, A&A, 649, A150, doi: 10.1051/0004-6361/202140440
  * Böttcher & Dermer (2010) Böttcher, M., & Dermer, C. D. 2010, ApJ, 711, 445, doi: 10.1088/0004-637X/711/1/445
  * Brunetti & Lazarian (2007) Brunetti, G., & Lazarian, A. 2007, MNRAS, 378, 245, doi: 10.1111/j.1365-2966.2007.11771.x
  * Bykov & Fleishman (1992) Bykov, A. M., & Fleishman, G. D. 1992, Monthly Notices of the Royal Astronomical Society, 255, 269, doi: 10.1093/mnras/255.2.269
  * Chang & Cooper (1970) Chang, J., & Cooper, G. 1970, Journal of Computational Physics, 6, 1 , doi: https://doi.org/10.1016/0021-9991(70)90001-X
  * Cielo et al. (2014) Cielo, S., Antonuccio-Delogu, V., Macciò, A. V., Romeo, A. D., & Silk, J. 2014, Monthly Notices of the Royal Astronomical Society, 439, 2903, doi: 10.1093/mnras/stu161
  * Comisso & Sironi (2018) Comisso, L., & Sironi, L. 2018, Phys. Rev. Lett., 121, 255101, doi: 10.1103/PhysRevLett.121.255101
  * Dermer et al. (1996) Dermer, C. D., Miller, J. A., & Li, H. 1996, ApJ, 456, 106, doi: 10.1086/176631
  * Donnert & Brunetti (2014) Donnert, J., & Brunetti, G. 2014, Monthly Notices of the Royal Astronomical Society, 443, 3564, doi: 10.1093/mnras/stu1417
  * Drury (1983) Drury, L. O. 1983, Reports on Progress in Physics, 46, 973, doi: 10.1088/0034-4885/46/8/002
  * Ellison & Double (2002) Ellison, D. C., & Double, G. P. 2002, Astroparticle Physics, 18, 213, doi: https://doi.org/10.1016/S0927-6505(02)00142-1
  * Ellison et al. (1990) Ellison, D. C., Jones, F. C., & Reynolds, S. P. 1990, ApJ, 360, 702, doi: 10.1086/169156
  * Ellison et al. (2013) Ellison, D. C., Warren, D. C., & Bykov, A. M. 2013, The Astrophysical Journal, 776, 46, doi: 10.1088/0004-637x/776/1/46
  * Fan et al. (2008) Fan, Z.-H., Liu, S., Wang, J.-M., Fryer, C. L., & Li, H. 2008, The Astrophysical Journal, 673, L139, doi: 10.1086/528372
  * Fermi (1949) Fermi, E. 1949, Phys. Rev., 75, 1169, doi: 10.1103/PhysRev.75.1169
  * Ferrand & Marcowith (2010) Ferrand, G., & Marcowith, A. 2010, A&A, 510, A101, doi: 10.1051/0004-6361/200913520
  * Giacalone & Ellison (2000) Giacalone, J., & Ellison, D. C. 2000, Journal of Geophysical Research: Space Physics, 105, 12541, doi: https://doi.org/10.1029/1999JA000018
  * Giacalone & Jokipii (1999) Giacalone, J., & Jokipii, J. R. 1999, ApJ, 520, 204, doi: 10.1086/307452
  * Jokipii et al. (2007) Jokipii, J., Giacalone, J., & Kóta, J. 2007, Planetary and Space Science, 55, 2267, doi: https://doi.org/10.1016/j.pss.2007.05.007
  * Kang et al. (2017) Kang, H., Ryu, D., & Jones, T. W. 2017, The Astrophysical Journal, 840, 42, doi: 10.3847/1538-4357/aa6d0d
  * Kardashev (1962) Kardashev, N. S. 1962, Soviet Ast., 6, 317
  * Katarzyński et al. (2006) Katarzyński, K., Ghisellini, G., Mastichiadis, A., Tavecchio, F., & Maraschi, L. 2006, A&A, 453, 47, doi: 10.1051/0004-6361:20054176
  * Kirk et al. (1996) Kirk, J. G., Duffy, P., & Gallant, Y. A. 1996, A&A, 314, 1010.  https://arxiv.org/abs/astro-ph/9604056
  * Kirk et al. (1994) Kirk, J. G., Melrose, D. B., & Priest, E. R. 1994, Plasma Astrophysics, ed. A. O. Benz, , & T. J.-L. Courvoisier (Springer Berlin Heidelberg), doi: 10.1007/3-540-31627-2
  * Kulsrud & Ferrari (1971) Kulsrud, R. M., & Ferrari, A. 1971, Ap&SS, 12, 302, doi: 10.1007/BF00651420
  * Lemoine & Pelletier (2003) Lemoine, M., & Pelletier, G. 2003, The Astrophysical Journal, 589, L73, doi: 10.1086/376353
  * Liu et al. (2004) Liu, S., Petrosian, V., & Melia, F. 2004, ApJ, 611, L101, doi: 10.1086/423985
  * Malkov & Drury (2001) Malkov, M. A., & Drury, L. O. 2001, Reports on Progress in Physics, 64, 429, doi: 10.1088/0034-4885/64/4/201
  * Marcowith & Casse (2010) Marcowith, A., & Casse, F. 2010, A&A, 515, A90, doi: 10.1051/0004-6361/200913022
  * Marcowith et al. (2020) Marcowith, A., Ferrand, G., Grech, M., et al. 2020, Living Reviews in Computational Astrophysics, 6, 1, doi: 10.1007/s41115-020-0007-6
  * Marcowith & Kirk (1999) Marcowith, A., & Kirk, J. G. 1999, A&A, 347, 391.  https://arxiv.org/abs/astro-ph/9905176
  * Massaro et al. (2006) Massaro, E., Tramacere, A., Perri, M., Giommi, P., & Tosti, G. 2006, A&A, 448, 861, doi: 10.1051/0004-6361:20053644
  * Massaro, E. et al. (2004) Massaro, E., Perri, M., Giommi, P., & Nesci, R. 2004, A&A, 413, 489, doi: 10.1051/0004-6361:20031558
  * Matthews, James H. et al. (2019) Matthews, James H., Bell, Anthony R., Araudo, Anabella T., & Blundell, Katherine M. 2019, EPJ Web Conf., 210, 04002, doi: 10.1051/epjconf/201921004002
  * Meli, A. & Biermann, P. L. (2013) Meli, A., & Biermann, P. L. 2013, A&A, 556, A88, doi: 10.1051/0004-6361/201016299
  * Melrose (1996) Melrose, D. B. 1996, Astrophysics and Space Science, 242, 209.  https://doi.org/10.1007/BF00645114
  * Mertsch & Petrosian (2019) Mertsch, P., & Petrosian, V. 2019, A&A, 622, A203, doi: 10.1051/0004-6361/201833999
  * Mignone et al. (2007) Mignone, A., Bodo, G., Massaglia, S., et al. 2007, The Astrophysical Journal Supplement Series, 170, 228, doi: 10.1086/513316
  * Mignone et al. (2018) Mignone, A., Bodo, G., Vaidya, B., & Mattia, G. 2018, The Astrophysical Journal, 859, 13, doi: 10.3847/1538-4357/aabccd
  * Mimica & Aloy (2012) Mimica, P., & Aloy, M. A. 2012, Monthly Notices of the Royal Astronomical Society, 421, 2635, doi: https://doi.org/10.1111/j.1365-2966.2012.20495.x
  * Miniati (2001) Miniati, F. 2001, Computer Physics Communications, 141, 17, doi: https://doi.org/10.1016/S0010-4655(01)00293-4
  * Miniati (2003) Miniati, F. 2003, MNRAS, 342, 1009, doi: 10.1046/j.1365-8711.2003.06647.x
  * Miniati et al. (2001) Miniati, F., Ryu, D., Kang, H., & Jones, T. W. 2001, The Astrophysical Journal, 559, 59, doi: 10.1086/322375
  * Mukherjee et al. (2021) Mukherjee, D., Bodo, G., Rossi, P., Mignone, A., & Vaidya, B. 2021, Monthly Notices of the Royal Astronomical Society, 505, 2267, doi: 10.1093/mnras/stab1327
  * Niemiec & Ostrowski (2006) Niemiec, J., & Ostrowski, M. 2006, The Astrophysical Journal, 641, 984, doi: 10.1086/500541
  * Nishikawa et al. (2007) Nishikawa, K.-I., Hededal, C. B., Hardee, P. E., et al. 2007, Astrophysics and Space Science, 307, 319, doi: 10.1007/s10509-006-9234-5
  * Ostrowski (1988) Ostrowski, M. 1988, Monthly Notices of the Royal Astronomical Society, 233, 257, doi: 10.1093/mnras/233.2.257
  * O’Sullivan et al. (2009) O’Sullivan, S., Reville, B., & Taylor, A. M. 2009, Monthly Notices of the Royal Astronomical Society, 400, 248, doi: 10.1111/j.1365-2966.2009.15442.x
  * Palmroth et al. (2018) Palmroth, M., Ganse, U., Pfau-Kempf, Y., et al. 2018, Living Reviews in Computational Astrophysics, 4, 1, doi: 10.1007/s41115-018-0003-2
  * Pareschi & Russo (2005) Pareschi, L., & Russo, G. 2005, Journal of Scientific Computing, 25, 129.  https://doi.org/10.1007/BF02728986
  * Park & Petrosian (1995) Park, B. T., & Petrosian, V. 1995, ApJ, 446, 699, doi: 10.1086/175828
  * Park & Petrosian (1996) —. 1996, ApJS, 103, 255, doi: 10.1086/192278
  * Parker (1955) Parker, E. N. 1955, Physical Review, 99, 241, doi: 10.1103/PhysRev.99.241
  * Perri & Zimbardo (2015) Perri, S., & Zimbardo, G. 2015, ApJ, 815, 75, doi: 10.1088/0004-637X/815/1/75
  * Petrosian (2012) Petrosian, V. 2012, Space Sci. Rev., 173, 535, doi: 10.1007/s11214-012-9900-6
  * Press et al. (1992) Press, W. H., Teukolsky, S. A., Vetterling, W. T., & Flannery, B. P. 1992, Numerical recipes in C. The art of scientific computing
  * Rajpurohit et al. (2020) Rajpurohit, K., Hoeft, M., Vazza, F., et al. 2020, A&A, 636, A30, doi: 10.1051/0004-6361/201937139
  * Rangelov et al. (2019) Rangelov, B., Montmerle, T., Federman, S. R., Boissé, P., & Gabici, S. 2019, The Astrophysical Journal, 885, 105, doi: 10.3847/1538-4357/ab43e5
  * Rieger & Duffy (2019) Rieger, F. M., & Duffy, P. 2019, The Astrophysical Journal, 886, L26, doi: 10.3847/2041-8213/ab563f
  * Sarazin (1999) Sarazin, C. L. 1999, The Astrophysical Journal, 520, 529, doi: 10.1086/307501
  * Schlickeiser (1989) Schlickeiser, R. 1989, ApJ, 336, 243, doi: 10.1086/167009
  * Schlickeiser (2002) —. 2002, Cosmic Ray Astrophysics
  * Schlickeiser & Dermer (2000) Schlickeiser, R., & Dermer, C. D. 2000, A&A, 360, 789.  https://arxiv.org/abs/astro-ph/0005601
  * Sironi & Spitkovsky (2011) Sironi, L., & Spitkovsky, A. 2011, ApJ, 726, 75, doi: 10.1088/0004-637X/726/2/75
  * Spitkovsky (2008) Spitkovsky, A. 2008, The Astrophysical Journal, 682, L5, doi: 10.1086/590248
  * Sturrock (1966) Sturrock, P. A. 1966, Phys. Rev., 141, 186, doi: 10.1103/PhysRev.141.186
  * Teraki & Asano (2019) Teraki, Y., & Asano, K. 2019, The Astrophysical Journal, 877, 71, doi: 10.3847/1538-4357/ab1b13
  * Tramacere et al. (2011) Tramacere, A., Massaro, E., & Taylor, A. M. 2011, The Astrophysical Journal, 739, 66, doi: 10.1088/0004-637x/739/2/66
  * Tverskoǐ (1967) Tverskoǐ, B. A. 1967, Soviet Journal of Experimental and Theoretical Physics, 25, 317
  * Vaidya et al. (2018) Vaidya, B., Mignone, A., Bodo, G., Rossi, P., & Massaglia, S. 2018, The Astrophysical Journal, 865, 144, doi: 10.3847/1538-4357/aadd17
  * Vaidya et al. (2017) Vaidya, B., Prasad, D., Mignone, A., Sharma, P., & Rickler, L. 2017, Monthly Notices of the Royal Astronomical Society, 472, 3147, doi: 10.1093/mnras/stx2176
  * Van Leer (1977) Van Leer, B. 1977, Journal of Computational Physics, 23, 276, doi: https://doi.org/10.1016/0021-9991(77)90095-X
  * van Weeren et al. (2017) van Weeren, R. J., Andrade-Santos, F., Dawson, W. A., et al. 2017, Nature Astronomy, 1, 0005, doi: 10.1038/s41550-016-0005
  * Vazza et al. (2021) Vazza, F., Wittor, D., Brunetti, G., & Brüggen, M. 2021, arXiv e-prints, arXiv:2102.04193.  https://arxiv.org/abs/2102.04193
  * Vurm & Poutanen (2009) Vurm, I., & Poutanen, J. 2009, The Astrophysical Journal, 698, 293, doi: 10.1088/0004-637x/698/1/293
  * Webb (1989) Webb, G. M. 1989, ApJ, 340, 1112, doi: 10.1086/167462
  * Winner et al. (2019) Winner, G., Pfrommer, C., Girichidis, P., & Pakmor, R. 2019, MNRAS, 488, 2235, doi: 10.1093/mnras/stz1792
  * Wolff & Tautz (2015) Wolff, M., & Tautz, R. C. 2015, A&A, 580, A58, doi: 10.1051/0004-6361/201525907
  * Wong et al. (2019) Wong, K., Zhdankin, V., Uzdensky, D. A., Werner, G. R., & Begelman, M. C. 2019, arXiv e-prints, arXiv:1901.03439.  https://arxiv.org/abs/1901.03439
  * Zimbardo & Perri (2017) Zimbardo, G., & Perri, S. 2017, Nature Astronomy, 1, 0163.  https://doi.org/10.1038/s41550-017-0163

##  Appendix A Analytical solution of Fokker-Planck Equation

Eq. (6) is very hard to solve for a proper general analytic solution. Various
work has been devoted to solve Eq. (6) for various transport coefficients
(e.g., Katarzyński et al., 2006; Park & Petrosian, 1995; Chang & Cooper, 1970;
Kardashev, 1962). Chang & Cooper (1970) solved Eq. (6) for the steady-state
solution and the solution could be written as,

$\displaystyle\chi_{\rm steady}(\gamma)=\chi_{\rm
0}\exp\Big{\\{}-\int^{\gamma}_{\rm
1}\Big{(}\frac{S(\gamma^{\prime},\tau)-D_{\rm
A}(\gamma^{\prime},\tau)}{D_{\gamma\gamma}(\gamma^{\prime},\tau)}\Big{)}d\gamma^{\prime}\Big{\\}}.$
(A1)

Katarzyński et al. (2006) solved Eq. (A1) for
$D_{\gamma\gamma}(\gamma,\tau)=D_{\gamma\rm 0}\gamma^{\rm 2}/2$ with
$D_{\gamma\rm 0}=1/t_{\rm A}$, $D_{\rm A}(\gamma,\tau)=\gamma/t_{\rm A}$ and
$S(\gamma,\tau)=S_{\rm 0}\gamma^{\rm 2}$. These form of the parameters are
typical for particles in plasma. The loss term $S(\gamma,\tau)$ gets a similar
form if Inverse-Compton radiation is taken in the Thompson limit with
Synchrotron radiation and the form for the diffusion coefficient
$D_{\gamma\gamma}$ which also matches the form from typical particle in cell
simulation as discussed above. The solution to Eq. (A1) with the above
mentioned parameters is,

$\displaystyle\chi_{\rm steady}(\gamma)=\chi_{\rm 0}\gamma^{\rm
2}\exp\\{-2S_{\rm 0}t_{\rm A}(\gamma-1)\\}.$ (A2)

Kardashev (1962) got a time-dependent solution for Eq. (6) without the loss
terms and showed the acceleration leads to a log-normal particle distribution
(similar to Eq. (36)).

So, if the particles only accelerate via STA the particle distribution follows
a log-normal form due to the fact that the STA process is a multiplicative
acceleration process (Tramacere et al., 2011). But if those particles loose
their energy via radiative means along with the acceleration the particle
distribution starts to follow an ultra-relativistic Maxwellian (Eq. (A2)),
which looks like a thermal or quasi-thermal spectrum with a scaled temperature
of $1/S_{\rm 0}t_{\rm A}$ which is also the value of $\gamma$ where, $t_{\rm
A}=t_{\rm L}$.