11institutetext: Discipline of Astronomy, Astrophysics and Space Engineering,
Indian Institute of Technology, Indore, Madhya Pradesh, India - 452020
22institutetext: Dipartimento di Fisica Generale, Universita degli Studi di
Torino, Via Pietro Giuria 1, 10125 Torino, Italy 33institutetext: Centre for
Astrophysics Research, Department of Physics, Astronomy and Mathematics,
University of Hertfordshire, College Lane, Hatfield AL10 9AB, UK

# A numerical study of the interplay between Fermi acceleration mechanisms in
radio lobes of FR-II radio galaxies

Sayan Kundu<EMAIL_ADDRESS>Bhargav Vaidya 11 Andrea Mignone 22
Martin J. Hardcastle 33

###### Abstract

Context. Radio-loud AGNs are thought to possess various sites of particle
acceleration, which gives rise to the observed non-thermal spectra. Stochastic
turbulent acceleration (STA) and diffusive shock acceleration (DSA) are
commonly cited as potential sources of high-energy particles in weakly
magnetized environments. Together, these acceleration processes and various
radiative losses determine the emission characteristics of these extra-
galactic radio sources.

Aims. The purpose of this research is to investigate the dynamical interplay
between the STA and DSA in the radio lobes of FR-II radio galaxies, as well as
the manner in which these acceleration mechanisms, along with a variety of
radiative losses, collectively shape the emission features seen in these
extra-galactic sources.

Methods. A phenomenologically motivated model of STA is considered and
subsequently employed on a magneto-hydrodynamically simulated radio lobe
through a novel hybrid Eulerian-Lagrangian framework.

Results. STA gives rise to a curved particle spectrum that is morphologically
different from the usual shock-accelerated spectrum. As a consequence of this
structural difference in the underlying particle energy spectrum, various
multi-wavelength features arise in the spectral energy distribution of the
radio lobe. Additionally, we observe enhanced diffuse X-ray emission from
radio lobes for cases where STA is taken into account in addition to DSA.

###### Key Words.:

Magnetohydrodynamics (MHD) – Methods: numerical – Acceleration of particles –
Radio continuum: galaxies – X-rays: galaxies – Turbulence

##  1 Introduction

Radio galaxies are considered one of the most energetic systems in the
universe. These extra-galactic objects are observed to possess a huge
reservoir of relativistic non-thermal particles, which collectively shape
their emission features (Blandford et al., 2019). Further, due to the
abundance of highly energetic particles, these galaxies are generally
considered a favourable site to study various high-energy phenomena
(Meisenheimer, 2003). In recent years, with the advent of multi-messenger
astronomy, different observations are uncovering various features and helping
us understand different micro-physical processes happening in these systems
(Marcowith et al., 2020).

Low frequency radio observations of these radio galaxies provide insights
about their morphological structures (see Hardcastle & Croston, 2020, for more
details), their magnetic field strength (Croston et al., 2005a) and their age
(Alexander & Leahy, 1987; Carilli et al., 1991; Mahatma et al., 2019). Based
on the brightness of these sources at $178$ MHz, they are classified as
Fanaroff Rilley (FR) class I (low-power) or II (high-power) (Fanaroff & Riley,
1974). These two classes of radio galaxies are observed to manifest different
morphological structures. While FR-II sources exhibit a one-sided smooth
spine-like structure with a bright termination point, FR-I sources show a two-
sided plume-like structure. Additionally, FR-II sources show prominent signs
of turbulent cocoons that have an extent of a few hundred kiloparsecs and are
often partly visible as lobes (Mullin et al., 2008; Hardcastle & Croston,
2020). These lobes are believed to be highly magnetized cavities of rarefied
plasma where most of the jet kinetic power is deposited. Radio lobes also have
a hotspot region near the jet termination region, responsible for accelerating
particles to high energies via diffusive shock acceleration (DSA) (Brunetti et
al., 2001; Prieto et al., 2002; Araudo et al., 2018). These freshly shock
accelerated particles further get mixed with the older plasma particles,
already residing in the lobe, which makes the lobe a turbulent playground for
various plasma waves to interact with the particles and eventually accelerate
them via stochastic turbulent acceleration (STA). Such a mechanism has also
been invoked to explain the particle acceleration in various astrophysical
systems such as 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; Tavecchio et al., 2022, and references therein), Fermi
bubbles (Mertsch & Petrosian, 2019), galaxy clusters (Brunetti & Lazarian,
2007; Donnert & Brunetti, 2014; Vazza et al., 2021). STA has been invoked as a
possible mechanism for producing Ultra High Energy Cosmic Rays (UHECRs) from
the radio lobe of Pictor A (Fan et al., 2008) and Cen A (Hardcastle et al.,
2009; O’Sullivan et al., 2009). Also, recently, it has been invoked as a
plausible candidate in explaining the spectral curvature usually observed in
FR-II radio lobes (Harris et al., 2019).

In addition to the radio observations, X-ray observations of these radio loud
AGNs have become popular due to the minimal contamination of the X-ray
radiation by non-AGN sources. Several components of these sources, such as
radio lobes, hotspots, and collimated radio jet spine, are observed to radiate
in the X-ray band (de Vries et al., 2018; Massaro et al., 2018). Additionally,
these lobes are often observed to give rise to diffuse X-ray emission from the
region between the host galaxy and the radio hot spot, which is usually
ascribed to the inverse Compton emission off the cosmic microwave background
radiation (IC-CMB) (Hardcastle et al., 2002; Croston et al., 2005b; Blundell
et al., 2006). Recent observations reveal that the non-thermal X-ray emission
from the radio lobe increases with red-shift, further supporting the IC-CMB
origin (Gill et al., 2021). Diffuse X-ray emission has also been reported in
the jets of the FR-I class of radio galaxies and has been ascribed to a
distributed particle acceleration mechanism (Hardcastle et al., 2007; Worrall
et al., 2008; Worrall, 2009). An IC-CMB model is also sometimes invoked to
explain X-ray emission from the jets of FR-II radio galaxies and quasars,
however such models require the jet to be highly relativistic and well aligned
with the line of sight and consequently tend to imply very large physical jet
lengths, sometimes in excess of several mega parsecs (Tavecchio et al., 2000;
Celotti et al., 2001; Ghisellini, G. et al., 2005). Further, recent
polarimetric studies and high-energy gamma-ray constraints provide evidence
supporting the synchrotron emission model as the origin of diffuse X-ray
emission from AGN jets (see Perlman et al., 2020, for a recent review). This
consequently requires particles with very high energies to be present in the
jet and also favours a distributed particle acceleration mechanism due to the
short synchrotron lifetime of the radiating particles.

The present work explores, for the first time, the interplay of vital particle
acceleration mechanisms in a weakly magnetised plasma environment such as the
radio lobes of FR-II radio galaxies and studies their effect on the emission
properties of these systems. Due to the complicated evolution of the dynamical
quantities as a result of non-linear plasma flow pattern inside these lobes,
we adopt a numerical approach for this work. In particular, we have employed
MHD simulations to produce radio lobes and analyse the emission features
caused by particle energization in the presence of shocks and underlying
turbulence. We have adopted our recently developed second-order accurate STA
framework (Kundu et al., 2021) for this purpose. Owing to the increased
computational complexity of the developed framework, this paper will focus on
a 2D axisymmetric MHD jet model only while leaving the more computationally
expensive 3D case to forthcoming works.

The paper is organised in the following way: we describe our numerical setup
for simulating a 2D axisymmetric AGN jet in section 2.1. Section 2.2 describes
the numerical model to compute the emission properties. In section 3, we
present the results of the simulations. In section 4 we summarise our findings
and discuss the limitations of our model.

##  2 Numerical Setup

In this section, we describe the numerical setup adopted for the present work.
The radio lobes are typically associated with the termination point of the AGN
jet, where the velocity of the jet material reduces considerably such that
relativistic effects become negligible (Huarte-Espinosa et al., 2011).
Further, as shown by Hardcastle & Krause (2013), numerical simulations of
realistic radio lobe require high Mach number flows as well as very high-
resolution meshes in order to have radio lobes in pressure equilibrium with
the surrounding medium and to resolve the transverse radial equilibrium.
Therefore, to investigate the emission profile of the radio lobes, we focus on
a non-relativistic scenario and perform a two dimensional axisymmetric ideal
MHD simulation using PLUTO code (Mignone et al., 2007). In particular, we
solve the following set of conservation equations,

$\frac{\partial\rho}{\partial t}+\nabla\cdot(\rho\boldsymbol{v})=0\,,$ (1)
$\frac{\partial\boldsymbol{v}}{\partial
t}+(\boldsymbol{v}\cdot\nabla)\boldsymbol{v}=-\frac{1}{\rho}\nabla
P+\frac{1}{\rho}(\nabla\times\boldsymbol{B})\times\boldsymbol{B}\,,$ (2)
$\frac{\partial P}{\partial t}+\boldsymbol{v}\cdot\nabla P+\Gamma
P\nabla\cdot\boldsymbol{v}=0\,,$ (3) $\frac{\partial\boldsymbol{B}}{\partial
t}=\nabla\times(\boldsymbol{v}\times\boldsymbol{B})\,,$ (4)

The quantities $\rho$, $P$, $\boldsymbol{v}$ and $\boldsymbol{B}$ represents
density, pressure, velocity and magnetic field respectively. The magnetic
field $\boldsymbol{B}$ further satisfies the constraint
$\nabla\cdot\boldsymbol{B}=0$. $\Gamma$ represents the ratio of specific heats
and its value is taken to be $5/3$, which is typically considered for
supersonic non-relativistic jets (Massaglia et al., 2016). Eqs. (1)-(4) is
solved with HLLC Riemann solver using piece-wise linear reconstruction, van
Leer flux limiter (Van Leer, 1977) and second order Runge-Kutta time-stepping.
Additionally, we consider divergence cleaning (Dedner et al., 2002) to satisfy
the solenoidal constraint of magnetic field.

###  2.1 Dynamical Setup

The two-dimensional axisymmetric simulations are carried out in a cylindrical
geometry $\\{r,z\\}$ such that the radial and vertical extends range from
$\\{0,0\\}$ to $\\{65L_{0},195L_{0}\\}$ with a resolution of $780\times 2340$.
The physical quantities defined in our simulations are appropriately scaled by
defining length, velocity and density scales. For the length, we define the
jet radius $r_{j}=L_{0}=2$ kpc as the scale length. The core density is
adopted as the scale for density such that $\rho_{0}=5\times 10^{-26}$ gm/cc.
Finally, for an ambient temperature $T_{a}=2$ keV, we define the sound speed
$c_{a}=v_{0}=730$ km/s as the scale velocity.

The ambient medium density is initialized with an isothermal King profile
(King, 1972),

$\rho_{a}=\frac{\rho_{0}}{\left(1+\left(\frac{R}{R_{c}}\right)^{2}\right)^{\frac{3\beta}{2}}},$
(5)

where $\rho_{a}$ is the ambient density that consists of a core with radius
$R_{c}=40L_{0}$ and $R/L_{0}=\sqrt{r^{2}+z^{2}}$ is the spherical radius. The
value of power law index is kept constant as $\beta=0.35$. Initially, the
ambient medium is set in a hydrostatic equilibrium using a gravitational
potential ($\Phi_{\rm k}$) (Krause, 2005),

$\Phi_{\rm k}=\frac{3\beta k_{B}T_{a}}{2\mu
m_{H}}\log\left(1+\left(\frac{R}{R_{c}}\right)^{\mathbf{2}}\right),$ (6)

where $k_{B}$, $\mu$ and $m_{H}$ are the Boltzmann constant, mean molecular
weight and mass of hydrogen atom respectively. The ambient pressure ($P_{a}$)
is computed as follows,

$P_{a}=\frac{\rho_{a}T_{a}k_{B}}{\mu}.$ (7)

The ambient medium is set to be non-magnetized initially with the expectation
that the magnetic field in the environment will have minimal impact on the
non-thermal particle transport and the subsequent emission features within the
lobe.

An under-dense beam of density $\rho_{j}=\eta\rho_{0}$ with velocity $v_{j}$
is continuously injected in the medium from a circular nozzle of radius
$r_{j}$, along the vertical direction ($\hat{z}$) at $t=0$, with $\eta=0.1$
being the density contrast. The nozzle is placed within the numerical domain
with a height of $0.5L_{0}$. The adopted resolution samples the jet nozzle
radius with 12 computational cells. The injection velocity ($v_{j}$) is
obtained by choosing the sonic Mach number $M$ such that,

$v_{j}=Mc_{a},$ (8)

with $M=25.0$. The injected beam includes a toroidal magnetic field ($B_{j}$)
with the following radial profile (Lind et al., 1989)

$B_{j,\phi}=\left\\{\begin{array}[]{cl}B_{m}\frac{r}{r_{m}}&\mbox{$\text{for}\,\,r\leq
r_{m}$}\\\ B_{m}\frac{r_{m}}{r}&\mbox{$\text{for}\,\,r_{m}\leq r\leq
r_{j}$}\\\ 0&\mbox{$\text{otherwise}$},\end{array}\right.$ (9)

where the value of $B_{m}$ is governed by the plasma-beta parameter and
$r_{m}$ is the magnetization radius. Such a magnetic field profile corresponds
to a uniform current density within the radius $r_{m}$, zero current density
between $r_{m}$ and $r$ and a return current at $r$. Further, such a
configuration also respects the symmetry condition on the $z$ axis ($B_{j}=0$
at $r=0$) (Komissarov et al., 2007). Additionally, a suitable gas pressure is
provided inside the jet to ensure radial balance between the hoop stress and
pressure gradient force,

$P_{j}=\left\\{\begin{array}[]{cl}\left(\delta+\frac{2}{\kappa}\left(1-\frac{r^{2}}{r_{m}^{2}}\right)\right)P_{e}&\mbox{$\text{for}\,\,r<r_{m}$}\\\
\delta P_{e}&\mbox{$\text{for}\,\,r_{m}\leq r<r_{j}$}\\\
P_{e}&\mbox{$\text{at}\,\,r=r_{j}$},\end{array}\right.$ (10)

where $\delta=1-\frac{r_{m}^{2}}{\kappa r_{j}^{2}}$ and
$\kappa=\frac{2P_{e}}{B_{m}^{2}}$ and $P_{e}$ is the pressure in units of
$\rho_{0}v_{0}^{2}$ at the nozzle radius, computed from the ambient medium
($P_{e}=P_{a}$ at $r=r_{j}$). Owing to the constraint imposed by 2D
axisymmetric geometry, the induction equation (Eq. 4) does not enable
conversion of the toroidal magnetic field ($B_{\phi}$) to a poloidal one. As a
result, we consider a minimal value of $B_{m}\sim 100$ $\mu$G, to avoid
significant amplification of the $B_{\phi}$ due to its continuous injection
into the computational domain over time. Further, the initial kinetic power of
the jet is calculated from the quantities defined at the jet nozzle (Massaglia
et al., 2016),

$W=\frac{\pi}{2}\left(\frac{\Gamma
k_{B}N_{A}}{\mu}\right)^{\frac{3}{2}}\eta\rho_{0}r_{j}^{2}M^{3}T_{a}^{\frac{3}{2}},$
(11)

where $N_{A}$ is Avogadro’s number. For the choices adopted in the present
work, we obtain $W\simeq 10^{45}$ erg/s corresponding to the FR-II class of
radio galaxies (Fanaroff & Riley, 1974).

For the boundaries, we employ axisymmetric boundary condition about the axis
for the inner $r$ boundary and free flow boundary conditions for all other
boundaries in the computational domain.

###  2.2 Numerical setup to compute emission

The non-thermal emission from the radio lobe is modelled using the Eulerian-
Lagrangian hybrid framework of the PLUTO code (Vaidya et al., 2018; Mukherjee
et al., 2021). It employs passive Lagrangian (or macro-) particles whose
dynamics is governed by the underlying fluid motion. Physically, these macro-
particles represent an ensemble of non-thermal particles (typically leptons)
residing very closely together in physical space with a finite energy
distribution.

The energy distribution of these macro-particles is evolved by solving the
following transport equation,

$\displaystyle\frac{\partial\chi_{p}}{\partial\tau}+\frac{\partial}{\partial\gamma}\left[(S+D_{A})\chi_{p}\right]=\frac{\partial}{\partial\gamma}\left(D\frac{\partial\chi_{p}}{\partial\gamma}\right)$
(12)

where $\tau$ is the proper time, $\gamma\approx p/m_{0}c$ is the Lorentz
factor of the electrons, with $m_{0}$ being the rest mass of the electron and
$c$ is the speed of light in vacuum. The dimensionless quantity
$\chi_{p}=N/n$, with $N(p,\tau)$ being the number density of the non-thermal
particles with momentum between $p$ and $p+dp$ and $n$ being the number
density of the fluid at the position of the macro-particle. The quantity $S$
represents various radiative and adiabatic losses. The acceleration due to
Fermi II order mechanism is given as $D_{A}=2D/\gamma$ with $D$ being the
momentum diffusion coefficient. We, for simplicity, neglect the source and
sink terms in the transport equation.

Eq. (12) is solved using a $2^{nd}$ order accurate finite-volume conservative
implicit-explicit (IMEX) scheme (Kundu et al., 2021). The radiative losses
considered include synchrotron, IC-CMB and adiabatic expansion to model the
cooling processes of relativistic electrons. Additionally, as the particle
spectra in the high energy region falls off rapidly due to various cooling
processes, we follow Winner et al. (2019) and set the values of $\chi_{p}=0$
beyond a threshold $\chi_{\rm cut}=10^{-21}$ . Note that Eq. (12) does not
include shock acceleration; instead, a separate sub-grid prescription is
employed to account for DSA (Vaidya et al., 2018; Mukherjee et al., 2021).

The micro-physics of turbulent acceleration is encapsulated in the diffusion
coefficient $D$. Typically, the empirical form of $D$ is given as an input in
numerical simulations (Donnert & Brunetti, 2014; Vazza et al., 2021), as its
quantification from first principles is complex particularly when applied to
study large scale astrophysical environments. In this work we opt for a
phenomenologically motivated ansatz of exponentially decaying hard-sphere
turbulence as a model of STA inside the radio lobe. We consider the
acceleration timescale ($t_{A}$) as follows (Kundu & Vaidya, 2022),

$t_{A}=\tau_{A}\exp\\{(t-\tau_{t})/\tau_{d}\\}$ (13)

where $\tau_{d}$ is the turbulence decay timescale, $\tau_{A}$ represents the
acceleration timescale when turbulence decay is absent (or
$\tau_{d}\to\infty$) and $t$ is the simulation time. $\tau_{t}$ is the
injection time of the macro-particle in a turbulent region. For a macro-
particle that encounters a shock, its value is set to the time at which the
last shock is encountered. While for those macro-particles that never undergo
any shock, the value of $\tau_{t}$ is set to initial injection time in the
computational domain.

This acceleration timescale has the capability to mimic the decay of
turbulence, generally observed in various astrophysical sources. The decay is
a consequence of the finite lifetime of the turbulence and prevents particles
from being continuously accelerated. For this work, we model $\tau_{A}$ and
$\tau_{d}$ as follows,

$\displaystyle\begin{split}\tau_{A}&=\frac{\tau_{c}(\gamma_{\rm
max}\rightarrow\gamma_{\rm min})}{\alpha},\\\ \tau_{d}&=\tau_{A}.\end{split}$
(14)

where $\tau_{c}(\gamma_{\rm max}\rightarrow\gamma_{\rm min})$ represents the
radiative loss time for a particle to cool from $\gamma_{\rm max}$ to
$\gamma_{\rm min}$ and $\alpha$ is the ratio of synchrotron cooling time and
acceleration time. It also controls the efficiency of STA. Larger $\alpha$
value corresponds to smaller $\tau_{A}$, STA timescale, and $\tau_{d}$,
turbulence damping timescale. Hence larger $\alpha$ indicates faster
stochastic acceleration and faster damping. Also, with lower values of
$\alpha$, the effect of STA asymptotically diminishes. It is a parametric
representation that models the turbulence that actually occurs in realistic
radio lobes of FR-II radio galaxies, which is unresolved in our simulation. In
this work, we vary its value and study how this affects the emission
signatures. The diffusion coefficient can subsequently be written as,

$D=\frac{\gamma^{2}\exp\\{-(t-\tau_{t})/\tau_{d}\\}}{\tau_{A}}.$ (15)

The $\gamma^{2}$ dependency of the diffusion coefficient is a characteristic
of the hard-sphere turbulence. Further, instead of such $\gamma^{2}$ dependent
diffusion coefficient, one could also explore alternative diffusion models.
For example, adopting Bohm type diffusion ($\propto\gamma$) could influence
the results, however such a study of varying dependence of diffusion
coefficient on $\gamma$ is beyond the scope of this paper. To explore the
ramifications of STA with varying efficiency on the emission of the simulated
radio lobe structure, we use two alternative values for $\alpha=10^{4}$ and
$10^{5}$ in this study. Further, to sample the jet cocoon uniformly, we inject
enough ($\sim 20$) macro-particles at every time step in the computational
domain. Initially, the normalized particle spectrum for each macro-particle is
assumed to be a power-law, defined as $\chi_{p}(\gamma)=\chi_{0}\gamma^{-9}$,
ranging from $\gamma_{\rm min}=1$ to $\gamma_{\rm max}=10^{5}$. The value of
$\chi_{0}$ is set by prescribing the energy density of the macro-particles to
be a fraction ($\approx 10^{-4}$) of the initial magnetic energy density. Note
that the initial spectral index has a negligible effect on the emission of the
system at later times as long as we consider a steep power-law.

To compute the emissivity, we convolve the instantaneous energy spectrum of
each macro-particle with the corresponding single particle radiative power and
extrapolate it to the nearest grid cells. In particular we solve the following
integral to compute the emissivity,

$j(\nu^{\prime},n^{\prime},\tau)=\int_{1}^{\infty}\mathcal{P}(\nu^{\prime},\gamma^{\prime},\psi^{\prime})N^{\prime}(\gamma^{\prime},\tau)d\gamma^{\prime}d\Omega^{\prime},$
(16)

where, $\mathcal{P}(\nu^{\prime},\gamma^{\prime},\psi^{\prime})$ is the power
emitted by a non-thermal particle per unit frequency ($\nu^{\prime}$) and unit
solid angle ($\Omega^{\prime}$) with lorentz factor $\gamma^{\prime}$, and
whose velocity makes an angle of $\psi^{\prime}$ with the direction
$n^{\prime}$. $N^{\prime}(\gamma^{\prime},\tau)$ is the number of micro-
particles between lorentz factor $\gamma^{\prime}$ and
$\gamma^{\prime}+d\gamma^{\prime}$ at time $\tau^{\prime}$. In the case of an
axisymmetric simulation, the magnetic field becomes independent of the polar
angle and therefore to consider the line of sight (LOS) effect in the
synchrotron emissivity, appropriate co-ordinate transformation is required
(Meyer et al., 2021). We transform the magnetic field from cylindrical to
Cartesian coordinates and compute the LOS effect by rotating the simulated
structure explicitly. The entire rotation (of $360^{\circ})$ is performed with
an interval of $5^{\circ}$. Subsequently, the intensity maps of the structure
are computed by doing a LOS integration of the calculated emissivity. Note
that all the emissivity calculation is performed by considering a viewing
angle of $\theta=90^{\circ}$ (i.e., along the $z=0$ plane in Cartesian co-
ordinates).

##  3 Results

Figure 1: Normalized density $\rho/\rho_{0}$ evolution of the simulated radio
lobe structure. The images depict a slice through the mid-plane of the
notional 3D volume; all images are reflection-symmetric around the jet axis,
and $z=0$ plane since the simulations are axisymmetric. The color-bar shows a
logarithmic scale of density. Figure 2: Temperature, pressure, absolute
velocity and plasma-beta maps of the simulated jet structure for time $t=117$
Myr. Temperature and pressure are shown in physical units; velocity is shown
in the units of $c$ and the color-bars are shown in logarithmic scale. The
average temperature of the radio lobe is of the order of $\sim 70$ keV,
average plasma-beta is $\sim 32$ and average velocity is $\sim 0.02c$. Figure
3: Evolution of the energy spectrum for a randomly chosen macro-particle for
all the cases described in Table 1. The macro-particle encountered shock at a
dynamical time of $t=25$ Myr. The color-bar shows how much time has elapsed
since the simulation began. The value of the lower end of the color-bar is set
to the time when the macro-particle encountered the final shock. Figure 4: PDF
of the cut-off energy for the entire macro-particle population. Left, Middle
and Right panel shows the PDF for case (a), (b) and (c) respectively.

We categorize the major results from our simulations in two parts. The first
part gives an overview of dynamical aspects of radio lobes and the second part
provides a detailed analysis of multi-wavelength emission signatures and
particle acceleration processes within these lobes.

###  3.1 Dynamics

We have carried out axisymmetric MHD simulations following the initial
conditions described in Section 2 using the relevant jet and ambient medium
parameters. The simulation is carried out up to a physical time of $\sim 120$
Myr. In Fig. 1, we show the density evolution of the injected jets at
different times, viz., t = $37$, $64$, $91$ and $117$ Myr. The density
structure at every time snapshot shows an expanding bi-directional under-dense
region which at a later time ($t=117$ Myr) can be identified as lobes (English
et al., 2016). Similar to Hardcastle & Krause (2013), we notice the formation
of a long, thin lobe initially and a transverse expansion afterward. This
subsequent expansion in the transverse direction is attributed to the
thermalization of the jet material by the shocks present in the lobe. Further,
we observe the formation of vortices at the lobe boundary, which are typically
attributed to Kelvin-Helmholtz instabilities originated from the velocity
shear between the lobe material and shocked ambient material. Moreover, the
entire structure is encapsulated within a forward-moving shock that can be
seen to propagate through the ambient medium. This shock remains in the
computational domain throughout the simulation time, preventing any mass,
energy, and momentum from escaping the domain.

In Fig. 2, we show the temperature (left panel), thermal pressure (second
panel), absolute velocity $|\boldsymbol{v}|$ (third panel) and plasma-beta
(right panel) maps of the bi-directional jet at time $t=117$ Myr. The
temperature of the lobe (average value of $\sim 70$ keV) is relatively higher
than the ambient medium ($T_{a}=2$ keV). This is expected given the presence
of a strong shock at the jet termination region, which is responsible for
heating the jet material in the cocoon. The existence of the strong shock can
further be seen from the pressure map as shown in the second panel of the
figure. The pressure map also provides evidence of multiple re-collimation
shocks along the jet axis. Such shocks are expected to be favourable sites for
accelerating particles via shock acceleration and are known to be a source of
localised high-energy emissions. Further, we observe the velocity of the jet
is within the non-relativistic limit with an average value of $\sim 0.02c$.
The plasma-beta map, as depicted in the right panel of the figure, shows that
the lobes are thermally dominated with an average lobe plasma-beta value of
$\sim 32$.

The under-dense lobes observed in 2D simulations resemble the radio galaxies
in a more consistent manner at later times (Hardcastle & Krause, 2013). In
particular, when the expansion results in the length of the under-dense region
being comparable to the core radius of the galaxy.

Therefore, in this work, for the emission studies, we adopt the dynamical
results at time $t=117$ Myr.

###  3.2 Emission

We now turn our attention to the emission signatures of our model. The
discussion will be based on the comparison of synthetic emission signatures
from different runs considered in our study. The parametric study focuses
mainly on the properties of the stochastic turbulent acceleration mechanism.
The details of these simulation runs are listed in Table 1, which considers
various acceleration scenarios corresponding to different turbulent
acceleration time scales $t_{A}$, while the background thermal fluid evolution
remains exactly the same.

Run ID | DSA | STA | Turbulent Decay | $\alpha$ |  Remarks
---|---|---|---|---|---
Case (a) | YES | NO | NO | 0 |  Energy spectrum exhibits power law with exponential cut-off; PDF of $\gamma_{\rm avg}$ shows power law; SED shows transient peaks.
Case (b) | YES | YES | YES | $10^{4}$ |  Individual macro-particle energy spectrum exhibits curvature; $\gamma_{\rm max}$ PDF indicates accumulation of particles around $10^{4}$; $\gamma_{\rm avg}$ PDF exhibits low-energy cut-off. Peak radiation from synthetic SED is $10^{10}$ Hz through synchrotron and $10^{19}$ Hz via IC-CMB.
Case (c) | YES | YES | YES | $10^{5}$ |  Individual macro-particle energy spectrum exhibits curvature. $\gamma_{\max}$ PDF shows particle accumulation around $10^{5}$; $\gamma_{\rm avg}$ PDF provides evidence of low-energy cut-off. Synthetic SED peak at $10^{13}$ Hz through syncrotron and $10^{21}$ Hz via IC-CMB.
Case (d) | YES | YES | NO | $10^{4}$ |  Individual macro-particle energy spectrum exhibits steady ultra-relativistic Maxwellian structure peaking at $\gamma\approx 10^{4}$.
Case (e) | YES | YES | NO | $10^{5}$ |  Individual macro-particle energy spectrum exhibits steady ultra-relativistic Maxwellian structure peaking at $\gamma\approx 10^{5}$.
Table 1: Properties of the different cases, considered in the present study
for calculating emission from the radio lobe. The first column contains labels
for cases for further reference. The second, third, and fourth columns
represent the presence or absence of DSA, STA, and turbulent decay effects on
the emission runs. The fifth column depicts the value of free parameter
$\alpha$ (Eq.14) chosen for different runs and the last column describes the
hey results for each of the cases.

The results obtained from cases (a) and (b) will be useful to comprehend the
impact of STA and its interplay with DSA. Cases (b) and (c) will highlight the
implications of having different turbulent decay timescales (see Eqs. 13, 14).
For cases (d) and (e), the turbulent decay is turned off by setting
$\tau_{d}\rightarrow\infty$ in Eq. (13). Comparing results from these cases
will demonstrate the effect of the turbulent decay process in our simulations.
In realistic astrophysical environments, we expect the turbulence to decay on
some time scale that is governed by the micro-physical properties of the wave-
particle interaction in that system. As the current work incorporates
turbulence via a sub-grid model, we have explored the implications of
different parameters through these 5 cases. All the results presented in this
section are for a dynamical time of $117$ Myr, unless specified otherwise.
Further, logarithmic binning has been adopted for all the histograms.

####  3.2.1 Effect of Turbulent acceleration on individual macro-particle
energy spectrum

In Fig. 3, we show the evolution of the energy spectra for all the cases
listed in Table 1 for a randomly chosen macro-particle which encountered final
shock at a dynamical time $t=25$ Myr. In the simulations presented in this
work, the majority of the Lagrangian macro-particles are observed to encounter
more than one shock. Among them, we selected this particular particle, which
had experienced multiple shocks only at earlier times, as a representative
candidate to demonstrate the effects of turbulent acceleration on the particle
energy distribution in the downstream of the shock for all the case scenarios.
The effect of multiple shocks on the energy spectrum of a Lagrangian macro-
particle without STA has already been investigated in the context of AGN jet
simulation (see Mukherjee et al., 2021; Giri et al., 2022, for example).

The spectral evolution of the macro-particle of case (a) is shown in the top
left panel. The spectrum exhibits a power-law with a high-energy cut-off which
gradually shifts to lower energy with time, owing to various energy losses.
Additionally, a small hump can be seen in the low-energy part of the spectrum
which is due to an excess of lower energy electrons arising from their higher
energy counterparts due to radiative cooling.

The shape of the spectrum changes considerably when STA is considered in
addition to DSA. For cases (d) and (e) (shown in the right plot of the middle
panel and left plot of the bottom panel respectively), the spectrum exhibits
an ultra-relativistic Maxwellian distribution at later times. This is a
consequence of a steady competition between stochastic acceleration and
radiative losses resulting in the acceleration of low-energy electrons towards
higher energies (Kundu et al., 2021). Moreover, the peak of the distribution
corresponds to the value of $\gamma$ at which acceleration and loss time
scales match, i.e., $\tau_{c}=t_{A}$. We find that the peak corresponds to
$\gamma\approx\alpha$ and it depends on the choice of the turbulent
acceleration timescale (see Eq. 13).

When turbulent decay is included (cases b and c) we observe flatness of the
spectrum in the lower energy regime, as compared to the power-law behavior
observed in case (a), along with a high energy cut-off. The flattening of the
lower energy component of the spectrum is a consequence of the fact that STA
provides a continuous acceleration to all the micro-particles, resulting in
their acceleration to higher energies, depopulating the low energy regime.

Further, we point out that for the macro-particles which have encountered a
shock, STA starts acting in the downstream and it modifies the energy spectra
on a time-scale that depends on $t_{A}$ (Eq. 13) which in turns is regulated
by the turbulent decay time-scale $\tau_{d}$ and consequently develops a cut-
off that moves towards lower energies.

In summary, the spectral evolution of a macro-particle, presented in Fig. 3
for different cases, clearly indicate that the presence of turbulent
acceleration significantly affects the spectral energy distribution and its
evolution. Our results indicate that, in the absence of turbulent decay,
spectral evolution eventually relaxes towards a steady-state configuration in
which energy losses are balanced by turbulent acceleration, while, on
accounting for the decay of turbulence, the energy spectrum exhibits a non-
stationary behaviour in time and the cut-off is governed by the radiative loss
time scale subsequent to the decay of turbulence. Further, the spectrum shows
flattening in the lower energy regime owing to the energization of low energy
micro-particles to higher energy by STA.

####  3.2.2 Effect of Turbulent acceleration on particle population

This section focuses on the effects of turbulent acceleration on the entire
macro-particle population in the lobe. In particular, we compute the effect of
the STA with turbulence decay on the cut-off energy ($\gamma_{\rm max}$) for
the macro-particle population. To compute the cut-off energy of a macro-
particle we consider a generic form of its energy spectrum,

$\gamma^{-m}\exp\left(-\frac{\gamma}{\gamma_{\rm max}}\right),$ (17)

where $m$ can be positive or negative depending on the macro-particle and
$\gamma_{\max}$ is the cut-off energy. The exponential decay term takes care
of the effects on the spectrum due to various radiative losses (see section
3.2.1). The value of $\gamma_{\max}$ is calculated by multiplying Eq. (17) by
a power-law profile, $\gamma^{10}$, and calculating the maximum point of the
resultant curve.

In Fig. 4, we show the probability distribution function (PDF) of the maximum
(or cut-off) energy ($\gamma_{\max}$) attained by individual macro-particles
for cases (a) (left panel), (b) (middle panel), and (c) (right panel). For
case (a), the distribution peaks around $\gamma_{\max}\approx 10^{2}$,
followed by a broken power-law like tail beyond that. The origin of this peak
can be attributed to the presence of various radiative losses in the system.
The peak is also observed to gradually move towards lower values of
$\gamma_{\max}$ with time. To support the above argument, we undertake the
following exercise: for a particle undergoing synchrotron cooling only, the
initial Lorentz factor $\gamma^{{}^{\prime}}$ after a time period of
$t^{{}^{\prime}}$ becomes,

$\gamma^{*}=\frac{1}{C_{0}B^{2}t^{{}^{\prime}}+\frac{1}{\gamma^{{}^{\prime}}}}$
(18)

with $C_{0}=1.28\times 10^{-9}$ is the synchrotron constant for electron, $B$
is the magnetic field. For our case, considering an averaged magnetic field of
$B=19.70\,\mu$G and $t^{{}^{\prime}}=117$ Myr we get $\gamma^{*}\approx
5.4\times 10^{2}$ for a range of $\gamma^{{}^{\prime}}$ values, which
correlates with the position of the peak. The break in the power-law around
$\gamma_{\max}\sim 10^{5}$ is attributed to the continuous injection of the
macro-particles in the computational domain with $\gamma_{\max}=10^{5}$ (see
section 2.2). The presence of an additional smaller peak around
$\gamma_{\max}\sim 10^{9}$ can also be observed. This smaller peak is a
transient feature, which arises from recently shocked macro-particles and is a
manifestation of the continuous injection of jet material along with the
Lagrangian macro-particles inside the computational domain. The presence of
this transient peak has been reported in earlier works as well (see for
example Borse et al., 2021). Further, the power-law trend of the tail of the
PDF is typically ascribed to the interplay between the continuous injection of
macro-particles in the computational domain and the shock acceleration of
these freshly injected particles. Such a power-law like behaviour of the
distribution in an AGN jet cocoon has also been reported in Mukherjee et al.
(2021).

The PDFs for cases (b) and (c) show some additional peaks, as compared to case
(a). The origin of the peak at $\gamma_{\max}\sim 10^{2}$ is similar to case
(a), while the high-energy one ($\gamma_{\max}\sim 10^{9}$) is again due to
recently shocked macro-particles. In addition, one can observe humps at
$\gamma_{\max}\sim 10^{4}$ (for case b) and at $\gamma_{\max}\sim 10^{5}$ (for
case c). Their presence is caused by particles undergoing turbulent
acceleration downstream of the shock, resulting in freezing the evolution of
the cut-off at $\gamma_{\max}\approx\alpha$ for some time due to the
competition between STA and radiative losses and afterwards, due to the decay
of turbulence, the cut-off continues to decrease towards lower energy as
dictated by loss processes.

Figure 5: PDF of the $\gamma_{\rm avg}$ (see Eq. 19) for the entire macro-
particle population. Left, Middle  and Right panel shows the PDF for case (a),
(b) and (c) respectively. Figure 6: Integrated spectrum of the entire macro-
particle population for all the three cases. The portion of the spectrum
highlighted in orange corresponds to the low energy break. The highlighted
portion of the spectrum by blue (and green) corresponds to the high energy
cut-off for case (b) (and case (c)).

To understand the distribution of electron energy within macro-particles, we
also estimate the average value of $\gamma$ (at the final simulation time,
$t=117$ Myr) denoted by $\gamma_{\rm avg}$ as:

$\gamma_{\rm avg}(t)=\frac{\int_{\gamma_{min}}^{\gamma_{\rm max}}\gamma
N(\gamma,t)d\gamma}{\int_{\gamma_{min}}^{\gamma_{\rm
max}}N(\gamma,t)d\gamma}\,.$ (19)

where $\gamma_{\max}$ and $\gamma_{\min}$ are given in section 2.2. In Fig. 5,
we plot the PDF of $\gamma_{\rm avg}$ for the entire macro-particle
population. In the left panel of the figure, we show the PDF for $\gamma_{\rm
avg}$ for case (a). The distribution exhibits a power-law tail
($\propto\gamma_{\rm avg}^{-q}$, with $q\approx 2.54$) beyond $\gamma_{\rm
avg}\sim 10^{2}$. For cases (b) and (c) (depicted in the middle and right
panel of the figure), the PDFs exhibit a power-law distribution starting from
$\gamma_{\rm avg}\sim 10^{3}$ with a small hump and an exponential cut-off.
The hump feature arises due to competition between STA and radiative losses
(see above). It is interesting to note that the slope of the power-law for
both the cases (b) and (c) ($q=0.29$, $0.38$ respectively) is flatter than
case (a). This is a consequence of the fact that STA continuously supplies
energy to the macro-particles by accelerating the low energetic micro-
particles to the higher energy, thus compensating for the radiative losses, as
opposed to the case with only DSA. Finally, in presence of both DSA and STA,
the $\gamma_{\rm avg}$ PDFs exhibit a low-energy break around $\gamma_{\rm
avg}\sim 10^{3}$ owing to the fact that STA boosts low-energy particles to
higher energies. This process is absent if only DSA is present, since there is
no selective mechanism to accelerate only the low-energy particles during
shock acceleration (which involve convolution of the entire upstream spectrum
of each macro-particles to downstream (Mukherjee et al., 2021)) and hence
$\gamma_{\rm avg}$ PDF cannot form a low-energy break.

Further, in Fig. 6, we present the integrated particle spectrum considering
the whole macro-particle population for each of the three case scenarios. The
integrated particle spectrum is calculated as follows:

$F(\gamma)=\sum_{i}\frac{\chi^{i}_{p}(\gamma)}{\mathcal{N}_{i}(\gamma)\int\chi^{i}_{p}(\gamma^{\prime})d\gamma^{\prime}}\,,$
(20)

where $i$ corresponds to individual macro-particles inside the computational
domain, $\chi^{i}_{p}(\gamma)$ is the distribution function of the $i^{th}$\-
macro-particle and $\mathcal{N}_{i}(\gamma)$ represents the number of macro-
particles with Lorentz factor $\gamma$. The DSA spectrum (case (a)) is in the
form of a broken power-law with the break at $\gamma\approx 5\times 10^{2}$
(region highlighted with orange color in the figure). Such a behaviour is
expected when computing a resultant distribution comprising all the macro-
particles, where the spectral evolution is mediated by shock acceleration and
radiative losses (Heavens & Meisenheimer, 1987). The position of the break has
a direct correspondence with the peak in the $\gamma_{\max}$ PDF for case (a)
and can be explained by the same reasoning (see Eq. 18). When STA is taken
into account (cases (b) and (c)), the spectrum exhibits an inverse power-law
behaviour for $\gamma\lesssim 4\times 10^{2}$, followed by a low energy break
and a power-law trend with a high energy cut-off highlighted in blue and green
for cases (b) and (c), respectively in the figure. The spectral behaviour in
the region $\gamma\lesssim 4\times 10^{2}$ is a manifestation of the low-
energy flattening in the individual macro-particle spectrum (see section
3.2.1) due to turbulent acceleration. The origin of the low energy break bears
a similar explanation as the case (a). However, for cases where STA is taken
into account, the cut-off is accompanied by piled up micro-particles (see case
(c) of Fig. 3) as opposed to case (a), which is why the break appears more
prominent in cases (b) and (c). The high energy cut-off in the integrated
particle spectrum (at $\gamma\approx 10^{4}$ for case (b) and $\approx 10^{5}$
for case (c)) is governed by the formation of the quasi-stationary cut-off in
the individual macro-particle spectrum due to the interplay of DSA and STA. As
a result, the position of these high energy cut-offs has an exact
correspondence with the peaks observed in Fig. 4 for the cases where STA is
taken into account. The power-law trend beyond $\gamma\gtrsim 10^{6}$ for all
the case scenarios is a consequence of the continuous macro-particle injection
in the computational domain and a fraction of them subsequently getting shock
accelerated.

In summary, turbulent acceleration with exponential decay modifies the macro-
particles’ maximum energy ($\gamma_{\rm max}$) distribution by presenting
additional hump to the PDFs. The location of the humps is closely connected to
the $\gamma$ of individual macro-particles where $\tau_{c}=t_{A}$. The PDF of
$\gamma_{\rm avg}$ for cases (b) and (c) exhibits a power-law trend with an
exponential cut-off and a low energy break. The integrated spectrum with only
DSA exhibits a low energy break, whereas with STA, an additional cut-off at
high energy is also seen.

Figure 7: Histogram of macro-particles with respect to $B_{\rm eq}/B_{\rm
dyn}$ to further study effect of STA on the macro-particle population.
Histograms are normalized and then scaled with the maximum value. Top panel:
Shows the histograms for 3 different cases at 4 different times (shown in
different colors). Top left, middle  and right panel shows the histogram for
case (a), (b) and (c) respectively. Bottom panel: 2D histograms showing
$\tau_{t}$ vs. $B_{\rm eq}/B_{\rm dyn}$ at the final time $t=117$ Myr for 3
cases. Bottom left, middle  and right panel shows the histogram for case (a),
(b) and (c) respectively. The colorbar at the bottom panel shows the number of
macro-particles.

####  3.2.3 Turbulent acceleration as a sustained acceleration process

In this section we examine how STA supports the macro-particles to sustain
their energy from extreme radiative losses. To properly characterize this
behaviour we consider an equivalent magnetic field for each macro-particle and
compare it with the dynamical magnetic field at the position of the macro-
particle. This is computed from the instantaneous single macro-particle energy
distribution as follows,

$\frac{B_{\rm eq}^{2}}{8\pi}=m_{0}c^{2}\int_{\gamma_{min}}^{\gamma_{\rm
max}}\gamma N(\gamma,t)d\gamma,$ (21)

where $B_{\rm eq}$ is the corresponding equivalent magnetic field.

Following Eq. (21), we compute $B_{\rm eq}$ for cases (a), (b), and (c) and
compare it with the corresponding dynamical magnetic field computed at the
local macro-particle position at each instant, $B_{\rm dyn}$. We plot the time
evolution of the histogram of the quantity, $B_{\rm eq}/B_{\rm dyn}$, on a
logarithmic scale, for all three cases in the top panel of Fig. 7, where
orange, blue, green, and black curves in each panel depict the histogram at
times $5$ Myr, $29$ Myr, $58$ Myr, and $117$ Myr, respectively. All the
histograms are normalized so that the maximum peak value is unity.

As shown in the top left panel, for case (a), the histogram gradually shifts
towards a state with, $B_{\rm dyn}\sim B_{\rm eq}$ as time progresses. The
other two cases (b) and (c) exhibit a similar pattern, where one observes a
broadening of the histogram as well. For case (a), the shape of the PDF can be
observed to evolves to a negatively skewed distribution on a logarithmic
scale. To analyze the reason behind this kind of evolution, we show (bottom
left panel) a 2D histogram depicting the value of $\tau_{t}$ with respect to
the magnetic field ratio, which indicates implying that the macro-particle
population with a larger magnetic field ratio has recently been shocked. This
should not be surprising since the shock acceleration energizes particles
thereby increasing $B_{\rm eq}$. The 2D histogram also shows that a relatively
small fraction of macro-particles has magnetic field ratios larger than unity,
due to the absence of any further acceleration process. As a result these
particles undergo strong cooling and quickly lose their energy, hence
featuring an exponential fall in the histogram beyond $B_{\rm eq}\sim B_{\rm
dyn}$ (top left plot of Fig. 7).

On the contrary, for cases (b) and (c) (top middle and right panels), the 1D
histogram evolves to a more extended structure, which closely resembles the
log-normal shape. This extended form of the histograms is ascribed to the
presence of STA which provides a continuous acceleration to the macro-
particles and helps them maintain their energy even in the presence of
radiative cooling. This is further confirmed by observing the corresponding 2D
histograms in the bottom panels (middle and right, respectively). In contrast
to case (a), both figures show more macro-particles in the region $B_{\rm
eq}/B_{\rm dyn}\gtrsim 1$. Also we can infer that even macro-particles that
were shocked earlier (smaller $\tau_{t}$) feature a higher value of $B_{\rm
eq}/B_{\rm dyn}$ as a result of the fact that with STA macro-particles can
sustain their energy for a longer amount of time.

In summary, for all the cases, one observes the distribution gradually evolve
towards a state where $B_{\rm eq}\sim B_{\rm dyn}$. Further, due to the
presence of STA, as compared to only DSA, the histogram manifests a more
extended structure that is evenly spread due to the macro-particles which were
shocked at earlier time but could sustain their energy from radiative losses
because of STA.

Figure 8: Synthetic spectral energy distribution (SED) for case (a) (shown in
green), case (b) (shown in red) and case (c) (shown in blue). The SED due to
synchrotron mechanism is shown in solid lines and the IC-CMB part is shown in
dashed lines. The vertical axis shows the value of $\nu F_{\nu}$ in arbitrary
units. Figure 9: Spatial distribution of the particles responsible for the
peaks in SED. Top panel: shows the position of the particle population with
$\gamma_{\rm max}\sim 10^{4}$ for case (a) (left), case (b) (middle) and case
(c) (right); Bottom panel: shows the position of the particle population with
$\gamma_{\rm max}\sim 10^{5}$ for case (a) (left), case (b) (middle) and case
(c) (right).

####  3.2.4 Synthetic Spectral Energy Distribution of Radio lobe

In Fig. 8, we present the spectral energy distribution (SED) for cases (a)
(green line), (b) (red line) and (c) (blue line). The SED is calculated by
integrating the emissivity (Eq. 16) along the line of sight (Vaidya et al.,
2018) with two different radiation mechanisms: synchrotron (solid lines) and
IC-CMB (dashed lines). The synchrotron SED shows, for case (a), enhanced
emission in the X-ray band with multiple peaks at $\nu\sim 10^{18}$ and
$\nu\sim 10^{21}$ Hz. These peaks are originated from freshly shocked macro-
particles (Borse et al., 2021; Mukherjee et al., 2021). This can be further
verified analytically using the relation between the critical (or cut-off)
frequency ($\nu_{c}$) of synchrotron radiation and the corresponding $\gamma$
(see for example Eqs. (5.80) from Condon & Ransom, 2016),

$\nu_{c}\approx\frac{\gamma^{2}eB}{2\pi m_{e}c},$ (22)

For instance, with an averaged magnetic field of the lobe $B=19.70$ $\mu$G and
$\nu_{c}\sim 10^{21}$ Hz, one obtains a corresponding value for $\gamma\sim
10^{9}$, which is consistent with the peak in the PDF of $\gamma_{\max}$ as
seen in Fig. 4.

For case (b), in addition to similar shock-induced transient signatures, the
synchrotron emission shows a distinct peak in the low energy GHz radio band
($\nu\sim 10^{10}$ Hz). The origin of such a low energy peak is a direct
evidence of turbulent acceleration and corresponds to the hump in the PDF at
$\gamma_{\max}\sim 10^{4}$ (see middle panel of Fig. 4). Likewise, the
synchrotron peak can also be observed for case (c) at a slightly higher
energy, $\nu\sim 10^{13}$ Hz. The macro-particles that are accelerated via STA
and give rise to the peak in PDF around $\gamma_{\rm max}\sim 10^{5}$ (right
panel of Fig. 4) are mainly contributing to the emission at this frequency
band. The macro-particle population that is stochastically accelerated in
cases (b) and (c) is not only responsible for synchroton emission but also
contributes to the distinct peaks in the IC-CMB spectral energy distribution
($\nu\sim 10^{19}$ Hz for case b, $\nu\sim 10^{21}$ Hz for case c). We have
verified that these correspond to the frequency of the photons scattered of a
population of electrons with energy $\gamma_{\max}\sim 10^{4}$ and
$\gamma_{\max}\sim 10^{5}$ for case (b) and (c), respectively. In fact, the
post-scattering frequency of the photons $\nu_{s}$ is related to the electron
energy as follows:

$\nu_{s}\approx\gamma_{\max}^{2}\nu_{0}\,,$ (23)

where $\nu_{0}$ is the frequency at which the cosmic microwave background
(CMB) radiates. Using $\nu_{0}=160$ GHz in Eq. (23), we find that an electron
population at $\gamma_{\max}\sim 10^{4}$ would scatter the CMB photons at a
frequency of $\sim 10^{19}$ Hz. A similar inference can be drawn for the
origin of the IC-CMB peak around $\sim 10^{21}$ Hz for case (c). Additionally,
the peaks in the $\gamma$-ray band ($\nu\sim 10^{27}$ Hz) for all cases
corresponds to the particles with $\gamma_{\max}\sim 10^{9}$.

After observing the SED and identifying the particle populations responsible
for the various peaks, we proceed to show the spatial distributions of these
particle populations in order to understand the resulting emission structure.
In Fig. 9 we show the spatial distribution of the particle populations
responsible for the aforementioned peaks. The top panels depict the particle
distributions with $\gamma_{\max}\sim 10^{4}$ for case (a) (left plot), case
(b) (middle plot) and case (c) (right plot). These particles are correlated to
the peak in the SED caused by IC-CMB at $\nu\sim 10^{19}$ Hz, as explained
earlier in this section. The macro-particles in case (a) can be seen to be
more confined around the shocks in the beam and, to a lesser extent, to the
cocoon region. The reason for this is that, after the shock acceleration, the
macro-particles energy evolution is governed by loss mechanisms only and as a
result, they lose a considerable amount of energy in a short distance. On the
contrary, when turbulent acceleration is included, the particle distribution
corresponding to $\gamma_{\max}\sim 10^{4}$ stretches over a wider area (see
the upper middle and right plot), since macro-particles can be re-accelerated
via turbulence, sustaining high-energy for a longer distance before losing a
substantial portion of their energy. In comparison to case (a), this extended
spatial distribution implies a more diffuse structure of X-ray radiation
attributable to IC-CMB. In the lower panel of Fig. 9 we show the spatial
distribution of the macro-particles with $\gamma_{\rm max}\sim 10^{5}$,
responsible for the peak in the IC-CMB SED at $\sim 10^{21}$ Hz. Similar to
the former scenario, the particle distribution shows an extended morphology
for case (c) as compared to other two cases, for the same reasons discussed
before. Interestingly, the spatial distributions for case (a) (left panel) and
case (b) (middle panel) have a very similar structure. The reason for this can
be further investigated by comparing the $\gamma_{\max}$ histograms for case
(a) and (b) (left and middle panels in Fig. 4), showing a similar behavior
(after the peak at $\gamma_{\max}\sim 10^{4}$ for the latter).

In summary, we showed that, in the presence of stochastic acceleration, the
emission from the radio lobe changes significantly compared to the case where
STA is neglected. With the inclusion of STA, the spatial distribution of the
X-ray emitting particles through IC-CMB exhibit a wider extent (see Fig. 9) as
compared to the only DSA case, indicating a emission structure which is
diffusive.

Figure 10: Spectral index map and spectral index distribution of the radio
lobe for cases (a), (b) and (c). The spectral index maps are drawn considering
two radio frequencies $1.5$ GHz and $15$ GHz. Top left panel: shows the
spectral index map for case (a), Top middle panel: for case (b) and Top right
panel: for case (c). Bottom panel: Spectral index distribution for all three
cases.

####  3.2.5 Spectral Index distribution

In this section, we focus on the effect of STA on the radio frequency regime
($\leq 15$ GHz). With the advent of several high-resolution low-frequency
telescope arrays, it is possible to quantify the distribution of the spectral
index in extended lobes (Alexander & Leahy, 1987; Harwood et al., 2013). In
this regime, the emission from astrophysical systems are dominated by
synchrotron radiation which follows a power-law relation with the frequency,
$I_{\nu}\propto\nu^{-\delta}$, with $\delta$ being the spectral index. In our
simulation we compute the intensity from the macro-particles energy
distribution (see section 2.2) and further calculate the spectral index
$\delta$ using the equation below:

$\delta=\frac{\log(I_{\nu_{2}})-\log(I_{\nu_{1}})}{\log(\nu_{1})-\log(\nu_{2})}$
(24)

In the top panel of Fig. 10 we show the Gaussian-filtered spectral index maps
of the radio lobe considering two frequencies, $\nu_{1}=1.5\,GHz$ and
$\nu_{2}=15\,GHz$ for cases (a), (b), and (c). All the spectral maps show -
from top to bottom - signs of spectral steepening from the outer regions of
the lobe (near the bow-shock) towards the inner part. Such a spatial
distribution can be further analyzed by observing the bottom panel of the
figure, where we plot the vertical distribution of the spectral index value on
the path depicted by the black dashed line shown in the corresponding top
panel, from the inner region of the lobe to the outer region. The spectral
index distribution behaves similarly for all three cases, showing a rapid
increase followed by a softer (or almost constant) one. By analyzing the slope
of this second part, we obtain an average value for case (a) as $-1.01$, while
for cases (b) and (c) it is $-0.80$ and $-0.49$ respectively. This implies
that the radiation spectrum becomes harder as one increases $\alpha$ in the
lobe. The spatial extent of the region with constant spectral index is larger
for case (b) as compared to cases (a) and (c). For case (a), this directly
follows from the absence of any continuous acceleration mechanism other than
shocks, and the ensuing radiative cooling of the macro-particles in the back
flow over a short time scale. In contrast, for cases (b) and (c), STA provides
additional continuous acceleration to the macro-particles. For this reason,
the macro-particles could be able to radiate for a longer amount of time and
the value of the spectral index could be maintained for a longer distance.
Additionally, due to faster turbulence decay, case (c) maintains the spectral
index for a shorter spatial extent as compared to case (b).

Our results have shown that the signature of continuous acceleration of
particles due to stochastic turbulence impact on several observables,
including the spectral index variation along the lobe. We also observed that
while, with increasing $\alpha$ the spectral index value inside the lobe
increases owing to the shorter acceleration timescale, the extent of the
region with constant spectral index decreases due to turbulence decay. We have
discussed the implications of the synthetic measures quantified in section 3
with multi-wavelength observational signatures in the next section.

##  4 Summary and Discussion

In this work we have presented 2D axisymmetric large-scale numerical
simulations of AGN jets using a fluid-particle hybrid approach, in order to
focus on particle acceleration processes and emissions from radio lobes. In
spite of their limitations, and owing to the prohibitive cost of 3D
computations, 2D models still provide fundamental insights on the interplay
between different acceleration mechanisms and their influence on emission
signatures.

Further owing to the multi-scale nature of the system, the underlying
turbulence is considered as a sub-grid manner and its effect on the cosmic ray
transport is modelled with a phenomenologically motivated ansatz for the
turbulent acceleration timescale that can mimic the turbulence decay process
usually observed in various astrophysical sources. By introducing this
timescale, we solve the cosmic ray transport equation to evolve their energy
distribution, accounting for diffusive shock acceleration (DSA), stochastic
turbulent acceleration (STA) as well as for radiative losses (synchrotron and
inverse Compton), as implemented in the PLUTO code by Vaidya et al. (2018). We
explore different scenarios by selectively including or excluding the
aforementioned acceleration mechanisms, and study their effects on the
emission signatures of the radio lobes.

Below, we summarise the primary results from this work.

  * •

We observe significant modification of the energy spectra of macro-particles
when turbulent acceleration is included in addition to DSA as compared to only
shock acceleration case. The interplay of DSA, STA and and turbulent decay
results in features such as flattening of the spectrum in the low-energy
region and a dynamically evolving high energy cut-off. These features produce
curvature in the particle spectrum which can further manifest in the emission
properties of the radio lobe (Duffy & Blundell, 2012).

  * •

The analysis of the maximum attainable energy results in a unimodal PDF with a
broken power-law tail when only shock acceleration is accounted for, while
when both DSA and STA are included, the PDF exhibits a bimodal structure.
Further, the PDF of the average energy ($\gamma_{\rm avg}$) for each macro-
particle shows a power-law profile with an exponential cut-off on inclusion of
STA. These distributions closely resemble the case in which STA is mediated by
continuous particle injection and escape (see Fig. 2b of Katarzyński et al.,
2006). Here particle injection due to shocks act as a source while the escape
is due to turbulence decay. The lobe integrated spectrum exhibits a broken
power-law like structure for DSA, whereas with STA it displays a high energy
cut-off in addition to the low energy break. The position of the low energy
break corresponds to the $\gamma$ where radiative loss time becomes equal to
the dynamical time. The integrated spectrum generated by including STA can be
utilised as a consistent input for one-zone radio lobe modelling that accounts
for particle acceleration due to turbulence.

  * •

Further analysis of STA and its effect on sustaining the particle’s energy
against radiative cooling is performed through the evolution of $B_{\rm
eq}/B_{\rm dyn}$ histograms, showing for all the three cases, that the system
evolves to a state where $B_{\rm eq}\sim B_{\rm dyn}$. However, with STA, the
corresponding distributions become wider when compared to the only shock
scenario, as a result of the additional energization.

  * •

The study of the synthetic SED of the simulated source demonstrates the
existence of additional peaks in the radio band due to synchrotron emission
and in the X-ray band through the IC-CMB mechanism when STA is taken into
account. Further analysis of the spatial distribution of the macro-particles
corresponding to these additional peaks implies a more extended and diffuse
emission in the X-ray band owing to the interplay of the two acceleration
mechanisms. The extent of the spatial distribution is further observed to be
modulated by changing the value of $\alpha$ (see Eq. 13). This implies that,
with an appropriate choice of $\alpha$, one might achieve diffuse emission
around localized regions inside the radio lobe (for example diffuse
synchrotron emission around the hotspot of 3C445, see Prieto et al., 2002).

  * •

The radio frequency spectral maps along with the spectral index profile inside
the lobe indicate a harder emission spectrum due to STA as compared to the DSA
case. The spectral index is observed to remain constant over a distance inside
the radio lobe whose length gets modulated with the efficiency of the
turbulent acceleration. The value of the spectral index in this region is
$\sim-0.49$ for case (c), for case (b) it is $\sim-0.8$ and for case (a) it is
$\sim-1.01$. This kind of behaviour has also been found in various
observations of radio lobes (Parma et al., 1999). Radio lobes of parsec-scale
AGN jets have been observed to exhibit similar characteristics (Hovatta et
al., 2014). However, it should be noted that from observation of radio lobes
there is no evidence of spectral index $\approx-0.5$ or higher. This,
consequently, may impose a limit to the extent and the effectiveness of STA in
the actual radio lobes.

### Present Limitations and Future Extension

The results shown in the present study represent a first step towards a more
realistic description of the complex interaction between the turbulent radio
lobe material and the non-thermal particles, and it is certainly limited by a
number of considerations.

2D axisymmetric models, for instance, are similar to 3D models only in the
case of stable jets and homogeneous media. Time-dependent jet propagation is
known to be prone to 3D instabilities (e.g., Kelvin–Helmholtz and current-
driven modes) that cannot be captured by axisymmetric models (see, e.g.
Mignone et al., 2010; Bodo et al., 2013, 2016). These instabilities are known
to have an effect on the jet emission (Acharya et al., 2021; Borse et al.,
2021) and can induce a range of non-axisymmetric structures, such as filaments
and shocks along jets and in the back-flowing zone (see for example Tregillis
et al., 2001; Matthews et al., 2019). Such non-axisymmetric structures are
known to enhance the turbulence inside the back-flowing region and hence would
strongly influence particle mixing (Jones et al., 1999).

Another potential issue with 2D axisymmetric simulations is that the induction
equation (Eq. 4) does not allow, because of the $\partial_{\phi}=0$ condition,
conversion of toroidal magnetic field ($B_{\phi}$) to poloidal field (Porth,
2013). This leads to the continuous amplification of the injected $B_{\phi}$
component in the computational domain over time, eventually affecting the jet
dynamics. However, for this work we consider a very small $B_{\phi}$ value to
lessen any substantial impact on the dynamics. Nevertheless, 2D computations
still allow our method to be tested with finer grid spacing providing better
resolution across shocked structures. This would be computationally expensive
in the fully 3D case.

Additionally, we also consider an un-magnetized ambient medium in the
expectation that the magnetic field in the ambient medium will have minimal
impact on the non-thermal particle transport within the lobe.

Our simulations describe the interaction between cosmic ray particles and jet
materials although the former behaves essentially as a passive scalar without
backreaction on the fluid. Future extension will consider more exhaustive two-
fluid approaches by also taking into account energy and momentum transfer
between the two components in a self-consistent way (Girichidis et al., 2020;
Ogrodnik et al., 2021). It should be emphasized that the employment of
parameters in our model is an unwanted, albeit necessary, consequence of the
fact that large scale simulations cannot possibly resolve (and therefore
sample) the small-scale turbulence regions. Sub-scale micro-physical processes
(such as turbulent acceleration timescale or MHD turbulence damping rate) must
therefore be encoded through a sub-grid recipe. In this work, in fact, we
consider a one parameter exponentially decaying hard-sphere turbulence as a
model of STA inside the radio lobe, with certain values for the parameter
($\alpha=10^{4},10^{5}$) and compute the emission signatures from the radio
lobes via synchrotron and IC-CMB processes.

Future extensions of this work will hopefully consider fully 3D
investigations, where the impact of non-axisymmetric plasma instabilities may
deeply affect the morphology. Additionally, the sub-grid prescription of
turbulence decay has a crucial role in governing some of the essential
properties of emission.

## References

  * Acharya et al. (2021) Acharya S., Borse N. S., Vaidya B., 2021, Monthly Notices of the Royal Astronomical Society, 506, 1862
  * Alexander & Leahy (1987) Alexander P., Leahy J. P., 1987, MNRAS, 225, 1
  * Araudo et al. (2018) Araudo A. T., Bell A. R., Blundell K. M., 2018, Nuclear and Particle Physics Proceedings, 297-299, 242
  * Asano & Hayashida (2018) Asano K., Hayashida M., 2018, ApJ, 861, 31
  * Belmont et al. (2008) Belmont R., Malzac J., Marcowith A., 2008, A&A, 491, 617
  * Blandford et al. (2019) Blandford R., Meier D., Readhead A., 2019, ARA&A, 57, 467
  * Blundell et al. (2006) Blundell K. M., Fabian A. C., Crawford C. S., Erlund M. C., Celotti A., 2006, ApJ, 644, L13
  * Bodo et al. (2013) Bodo G., Mamatsashvili G., Rossi P., Mignone A., 2013, Monthly Notices of the Royal Astronomical Society, 434, 3030
  * Bodo et al. (2016) Bodo G., Mamatsashvili G., Rossi P., Mignone A., 2016, Monthly Notices of the Royal Astronomical Society, 462, 3031
  * Borse et al. (2021) Borse N., Acharya S., Vaidya B., Mukherjee D., Bodo G., Rossi P., Mignone A., 2021, A&A, 649, A150
  * Brunetti & Lazarian (2007) Brunetti G., Lazarian A., 2007, MNRAS, 378, 245
  * Brunetti et al. (2001) Brunetti G., Bondi M., Comastri A., Pedani M., Varano S., Setti G., Hardcastle M. J., 2001, ApJ, 561, L157
  * Bykov & Fleishman (1992) Bykov A. M., Fleishman G. D., 1992, Monthly Notices of the Royal Astronomical Society, 255, 269
  * Carilli et al. (1991) Carilli C. L., Perley R. A., Dreher J. W., Leahy J. P., 1991, ApJ, 383, 554
  * Celotti et al. (2001) Celotti A., Ghisellini G., Chiaberge M., 2001, MNRAS, 321, L1
  * Condon & Ransom (2016) Condon J. J., Ransom S. M., 2016, Essential Radio Astronomy
  * Croston et al. (2005a) Croston J. H., Hardcastle M. J., Harris D. E., Belsole E., Birkinshaw M., Worrall D. M., 2005a, ApJ, 626, 733
  * Croston et al. (2005b) Croston J. H., Hardcastle M. J., Harris D. E., Belsole E., Birkinshaw M., Worrall D. M., 2005b, The Astrophysical Journal, 626, 733
  * Dedner et al. (2002) Dedner A., Kemm F., Kröner D., Munz C. D., Schnitzer T., Wesenberg M., 2002, Journal of Computational Physics, 175, 645
  * Dermer et al. (1996) Dermer C. D., Miller J. A., Li H., 1996, ApJ, 456, 106
  * Donnert & Brunetti (2014) Donnert J., Brunetti G., 2014, Monthly Notices of the Royal Astronomical Society, 443, 3564
  * Duffy & Blundell (2012) Duffy P., Blundell K. M., 2012, MNRAS, 421, 108
  * English et al. (2016) English W., Hardcastle M. J., Krause M. G. H., 2016, MNRAS, 461, 2025
  * Fan et al. (2008) Fan Z.-H., Liu S., Wang J.-M., Fryer C. L., Li H., 2008, ApJ, 673, L139
  * Fanaroff & Riley (1974) Fanaroff B. L., Riley J. M., 1974, MNRAS, 167, 31P
  * Ferrand & Marcowith (2010) Ferrand G., Marcowith A., 2010, A&A, 510, A101
  * Ghisellini, G. et al. (2005) Ghisellini, G. Tavecchio, F. Chiaberge, M. 2005, A&A, 432, 401
  * Gill et al. (2021) Gill A., Boyce M. M., O’Dea C. P., Baum S. A., Kharb P., Campbell N., Tremblay G. R., Kundu S., 2021, ApJ, 912, 88
  * Giri et al. (2022) Giri G., Vaidya B., Rossi P., Bodo G., Mukherjee D., Mignone A., 2022, A&A, 662, A5
  * Girichidis et al. (2020) Girichidis P., Pfrommer C., Hanasz M., Naab T., 2020, MNRAS, 491, 993
  * Hardcastle & Croston (2020) Hardcastle M. J., Croston J. H., 2020, New A Rev., 88, 101539
  * Hardcastle & Krause (2013) Hardcastle M. J., Krause M. G. H., 2013, MNRAS, 430, 174
  * Hardcastle et al. (2002) Hardcastle M. J., Birkinshaw M., Cameron R. A., Harris D. E., Looney L. W., Worrall D. M., 2002, The Astrophysical Journal, 581, 948
  * Hardcastle et al. (2007) Hardcastle M. J., et al., 2007, ApJ, 670, L81
  * Hardcastle et al. (2009) Hardcastle M. J., Cheung C. C., Feain I. J., Stawarz Ł., 2009, MNRAS, 393, 1041
  * Harris et al. (2019) Harris D. E., et al., 2019, ApJ, 873, 21
  * Harwood et al. (2013) Harwood J. J., Hardcastle M. J., Croston J. H., Goodger J. L., 2013, Monthly Notices of the Royal Astronomical Society, 435, 3353
  * Heavens & Meisenheimer (1987) Heavens A. F., Meisenheimer K., 1987, MNRAS, 225, 335
  * Hovatta et al. (2014) Hovatta T., et al., 2014, AJ, 147, 143
  * Huarte-Espinosa et al. (2011) Huarte-Espinosa M., Krause M., Alexander P., 2011, MNRAS, 417, 382
  * Jones et al. (1999) Jones T. W., Ryu D., Engel A., 1999, The Astrophysical Journal, 512, 105
  * Katarzyński et al. (2006) Katarzyński K., Ghisellini G., Mastichiadis A., Tavecchio F., Maraschi L., 2006, A&A, 453, 47
  * King (1972) King I. R., 1972, ApJ, 174, L123
  * Kirk et al. (1996) Kirk J. G., Duffy P., Gallant Y. A., 1996, A&A, 314, 1010
  * Komissarov et al. (2007) Komissarov S. S., Barkov M. V., Vlahakis N., Königl A., 2007, MNRAS, 380, 51
  * Krause (2005) Krause M., 2005, Astronomy & Astrophysics, 431, 45
  * Kundu & Vaidya (2022) Kundu S., Vaidya B., 2022, Interplay of various particle acceleration processes in astrophysical environment, to appear in ”Proceedings of the International Astronomical Union (IAU)”
  * Kundu et al. (2021) Kundu S., Vaidya B., Mignone A., 2021, The Astrophysical Journal, 921, 74
  * Lind et al. (1989) Lind K. R., Payne D. G., Meier D. L., Blandford R. D., 1989, ApJ, 344, 89
  * Liu et al. (2004) Liu S., Petrosian V., Melia F., 2004, ApJ, 611, L101
  * Mahatma et al. (2019) Mahatma V. H., Hardcastle M. J., Croston J. H., Harwood J., Ineson J., Moldon J., 2019, Monthly Notices of the Royal Astronomical Society, 491, 5015
  * Marcowith & Casse (2010) Marcowith A., Casse F., 2010, A&A, 515, A90
  * Marcowith et al. (2020) Marcowith A., Ferrand G., Grech M., Meliani Z., Plotnikov I., Walder R., 2020, Living Reviews in Computational Astrophysics, 6, 1
  * Massaglia et al. (2016) Massaglia S., Bodo G., Rossi P., Capetti S., Mignone A., 2016, A&A, 596, A12
  * Massaro et al. (2018) Massaro F., et al., 2018, The Astrophysical Journal Supplement Series, 234, 7
  * Matthews et al. (2019) Matthews J. H., Bell A. R., Blundell K. M., Araudo A. T., 2019, Monthly Notices of the Royal Astronomical Society, 482, 4303
  * Meisenheimer (2003) Meisenheimer K., 2003, New Astronomy Reviews, 47, 495
  * Mertsch & Petrosian (2019) Mertsch P., Petrosian V., 2019, A&A, 622, A203
  * Meyer et al. (2021) Meyer D. M. A., Pohl M., Petrov M., Oskinova L., 2021, MNRAS, 502, 5340
  * Mignone et al. (2007) Mignone A., Bodo G., Massaglia S., Matsakos T., Tesileanu O., Zanni C., Ferrari A., 2007, The Astrophysical Journal Supplement Series, 170, 228
  * Mignone et al. (2010) Mignone A., Rossi P., Bodo G., Ferrari A., Massaglia S., 2010, Monthly Notices of the Royal Astronomical Society, 402, 7
  * Mukherjee et al. (2021) Mukherjee D., Bodo G., Rossi P., Mignone A., Vaidya B., 2021, MNRAS, 505, 2267
  * Mullin et al. (2008) Mullin L. M., Riley J. M., Hardcastle M. J., 2008, MNRAS, 390, 595
  * O’Sullivan et al. (2009) O’Sullivan S., Reville B., Taylor A. M., 2009, MNRAS, 400, 248
  * Ogrodnik et al. (2021) Ogrodnik M. A., Hanasz M., Wóltański D., 2021, ApJS, 253, 18
  * Parma et al. (1999) Parma P., Murgia M., Morganti R., Capetti A., de Ruiter H. R., Fanti R., 1999, A&A, 344, 7
  * Perlman et al. (2020) Perlman E. S., Clautice D., Avachat S., Cara M., Sparks W. B., Georganopoulos M., Meyer E., 2020, Galaxies, 8
  * Petrosian (2012) Petrosian V., 2012, Space Sci. Rev., 173, 535
  * Porth (2013) Porth O., 2013, MNRAS, 429, 2482
  * Prieto et al. (2002) Prieto M. A., Brunetti G., Mack K.-H., 2002, Science, 298, 193
  * Schlickeiser & Dermer (2000) Schlickeiser R., Dermer C. D., 2000, A&A, 360, 789
  * Tavecchio et al. (2000) Tavecchio F., Maraschi L., Sambruna R. M., Urry C. M., 2000, The Astrophysical Journal, 544, L23
  * Tavecchio et al. (2022) Tavecchio F., Costa A., Sciaccaluga A., 2022, Monthly Notices of the Royal Astronomical Society: Letters
  * Tregillis et al. (2001) Tregillis I. L., Jones T. W., Ryu D., 2001, ApJ, 557, 475
  * Vaidya et al. (2018) Vaidya B., Mignone A., Bodo G., Rossi P., Massaglia S., 2018, The Astrophysical Journal, 865, 144
  * Van Leer (1977) Van Leer B., 1977, Journal of Computational Physics, 23, 276
  * Vazza et al. (2021) Vazza F., Wittor D., Brunetti G., Brüggen M., 2021, A&A, 653, A23
  * Vurm & Poutanen (2009) Vurm I., Poutanen J., 2009, The Astrophysical Journal, 698, 293
  * Winner et al. (2019) Winner G., Pfrommer C., Girichidis P., Pakmor R., 2019, Monthly Notices of the Royal Astronomical Society, 488, 2235
  * Worrall (2009) Worrall D. M., 2009, A&A Rev., 17, 1
  * Worrall et al. (2008) Worrall D. M., et al., 2008, ApJ, 673, L135
  * de Vries et al. (2018) de Vries M. N., et al., 2018, MNRAS, 478, 4010