Particles in Relativistic MHD Jets. I. Role of Jet Dynamics in Particle Acceleration

Relativistic jets from (supermassive) black holes are typically observed in nonthermal emission, caused by highly relativistic electrons. Here, we study the interrelation between three-dimensional (special) relativistic magnetohydrodynamics, and particle acceleration in these jets. We inject Lagrangian particles into the jet that are accelerated through diffusive shock acceleration and radiate energy via synchrotron and inverse Compton processes. We investigate the impact of different injection nozzles on the jet dynamics, propagation, and the spectral energy distribution of relativistic particles. We consider three different injection nozzles—injecting steady, variable, and precessing jets. These jets evolve with substantially different dynamics, driving different levels of turbulence and shock structures. The steady jet shows a strong, stationary shock feature, resulting from a head-on collision with an inner back-flow along the jet axis—a jet inside a jet. This shock represents a site for highly efficient particle acceleration for electrons up to a few tens of TeV and should be visible in emission as a jet knot. Overall, we find that the total number of shocks is more essential for particle acceleration than the strength of the shocks. The precessing jet is most efficient in accelerating electrons to high energies reaching even few hundred TeVs, with power-law index ranging from 2.3 to 3.1. We compare different outflow components, such as the jet and the entrained material concerning particle acceleration. For the precessing nozzle, the particle acceleration in the entrained material is as efficient as that in the jet stream. This is due to the higher level of turbulence induced by the precession motion.


INTRODUCTION
Astrophysical jets are often observed as relativistic, well-collimated outflows of plasma and fields originating in a region near a compact object accreting matter.These jets are present in a number of astrophysical systems on various scales including Active Galactic Nuclei (AGNs), high mass X-ray binaries or microquasars, gamma-ray bursts (GRBs) and young stellar objects (YSOs).Jets originate from a deep gravitational potential well and from sources hosting strong magnetic fields and an accretion disk (Hawley et al. 2015).
AGN jets travel long distances, reaching hundreds of kpc from the supermassive black hole.The parent galaxies of these sources (called radio galaxies), can be categodubey@mpia.de,fendt@mpia.de,bvaidya@iiti.ac.in * Fellow of the International Max Planck Research School for Astronomy & Cosmic Physics at the University of Heidelberg rized into two classes (Fanaroff & Riley 1974) depending on their radio jet kinetic power and morphology.
Although the processes responsible for formation and launching of these jets are not yet completely understood, the power released from the rotational energy of the black hole (Blandford & Znajek 1977) or the accretion flow (Blandford & Payne 1982) are thought to be the sources of jet energy.This power is ultimately transferred from the gravitational energy to the jet kinetic energy via the magnetic field.The magnetic field can be deduced from observations of power-law spectra, polarization spectral distribution and polarisation.The magnitude of the magnetic field is typically obtained from the cutoff in power-law spectra, whereas the direction is determined via polarisation of these systems, (Meisenheimer et al. 1997;Carilli et al. 1999;Heavens & Meisenheimer 1987;Brunetti et al. 2003) that is due to non-thermal particles that undergo synchrotron and inverse Compton losses.
AGN jets are typically super-(magneto)sonic, and hence produce strong shocks leading to a turbulent magnetic field and velocity component.This makes these sources an ideal site for accelerating particles to high energies (Hargrave & Ryle 1974).Strong shocks, through diffusive shock acceleration, play a vital role in accelerating particles (Krymskii et al. 1978;Blandford & Ostriker 1978;Bell 1978a,b).Additionally, various other mechanisms of particle acceleration are also at play in these sources.These include second-order Fermi acceleration (Fermi 1949;Kundu et al. 2021) particularly in the radio lobes, magnetic reconnection in highly magnetized jet regions (Giannios et al. 2009) and shear acceleration (Rieger et al. 2007;Sironi et al. 2021).
The dynamics of the jet seems to plays an essential role in the emission and subsequently, in studying the shock structure and particle acceleration.In particular, jets that undergo a time-dependent injection of material (and magnetic field), such as intermittent jets and precessing jets, are able to drive a high excess of turbulence or instabilities in the jet flow.This leads to more shocks, and hence, more sites of particle acceleration, and also to an increased efficiency of particle acceleration.(Yates et al. 2018;Giri et al. 2022).
The AGN jet sources are also thought to be one of the sites for production of ultrahigh energy cosmic rays (Pierre Auger Collaboration et al. 2017;Aab et al. 2018;Matthews et al. 2018).Although several studies have pointed out the ineffectiveness of highly relativistic shocks to accelerate particles to such high energies (Reville & Bell 2014;Bell et al. 2018), other investigations have shown that non-relativistic or mildly relativistic shocks can accelerate particles to high energies (Kirk & Duffy 1999;Marcowith et al. 2016;Matthews et al. 2019;Marcowith et al. 2020;Araudo et al. 2021;Ortuño-Macías et al. 2022).Such shocks are capable of accelerating both protons and electrons.In particular, existence of high energy protons is crucial for understanding the origin of neutrino emission associated with blazar jets (Mannheim 1993).Even at large scales, these accelerated particles play a crucial role in governing the heating due to feedback.Simulations of jet inflated bubbles in presence of cosmic ray transport have shown the dominant role of cosmic ray mediated heating at large cluster scales (Ehlert et al. 2018).
As the propagation of astrophysical jets is highly dynamical due to the various non-linear physical processes at work, numerical simulations play an essential role for studying jets.A number of numerical simulations have been performed to study various aspects of jet propagation (for relativistic jets see e.g.Leismann et al. 2005;Mizuno et al. 2007;Keppens et al. 2008;Rossi et al. 2008;Meliani & Keppens 2009;Mignone et al. 2010).In particular, shock diagnostics and particle acceleration in the jet are studied in Mukherjee et al. (2021) for a straight jet with small velocity perturbations to induce turbulence through three-dimensional (3D) relativistic magnetohydrodynamical (RMHD) simulations.
In our paper, we present 3D (special) RMHD simulations of jets injected into a uniform ambient medium.We inject Lagrangian macro-particles into the propagating jet that are advected along with the flow, each of them representing an ensemble of electrons.We then look how these particles are accelerated along the jet flow by internal and external shocks (Vaidya et al. 2018).We study different kinds of time-dependent injection nozzles and investigate how injection from these nozzles affects the dynamics of the jets -its turbulence, shock structure, symmetry -and, as a result, how particle acceleration varies for these nozzles.More specifically, we use three different nozzles -injecting (i) a straight jet, (ii) a precessing jet and (iii) an intermittent jet into the ambient medium.
The different nozzles we apply are thought to be generated by different jet launching conditions at the jet origin.Jet launching itself -the transition from accretion to ejection -is not further considered in our simulations.We may distinguish between jets being launched from resistive accretion disks, as modeled in the non-relativistic case in axisymmetry (see e.g.Casse & Keppens 2002;Zanni et al. 2007;Sheikhnezami et al. 2012;Stepanovs & Fendt 2016).Fully 3D simulations e.g. for binary systems in a Roche potential (Sheikhnezami & Fendt 2015) indicate on precession effects on the accretion disk and thus the jet launched from that disk.
Jet launching from the black hole ergosphere (the BZ process) has been numerically modeled since two decades, mostly applying the evolution from an initial torus surrounding the black hole (McKinney & Gammie 2004;Tchekhovskoy et al. 2011).More recent work applying resistive GR-MHD has been able to consider in addition the evolution of disk winds and jets from thin accretion disks and comparing them to the central spine jet (Qian et al. 2018;Vourellis et al. 2019;Dihingia et al. 2021), even including a disk dynamo (Vourellis & Fendt 2021).We note that the kinematics of the BZ jets obtained in these simulations is biased by the floor model for density and pressure that needs to be considered.
The observation of jet knots hints to a time-dependent mechanism behind their origin, probably a timedependent injection (nozzle).Some recent observations indeed provide indication for precessing AGN jets (see e.g.Britzen et al. 2018;von Fellenberg et al. 2023), and their connection to binary black holes (see e.g.Kharb et al. 2017;Krause et al. 2019).
Our paper is structured as follows.In section 2 we detail our model approach and the numerical methods.In section 3 we discuss the dynamical evolution of the different jets, their energetics and turbulence level.In section 4 we finally present energy spectra of particles moving with the jet and being accelerated by shocks arising in the turbulent jet motion.We conclude with a summary in section 5.
Here, D is the laboratory density, m is the momentum, B is the magnetic field in the lab frame, E t is the total energy density, v is the velocity, γ f is the Lorentz factor of the fluid, and Ī is the diagonal tensor.The magnetic field can be represented by the 4-vector with the magnetic energy density the momentum vector m i = w t γ f v i −b 0 b i , the relativistic total enthalpy w t = ρh + b 2 , and the total pressure of the fluid p t = p + (b 2 /2) where p is the gas pressure (Del Zanna et al. 2003;Mignone & McKinney 2007).Further, the specific enthalpy h is related to the internal energy ε of the gas as where ρ = D/γ f is the density in the fluid frame and p is the gas pressure.
In order to close the aforementioned equations, an equation of state relating h with ρ is further needed (which we discuss in Section 2.1).With this, the total energy density in the lab frame is given as (5) (Mignone & McKinney 2007), where the first two terms represent the kinetic, thermal and the rest mass energy density while the last two terms represent the electromagnetic energy density.Specifically, the third term represents the magnetic energy density, whereas the last term represents the electric energy density as a result of an electric field E = −v×B (Keppens et al. 2008), altogether considering the contribution by Poynting flux.
To study the particle acceleration, we use the Lagrangian particle module of the PLUTO code (Vaidya et al. 2018), which uses passive Lagrangian particles moving in an Eulerian grid to model the non-thermal spectral signatures from relativistic jets.Essentially, we consider particle acceleration due to diffusive shock acceleration (DSA; or first-order Fermi acceleration) following the sub-grid approach developed in Vaidya et al. (2018); Mukherjee et al. (2021).The Lagrangian macroparticle, in this approach, represents an ensemble of particles (electrons) following the fluid, initially distributed as following a power-law energy spectrum.This initial spectrum subsequently evolves for each macro-particle in time and momentum space, taking into account acceleration due to DSA and cooling as a result of adiabatic expansion, synchrotron radiation processes, and inverse Compton (IC) scattering of particle in a background of cosmic microwave background (CMB) photons.For this purpose, the particle code solves the relativistic transport equation for cosmic rays in scattering medium (Webb 1989) which we discuss in Section 2.4.

Equation of State
In order to solve the RMHD equations mentioned above, we need a proper closure provided by an additional equation, that is the equation of state (e.o.s.), relating two thermodynamic parameters (e.g.density and internal energy).
For a perfect gas, this relationship is derived using the relativistic theory (Synge 1957).However, applying this general relationship in a numerical code is too time-consuming and for saving computational resources, a constant-Γ e.o.s. is often used.We apply an approximate e.o.s. which follows a simple analytical form and is fast and suitable to be adopted for numerical studies, (Mathews 1971;Mignone et al. 2005) where Θ = p/ρ is the temperature.This e.o.s., widely denoted as Taub-Mathews (TM) e.o.s., differs from the theoretical Synge e.o.s.only by a few percent, and provides thermo-dynamical variables that are very similar to the constant-Γ e.o.s. in the limit of a cold gas (Γ = 5/3, Θ → 0) and a hot gas (Γ = 4/3, Θ → ∞).For intermediate temperatures, the respective values vary smoothly between the two limiting cases.For the TM e.o.s., the adiabatic index is Γ = (h − 1)/(h − 1 − Θ), whereas the sound speed relative to the fluid c s is defined as Applying c s we define the (ordinary) relativistic sonic Mach number in the lab frame M s = v/c s and the proper sonic Mach number with the proper speed of the fluid u and the proper sound speed u s relative to the fluid (and similarly for the Lorentz factors γ f and γ fs , respectively).Here we follow the arguments by Konigl (1980) and others stating that the proper sonic Mach number M s is the relativistic generalization of the Newtonian Mach number.
Similarly, the Alfvén velocity c a , the magnetisation σ and the plasma β are defined as and similarly for the related Alfvén Mach numbers M a and M a (Mizuno et al. 2014).
As extractable energy E x of the fluid in the lab frame, we define the total energy E t (Equation 5) subtracted by the rest mass energy density Dc 2 .Substituting the specific enthalpy h for the TM e.o.s we obtain where the first and second terms on the right hand side represent the kinetic energy density E kin and the thermal energy density E th , respectively, and the last two terms represent the electromagnetic energy density E em .
The extractable power of the jet P x at time t is calculated by integrating the extractable energy density E x over the volume V and dividing by the time t.Thus, is the (volume-weighted) mean extractable energy density over the outflow volume V (considering cells, where v > 0).We have thoroughly tested both ideal and TM (implemented in PLUTO as Taub e.o.s.) e.o.s. in our simulations in order to compare their effects on the flow dynamics and internal structure.We came to the conclusion that the TM e.o.s. is a much better approach when treating (relativistic) gases of different temperature, such as the low-density, hot gas in the jet, and the high-density, cool gas of the ambient medium.The TM approach is particularly interesting when it comes to the dynamics of the jet cocoon or back-flow material consisting of shocked and entrained gas, thus a mixture of jet gas and ambient gas.

Numerical Specifics
Applying a 3D Cartesian grid, we model magnetized, rotating, one-sided relativistic jets injected from different injection nozzles.
In order to be able to study the turbulence and shock properties properly, it is necessary to resolve small-scale motions as good as possible.We have thus performed a resolution study by modeling and comparing the evolution and structure of the jet at different grid resolution (see Appendix B) applying the same input profiles and parameters values as specified for the steady, straight jet setup in the next section.For the resolution study we have compared simulations with a resolution of 5, 10, 15, 20, 25, and 30 cells per jet radius r j .We find that while the jet dynamical variables as well as particle energy spectrum are sufficiently converged for a resolution of 25 grid cells per jet radius.
We use divergence cleaning to implement the solenoidal condition for the magnetic field ∇•B = 0. We adopt a multidimensional shock flattening algorithm, a second order Runge-Kutta time-stepping, a Courant-Friedrichs-Levy number of 0.25, and a HLL Riemann solver with linear reconstruction that is second order accurate in space.
While the MHD simulations are basically scale-free and could be applied to a variety of sources of different size and energy output, the treatment of radiation as well as of particles requires a proper astrophysical scaling, thus adapted to individual sources of interest.Radiation effects, such as particle acceleration and cooling as well as determining emissivities and intensities, depend on actual density, pressure (temperature) or magnetic field strength.For the present paper, since we are mostly interested in the dynamics and the particle acceleration, we work in dimensionless code units for most of the time (except where we explicitly mention the physical units).In order to convert from code units to astrophysical units, a suitable scaling factor must be applied.However, the synchrotron cooling modeled in our simulations depends on the physical value of the magnetic field.In order to deal with that, we have chosen one set of such normalization units in our simulations (see Table 1).The different variables evolved in our code are normalized based on the choice of three basic Note-Conversion factors from code units to physical scales.Shown for the three basic parameters: length l0, speed v0 and density ρ0 along with the derived normalization factors for time t0, magnetic field B0, pressure p0, temperature T0, and power P0.
unit parameter, the unit length l 0 = r j , the unit speed v 0 = c and the unit density ρ 0 .Other variables or the normalized accordingly as where t 0 , B 0 , p 0 and T 0 represent the unit time, the unit magnetic field, the unit pressure and temperature, respectively, while m u , µ and k B denote the atomic mass unit, the mean molecular weight and the Boltzmann constant, respectively1 .The unit for the power is then defined as P 0 = l 3 0 ρ 0 v 2 0 /t 0 .

Initial Conditions
We employ a constant density profile across the domain except within a cylindrical region inside the domain, with a radius r j = 1, a height z j = 1, centered at x = y = z = 0.This is the injection nozzle.Initially, the ambient medium is defined as the region outside the injection nozzle with a density ρ a = 1000, and being at rest, v a = 0.The initial gas pressure p = 0.1 is constant throughout the domain (including the injection nozzle).
The magnetic field in the ambient medium is purely vertical, B z,a = 0.176, thus B r,a = B φ,a = 0. Here, B r,a , B φ,a and B z,a are the components of the magnetic field in the ambient medium in cylindrical coordinates r, φ and z, respectively.

Boundary Conditions -Steady Jet Equilibrium
We investigate three different injection nozzles -a steady injection, and two nozzles with time-dependent injections giving -a jet with a variable velocity and a precessing jet.Here, we first describe the boundary con-ditions for the steady jet injection.A time variation of these will be applied for the other two nozzles (see next sub-section).
In all our simulations, we inject an under-dense jet with density ρ j = 1 from an injection nozzle of radius r j = 1 and height z j = 1, centered at x = y = z = 0.The nozzle is prescribed by a PLUTO user-defined internal boundary condition.
We have implemented outflow boundary condition at all other boundaries of the domain for all setups, applying a zero-gradient condition for all variables.Through this, we make sure that the material at the boundary can leave the domain freely.Note that using such boundary condition may give rise to an inward Lorentz force in case of sub-fast magnetosonic flows, resulting in artificial collimation of the jet (Porth & Fendt 2010).To avoid this, we restrict the outflow boundaries sufficiently far from the jet stream, not allowing the jet material to leave the domain during the simulation.
Since our main goal is to investigate how the physics and geometry of the nozzle impacts the jet propagation and subsequent high energy particle acceleration, we find it essential to apply a verified equilibrium solution for the gas and the magnetic field that is injected through the nozzle.Any internal instability of the injected gas may lead to an (artificial) modification of the internal equilibrium.
As a proper equilibrium condition for the jet injection we apply the equilibrium solution derived by Bodo et al. (2019) inside the jet injection nozzle.For convenience of the reader, we show the corresponding equations for the leading variables (e.g.velocity, magnetic field) along with plots of the injection profiles in Appendix A. Applying these profiles for our simulations, we obtain a stable jet injection with the bulk motion in the selected injection direction (z-axis) and a rotation around the jet axis.We refer to this simulation setup with steady injection as std throughout the text.

Boundary Conditions -Time Variation
We now describe the time-dependent nozzles for a straight jet with a time-variable injection and a precessing jet.For the time-variable injection, we choose to vary the velocity of the jet in the z−direction, applying a sinusoidal variation above a floor velocity, where v z,flr is a floor velocity (comparable to the other two setups), v 0 is the amplitude of the variable component of v z , and ω var is the frequency of the variation.
Here, we set v z,flr and v 0 at a value of 80 and 20 percent of v z , the velocity of the steady jet in z−direction, respectively.The maximum velocity (in z−direction) of the injected jet in variable nozzle setup is then same as that of the steady nozzle setup, corresponding to a Lorentz factor γ f = 10 at the jet axis.Through this we get a straight jet with a velocity varying sinusoidally with a period P var = 2π/ω var = 5 in the z−direction.
We refer to this simulation setup with variable injection as var throughout the text.
The setup for the precessing nozzle is more complex as we need to change that vector components of velocity and magnetic field periodically, while maintaining the equilibrium solution across the nozzle.We thus need to introduce additional terms in the velocity vector of the rotating jet described in previous section.We do this by using transformation matrices T x , which represents a rotation by an angle ψ about the x−axis, and T z , which represents a time dependent rotation at an angular velocity ω prc about the z−axis.Using these transformation matrices, we get the precessing velocity vector with and Here, the un-primed velocity components correspond to the non-precessing jet described in the previous section and t is time.Multiplying this unprimed velocity vector gives us a jet with a velocity vector having the same magnitude as the unprimed velocity vector, but now inclined by an angle ψ from the z−axis.Subsequently multiplying this resultant tilted velocity vector by T z , we rotate this velocity vector with an angular velocity ω prc in the x − y plane.Hence, we get a rotating jet which is precessing at a precession angle ψ = 10 o about the z−axis with a precession period P prc = 2π/ω prc = 5.We refer to this simulation setup with precessing injection as prc throughout the text.Note that the angular frequency of ω var corresponds to a sinusoidal variation of velocity in a particular direction i.e the jet axis.On the other hand ω prc refers to a precessional motion of the injected velocity vector.We present the input values of selected parameters which we have chosen for jet profiles as mentioned in Appendix A, along with other initial values for some jet parameters applied in the jet nozzle at t = 0, as well as the initial parameters for particle injection (discussed in the next section) in Table 2.Note that these values are common to all three injection nozzles as we prescribe the same initial boundary conditions for all the nozzles.

Particle Injection
So far we have discussed how we numerically govern the dynamics of the jet propagation.Here, we present how we implement particles into the jet, how we follow them in a Lagrangian approach, and how we accelerate them.
Overall, we apply the particle module for PLUTO invented by Vaidya et al. (2018).From the jet nozzle, along with the fluid material, we also inject Lagrangian macro-particles into our simulation domain, which, at each spatial point, always have the same velocity as the fluid.
These macro-particles should be understood as ensembles of micro-particles (electrons) that evolve and flow with the jet material and subsequently distribute over the jet volume.The energy distribution of these micro particles establishes the energy spectrum of the macroparticle.
In order to avoid ambiguity in notation, from now on in this paper we denote the macro-particle (the ensemble of micro-particles as particle and the constituent microparticles as electrons. For the particle injection we define a region of circular cross-section with radius r j and at a height of one pixel above the injection nozzle.This cross-section is divided into a cylindrical grid with 25 uniformly distributed radial divisions (along r) and 36 randomly distributed angular divisions (along φ) with 0 < φ < 2π, that are fixed in time.Here r and φ are usual cylindrical coordinates.Hence, we get a collection of 25 × 36 randomly distributed points from which we inject Lagrangian macro-particles.This injection is done after every two (numerical) time steps so that for a simulation run lasting 50 physical time steps, we get 10 6 such particles injected into the domain, ensuring a proper sampling of the jet material.
Each macro-particle represents an ensemble of nonthermal electrons with an energy distribution following a power-law with a chosen initial power-law index α = 6.We apply initially a power-law cut off at γ min = 10 2 and γ max = 10 8 .Obviously, when the jet evolves, particle acceleration and non-thermal cooling will affect the cutoff energies as well as the particle distribution evolves Note-Input parameters for all three injection nozzles (in code units).Note that for t = 0 these are identical, while for the time-varying nozzles the boundary conditions change over time.Shown in columns 2-8 are the input values of jet density ρj, ambient density ρa, gas pressure p, the magnetic field z-component Bzc, the magnetic field pitch angle δc, the angular velocity Ωc, and the jet Lorentz factor γc.The subscript c denotes the respective values at the z-axis, thus i.e. at r = 0. Also shown are (columns 9-14) the initial values of selected physical parameters (averaged over the nozzle), such as the proper sonic Mach number Ms, the proper Alfveénic Mach number Ma, the magnetization σ, the plasma β parameter, and the extractable power Px.Further, we show the total densities of kinetic energy E kin , the magnetic energy Emag, and the thermal energy E th .Columns 18-21 show the initial particle parameters such as the initial power-law index α, minimum and maximum energy cutoffs γmin and γmax, respectively, and the equipartition deviation factor concerning magnetic energy and particle energy.
away from the initial power law.With that, N (γ) measures the number of particles per unit volume with a Lorentz factor γ, while N 0 is defined by the number density of electrons N e as γmax γmin With a chosen normalization of m e c 2 = 1, the Lorentz factor of the electron γ then gives us the energy of the electron, where m e is the mass of the electron.Hence, in order to convert the particle energy into physical scales, the Lorentz factor of the electron γ must be multiplied with m e c 2 = 0.511 MeV.Finally, we quantify N e by using assumption of equipartition of energy density between magnetic field and radiating electron giving where U e is the energy density of electron and B eq represents the magnetic field corresponding to the equipartition.The magnetic field B dyn in our simulation, varies from B eq by a factor , which can be used to vary the ratio of the energy density between the magnetic field and the radiating electrons in our simulations.The value of N e can now be calculated by substituting Equation 15and Equation 16in Equation 17as Since the magnetic field in the jet injection nozzle varies following the profiles mentioned in Appendix A, we define B dyn as the average total field strength in the injection nozzle, giving B dyn = 6.27.With a choice of 2 = 8 × 10 −5 in our simulations, such that the particle energy is initially in sub-equipartition with the energy in the magnetic field, we get N e ∼ 0.29 and N 0 ∼ 1.46 × 10 10 .
Subsequently, these particles are advected along with the flow and fill the outflow volume.

Particle Acceleration & Evolution
After being injected at the jet base with an initial energy distribution, the particles evolve in the time and momentum space following a Lagrangian approach in which the particles follow the fluid velocity streamlines.As a result, the position of the particle r p is updated with time t according to the equation dr p /dt = v p = v f , where v p is the velocity of the particle and v f is the velocity of the fluid interpolated from the underlying Eulerian grid at the position of the particle.It must be noted that the particles here do not have any feedback on the fluid i.e. they do not change the fluid properties (e.g.density and velocity).
This procedure takes into account the acceleration of the particle through diffusive shock acceleration.Additionally, we take into account energy losses due to the adiabatic expansion of the jet, the synchrotron cooling of the particles by acceleration in the jet magnetic field, and IC scattering of CMB photons by the particles.
The DSA is often invoked to explain the presence of non-thermal particles in astrophysical systems.It results from repeated scattering of a particle off a coherent in-homogeneous magnetic field e.g. at a shock.To model DSA, we follow the novel approach developed in Vaidya et al. (2018).Here, the Lagrangian particles are flagged when entering a shocked region, defined as the area of cells with negative velocity divergence, together with a pressure gradient ∇p/p above some (normalized) threshold of 2. Following these conditions, the shocks in PLUTO are resolved with 3 grid cells.
Particles are followed carefully as they travel through the shocked region, and the pre-shock and post-shock states are defined as the states with the minimum and maximum value of pressure, respectively.Based on these states, we calculate the orientation of the magnetic field with respect to the shock normal and the strength of the shock, which is quantified by calculating the compression ratio in the shock rest frame for relativistic case (for a detailed discussion see Guess 1960;Lichnerowicz 1976).
Here, the subscripts 'u' and 'd' denote the upstream and downstream values, respectively.Using these parameters, the post-shock distribution of particles is updated as per the theory of DSA (see Vaidya et al. (2018) and references therein).
Radiative losses in our simulations occur due to synchrotron process, as well as up-scattering of the surrounding CMB radiation through the IC process.These processes are implemented in the particle module of PLUTO by solving the time-dependent particle transport equation (Webb 1989;Vaidya et al. 2018) where τ is the proper time, the first term inside the round bracket represents the losses from adiabatic expansion, γl represents radiative loss due to synchrotron and IC processes, and u f is the proper fluid velocity.The above equation is solved for each macro-particle advecting spatially with the fluid considering that the spectral distribution N (γ) of constituent micro-particles (i.e., electrons) is evolved at every step accounting for the above loss and shock acceleration processes based on local conditions interpolated from the fluid grid at particle's position.

JET DYNAMICS AND ENERGETICS
We first discuss general features of the jet evolution before we compare specifics of how the different jet nozzles affect the structure of the propagating jet.

General Jet Evolution
In Figure 1, we show the evolution of the dynamical structure of the jet for the different approaches of a time-independent injection (steady nozzle, std) in the top panel, a time-variable injection (var) in the middle panel, and a precessing nozzle (prc) in the bottom panel.We display 2D slices of the 3D density distribution superimposed by (projected) magnetic field lines in the x − z plane at y = 0 at dynamical times t = 20, 30, 40, and 50 (from left to right).
Following the injection from the nozzle, the steady jet (std) propagates at high speed in z−direction, filling the domain with low-density jet material before terminating at the jet head with a termination shock.Subsequently, the jet inflates a low-density cocoon, thereby forming shocks as a result of interaction with ambient medium and also due to interaction within the jet.The steady jet also encounters re-collimation shocks as it evolves.
Note, however, that the strength and even the existence of these re-collimation shocks depends critically on the numerical resolution applied (see Appendix B), which emphasizes the need of high-resolution simulations.The turbulent nature of the jet is evident from the small-scale fluctuations of magnetic field in the jet as compared to a smoother structure around the jet.Despite of its turbulent structure, the jet remains rather symmetric about the z-axis (but see our discussion in Section 3.2 concerning the lateral outflow structure).
We may distinguish four dynamically distinct regions in the simulation domain.First is the high velocity jet, producing strong shocks.Then there is the back-flow of shocked jet gas around the jet that also carries gas entrained from the ambient medium.Further out comes a cocoon structure with shocked ambient gas.All this is enclosed by the original ambient gas that is not yet affected by the jet flow and remains at rest.
An interesting feature, again a feature that is not seen in low-resolution runs, is the existence of a sharp discontinuity at z ∼ 15 in the jet.We call this feature the steady internal shock.This shock surface starts to develop around time t 30 and does not evolve much with time.Along the jet axis, this steady internal shock is located at z = 14.8 at t = 50.Essentially, it does not propagate much, as long as the termination shock at the jet head is within the computational domain.We understand this as resulting from the interaction between the jet beam flowing forward and the presence of an inner back-flow2 .To understand quantitatively the kinematics of the steady internal shock, we find that between time t = 40 and 50, the jet head advances from z = 26.5 to z = 32.5 i.e. a distance 6 (in code units).In the same time interval, the steady internal shock advances from z = 13.9 to z = 14.8, a distance of 0.9 (in code units).
Our explanation is supported by the observation that when running the simulation for longer times, such that the termination shock at the jet head leaves the domain, there is no material deflected from the jet head into the original jet beam.As a result, the steady internal shock  surface starts to evolve with time and fades away.An interesting manifestation of this strong, steady internal shock feature could be the formation of a knot in the emission map.Prime observational examples of such a steady knot feature is HST-1 in M87 (Biretta et al. 1999;Harris et al. 2006), and other knot-like features in AGNs discussed in Hervet et al. (2016); Lister et al. (2021).A detailed investigation of this feature will be presented in a subsequent study.
The jet with a variable velocity (var) forms a number of bow-shocks as it propagates, resulting in a skeletonlike appearance.These bow-shocks can be important for particle acceleration in the jet as we will investigate in further sections.This jet remains more collimated than the jet with steady injection, but shows a similar linear extent along the propagation direction.
The precessing jet (prc) is, obviously, highly asymmetric.Essentially, it shows a smaller linear extent in the propagation direction as compared to the steady jet at the same time, a natural consequence of the smaller momentum in this direction.We note, however, that the magnitude of the total velocity in both cases (precessing and steady jet) is same and hence, we inject same total kinetic energy density ρv 2 in the two cases.As a result of precessing velocity vector of the injected material, we see multiple finger-like extensions of the jet head.Note that this is a truly 3D feature, while the figures show slices of the 3D structure.This can also result in enhanced interaction with the ambient medium, which can result in more shocks and turbulence as we discuss in the next section.
We already mentioned the high degree of symmetry we observe in the jets injected without precession.This symmetry shows up despite the 3D treatment of the simulation and also despite the application of Cartesian coordinates, altogether demonstrating the quality of our numerical setup.Note however, that there are locations around the jet where the axial symmetry is partly broke (see discussion in next section).

Lateral Jet Structure
We further investigate the lateral structure of the jet by showing slices of the density distribution across the jet thus in the x−y plane along different distances along the jet and also for a chosen time t = 50 (see Fig. 2).
When we look at the full extent of the outflow structure at z = 15 (left panel), wee see a low-density (in blue) structure in the center -the high velocity jetsurrounded by a structure of intermediate density (yellow), that is moving with lower (and negative) speed, which we interpret as back-flow of shocked jet material.Both are surrounded by a high density and low velocity (in arbitrary direction) cocoon (in red) that contains shocked material from the ambient medium.Outside the cocoon we find the uniform ambient medium (in orange) -still in its initial state.
The cocoon shows a perfectly circular shape, implying that we have sufficient resolution to avoid instabilities generated as a result of injecting a cylindrical jet in a Cartesian grid.The other panels show a zoomed-in version of the density distribution at different values of z, all at the same time.These panels resolve the innermost jet structure close to the jet axis.
At z = 5, thus close to the injection nozzle, we see features resulting from active mixing of jet material with the back flow material.At the very center we see the original jet beam with density ρ = 1 (i.e.log ρ = 0).This is surrounded by a denser region formed as a result of shocks at the contact discontinuity at the jet beam boundary and the back-flowing jet material.
Note the elliptical region of small substructures that has a factor ten time higher density as compared to the jet beam .Interestingly the orientation of this ellipse is not aligned with the numerical grid.Also the number of these features is six, and not four as maybe expected from the Cartesian grid.We may understand this substructure as braids.Just outside these elliptically aligned braids we see two features ("red") that belong to a spiral structure.Indeed, further out the lower density material (greenish) seems also aligned in a spiral shape.Even further out, at r 3, again six features appear, in similar order as the inner braids.The physical reason of this ordered substructure is yet unknown to us.In principle, however, we may think about this as longitudinal re-collimation shock forming a braided structure.The exact origin of this ordered substructure may also be related on the profiles of variables injected from the nozzle, in particular from the profile of the jet rotation.However, we note that these braids are located outside of the jet stream (even at the location of the re-collimation shock), where there can hardly be shear between the jet rotation and the surrounding gas.Interestingly, while these braids break the toroidal symmetry of the outflow, their alignment in φ-direction shows still a degree of symmetry.Furthermore, further downstream, the jet flow recovers again a good degree of axisymmetry.
Further downstream, we show the zoomed-in version of the first plot of Figure 2 at z = 15.We see a denser region in the center, implying that we have captured the shocked area just above the sharp gradient in density and pressure at the steady internal shock.The substructure of two times six braids has disappeared, the jet forms a ring-like structure of rings of different density, with the central jet being densest.
Even further downstream from the nozzle, at z = 25, the jet approaches the jet head where we observe a high density due to the presence of the termination shock.

Further Jet Diagnostics
In order to see how the other dynamical properties of the jet vary over the domain, we show 2D slices of various jet variables at time t = 50, again concentrating on the steady jet nozzle, simulation std (see Fig 3).The pressure distribution shows, like the density distribution, a number of re-collimation shocks in the jet, as well as the strong steady internal shock feature at z ∼ 15, and the termination shock near the jet head.
These features are also visible in the plasma β = 2p/B 2 , the temperature T , and the proper sonic Mach number M s .We show the toroidal components of velocity and magnetic field, by plotting the y−component of v φ and B φ .Through this, we can visualize v φ and B φ coming out (in blue), and going in (red) the x − z plane.We note that by this we also see the axisymmetric nature of the fields demonstrated, strongly indicating that we have applied a proper resolution when injecting a cylindrical jet into a Cartesian domain.We also note the opposite signs of toroidal velocity and magnetic fields in the injection profiles (Equation A1).This reversal of sign is a typical consequence of jet launching and acceleration models (Blandford & Payne 1982;Ferreira & Pelletier 1995;Zanni et al. 2007).
The diagram for the vertical velocity v z shows an interesting feature of jet dynamics.As the jet develops, it expands and eventually undergoes a bifurcation at z ∼ 15.This results from the interaction with backflowing material that was reflected by the jet head.Subsequently, the propagating jet material flows around this inner back-flowing material.Note that the existence of such a feature is not straightforward, as the backflowing material can be expected to bypass the jet material around the jet.Here, the jet reflection is so strong, and of such a reflection angle that it comes back right towards the jet flow.We believe that the reason for this flow structure is the re-collimation shock right behind the jet head, leading to a reflection along the jet axis.
After all, this particular flow feature leads to the formation of the strong steady internal shock structure that we see in the distribution of various physical variables.The essential condition for this steady internal shock to form, we think is the existence of the strong shock at the jet head.In fact, if we let the jet head move out of our computational domain, such that there is no jet head to reflect the jet material (and hence, no back-flow), the strong steady internal shock fades away.
The strong head-on collision between the jet material and the inner back flow does not only provide a site for strong particle acceleration (see below), but also for neutrino production (see e.g.Britzen et al. 2019 for an observationally motivated scenario).From our modeling we indeed expect hadronic material to be entrained into the strong inner back-flow, as it carries shocked material.
Additionally, we show the distribution of a tracer for the origin of the material, where a value of unity tracks the material that originates from the jet nozzle, and a value of zero for material that does not belong to the jet, rather to the initial ambient medium.This enables us in principle to disentangle the jet material from the ambient gas.The intermediate values of the tracer, 0 < tr < 1, basically quantify the mixing state of two media.We see that as the jet propagates, it is expanding, also as a result of encountering the inner back-flowing material at z ∼ 15.Note that the back-flowing material (with negative velocity v z ) is indeed a mixture of jet material and ambient gas and and hence has a tracer value < 1.

Jet Nozzles
Before we compare the particle acceleration and nonthermal cooling processes in our simulations applying different jet nozzles, here we first investigate the differences in the dynamical features of jet propagation that are governed by the nozzle physics.Eventually, the jet dynamics will directly determine the acceleration and radiation properties of the jet.
To this end, we consider the contribution of various dynamical parameters only from the jet-cocoon structure and ignore the ambient medium.For this purpose, we define the jet-cocoon region as being composed by the grid cells where the speed v > 0.Then, the value of a certain parameter at a particular time is calculated by summing up the contribution from all the grid cells in the jet-cocoon region at that time.In Figure 4 we show for all the three nozzles, the time evolution of total energy density in different energy channels, viz the kinetic energy density E kin , the thermal energy density E th and the electromagnetic energy density E em , each defined as in equation 10.
We first see that the outflow dynamics produced by the std and prc setup nozzles show similar kinetic energy density, whereas the var setup nozzle evolves the outflow into a state of somewhat lower kinetic energy density.This is expected as the kinetic energy density we inject depends on the velocity (or Lorentz factor) injection of the jet, that is on average lower in the variable injection nozzle var as compared to the other two nozzles, which have the same velocity injection.Comparing the injection nozzles std and prc, both of which have the same kinetic energy density at initial times, we find that the precessing outflow has a higher kinetic energy at later times.Here, the jet is constantly changing its direction, and the interaction with the ambient medium is weaker.As a result, less of kinetic energy is converted to other energy channels for the case of the precessing outflow.
For the magnetic energy density we find that the steady jet approach leads to the highest level, followed by the variable and the precessing jet.This is a consequence of our injection profiles, as the magnetic field which is injected along the jet depends on the jet velocity (see Equation A1).Thus, the variable jet outflow is expected to be less magnetically energized than the steady jet.The precessing jet on the other hand, produces and contains jetlets 3 moving in different directions.These jetlets, on interacting, may lead to (numerical) reconnection of magnetic field.This is interesting from the point of energetics, as this will also lead to (numerical) heating of the jet.Due to the high resolution we apply, numerical diffusion and reconnection does play a minor role only.However, this energy will have in principle implications on the particle acceleration in such jets.In our approach of ideal MHD we cannot not model physical magnetic reconnection.
Another interesting point to note here is that although the steady and the variable jet are in equipartition with respect to the kinetic, thermal and electromagnetic energy densities, the precessing jet with similar injection as the steady jet is less magnetically dominant than the other jets we discuss.

Outflow Turbulence Level
In the approach of the present paper, acceleration of particles is achieved by shocks.Shocks are typically generated by a misaligned velocity field.Such velocity fields are a typical signature of turbulence.From Figure 1, we (literally) see that the level of turbulence looks different for the different jet nozzles we apply.We thus expect that the various nozzles we consider in this work drive turbulence in the domain differently based on their dynamics.We now want to quantify the level of turbulence generated by the various injection nozzles.
3 With jetlets we denote the jet portions that are injected from the nozzle at each time.These jetlets can also be understood as small volumes of material traveling together with the same bulk velocity.
In this context, we define the level of turbulence in a variable X at a certain grid cell 'i' located at (x, y, z) as where X i is the mean value of the variable X over all grid cells inside a cube of 20 3 grid cells, centered at (x, y, z).
Note that for calculating the local turbulence level δX T i , we only consider grid cells i where the tracer tr > 0.1.So, we do not take into account the grid cells that are in the ambient medium for calculating turbulence as well as for averaging.Also, we consider grid cells only above z = 2, thus above the injection nozzle only.As a measure of the total level of turbulence of the jet we then sum over all local measures of turbulence X T = Σ i δX T i (independently for all components x, y, z if X represents a vector field).
We can now compare the turbulence level in the velocity and the magnetic field for the different jet nozzles by calculating spatial averages of v T /v and B T /B, respectively.We find that both, the variable jet nozzle and the precessing jet, drive a higher level of turbulence in the velocity with an average v T /v 0.21, while for the steady jet we find average v T /v 0.15.Similarly, for the turbulence level in magnetic field we find B T /B 0.10, 0.17, 0.33 for the steady, variable, and precessing jet, respectively.
Another way of quantifying the level of turbulence is to compare the different energy channels viz the kinetic, thermal and electromagnetic energy densities E T kin , E T th , and E T em , respectively, as defined in equation 10.We show the exact values of these terms in Table 3.These can then be compared to the corresponding (integrated) energy terms (see Equation 10), which we also show in Table 3, providing the fraction of the total energy trans-formed into turbulent energy.After all, we expect that different nozzles govern a different fraction of the injected energy that is turned into turbulent energy.
Computing these fractions we find that 1.12% of the injected kinetic energy is turned into turbulent kinetic energy for the std nozzle setup.Similarly, the var and prc nozzles turn 1.81% and 1.12% of the injected kinetic energy into the turbulent kinetic energy, respectively.On the other hand, the std, var and prc nozzle convert 1.62%, 0.89% and 1.71% of their injected magnetic energy into the turbulent magnetic energy, and 0.23%, 0.22% and 1.26% of the injected thermal energy into the turbulent thermal energy, respectively.From this comparison, we conclude that the time-variation of the jet injection, in particular the precessing jet, leads to enhanced level of turbulence in the domain.
As a turbulent jet will lead to the formation of more number of shocks, turbulence has an implication on the acceleration of particles as well.Thus, from investigating the outflow turbulence levels we predict that the precessing jet will be most efficient in accelerating particles.
We show values of some characteristic dynamical output parameters at time t = 50 for the different injection setups, along with the parameters connected to particle acceleration (see next section) in Table 3.

PARTICLE ACCELERATION
In our simulations we inject Lagrangian macroparticles at the base of the jet which follow the flow of the fluid.

Propagation of Particles
In Fig. 5 we show the distribution of these particles in the jet for the simulation run std at time t = 50.All the particles injected are projected on the x − z plane, but in reality follow a 3D distribution.In addition, we color code these particles denoting the compression ratio η they have experienced (in the left panel) and the shock velocity in z−direction at that location (in the middle panel).
As can be seen, particles just below the strong steady internal shock at z = 15 have a high compression ratio η, suggesting the presence of very strong shocks at this position.The distribution of the vertical velocity v z for the particles follows exactly that for the jet material.Obviously this is expected for Lagrangian particles.Hence, we see a forward going jet with positive v z (in red) along with the back-flow where particles have a negative v z (in blue).
Figure 5 shows also the trajectories of few selected particles in the jet.We may follow their path from injection evolving over time.It can be seen that after injection at the jet base, the particles interact with the local fluid at their present location.Depending on the interaction, a particle may traverse over the whole jet length or may end up in the jet back flow or cocoon upon reflecting off a shock.
When a particle enters a shock, the strength of the shock is given by its compression ratio η as calculated by Equation 19.When a particle leaves the shock, this value of η is retained by the particle until it meets another shock with a different compression ratio.
Figure 6 shows the statistics of the compression ratio η at time t = 50 over the whole numerical domain for simulations applying different injection nozzles.Note that we detect compression ratios beyond the theoretical limit of 4 for RMHD shocks.We note here that this classical limit of 4 for MHD shocks is enhanced by the Lorentz factor in case of RMHD shocks and can have large values (Blandford & McKee 1976).However, the compression ratio η calculated in the shock rest frame, as is done in our study for a particular shock rest frame (Normal Incidence Frame), is limited by a values of 4 (Guess 1960;Appl & Camenzind 1988;Kirk & Duffy 1999;Summerlin & Baring 2012).
The number of shocks with compression ratio η > 4 in our simulations is low (<10%) and primarily arises due to numerical artifact when the macro-particle passes through a complex network of turbulent/interacting shocks in a short time and is unable to resolve multiple shocks.As the fraction of such particles is much less than those that sample the shocks accurately, such an update will have no effect on the quantitative comparison of compression ratio of shocked particles from different nozzles as described below.
We see that the variable jet nozzle var produces more strong shocks, compared to the precessing nozzle prc, which is similar to var concerning strong shocks.This can be explained as being the result of different dynamics of two injection nozzle setups.In the var nozzle setup, the slow moving jet material injected at an earlier time is repeatedly hit by the fast moving material injected at later time during various jet injection cycles.As a result, there are more number of such back-on collisions in the var setup as compared to the prc setup, where the injected jet continuously changes the direction of its motion.This results in more of the stronger shocks being formed in the var setup.On the other hand, the standard injection std simulation has shocks comparatively weaker than for the other two nozzles as a result of its time-independent injection.However, note that the number of shocks N shk a particle has encountered along its path is not represented by the compression ra- Note-Output parameters of simulation runs for the steady outflow setup std, the variable outflow setup var and the precessing outflow setup prc(in code units) at t = 50.Shown are the average values for the proper sonic Mach number Ms, the proper Alfvénic Mach number Ma, the magnetization σ, the plasma β, and the jet Lorentz factor γ.This is followed by the extractable power Px ≡ ExV /t, the outflow volume V (for cells where v > 0), the total energy density E V tot over the whole outflow volume, the integrated values of the kinetic energy density E kin , the thermal energy density E th , the electromagnetic energy density Eem, and the turbulent kinetic, thermal and electromagnetic energy densities E kin , E th , and Eem, respectively.The averaging and the integration is done for the cells where the tracer tr > 0.1.We also show for all the particles in the domain their total effective energy E ef p (i.e. total particle energy above γ = 10 3 ), their minimum Lorentz factor γmin, and their maximum Lorentz factor γmax. tio η, but by the number of times the η has changed its value.
In Figure 6 we show the histogram of the number of shocks the particles have encountered.We see that although the precessing nozzle prc leads to somewhat weaker shocks than the variable nozzle var, the overall number of shocks is substantially larger.In fact, the total number of shocked particles in the var setup is even lower than for the std setup.
Since both, the strength of the shock and also the number of shocked particles, are crucial in regard to the energetics of the particles, it is interesting to study which of these two is dominant.This can be investigated by comparing the particle spectra in different jets, as we will discuss in the next section.

Particle Energetics
The Lagrangian macro-particles (ensemble of electrons) which are injected at the jet base follow a powerlaw spectrum for the constituent electrons given by Equation 15.As the particles evolve with time (along their path), they may encounter shocks and become ac- celerated depending on the compression ratio they experience and the orientation of magnetic field with respect to the shock normal (see Vaidya et al. 2018).It should be noted that spectra of particles with compression ratio beyond a value of 4 (< 10% of all particles) is updated assuming an asymptotic limit of the power-law index α = 2.23 for ultra-relativistic shocks (see e.g.Kirk et al. 2000;Huang et al. 2023).
In addition to being accelerated while crossing shocks, the particles cool down due to various processes discussed in section 2.4.As a result, the overall spectrum of injected particles changes over time, from the initial steep power law to an energy spectrum indicating particle energy beyond and below the initial cut-off energies and with varying slope.
In Figure 7 we show the spectrum of one example particle (labeled P3 in the right panel of Fig. 5) for some simulation time steps.A good presentation of the energy spectra is to plot the particle energy density distribution N (γ)γ 2 as a function of particle energy γ.
When the particle is injected at time t 15, it has the initial spectral slope α = 6 and a maximum (thermal) Lorentz factor γ max = 10 8 .Note that the jet bulk motion (in z-direction) injected has a (maximum) Lorentz factor γ f of 10.
As the particle evolves, following the fluid, it cools down -with higher energy electrons of the macro particle cooling faster than the lower energy electrons.As a result, we see the high-energy tail of initial spectrum decaying and a subsequent decrease of γ max during t = 20.Before t 22, the particle encounters4 a shock leading to a change in the slope of the spectrum and an increase in γ max showing acceleration.The particle again enters a shock region at t 28 and 45, again modifying the spectral slope.However, as a consequence of rapid cooling of high-energy electrons, overall we see a decrease in the maximum Lorentz factor γ max .Finally, at t 50, The blue and green dashed lines show the slope of spectra with power-law index α = 2.6 and 3.1, respectively, for comparison with the particle spectra.The magenta dashed line shows the asymptotic limit of α = 2.23 for ultra-relativistic shocks.
the particle encounters a strong shock resulting in a flat spectrum and a sudden increase in γ max .Now we can compare the effect of the different injection nozzles on the particle population in the different jets they generate.A first point of view is to study the full spectrum of the jets from the different nozzles.We thus plot the combined spectrum of all the particles in the domain.This is done by adding the contributions of each of the particles to an overall spectrum.In Fig. 8 we compare the result for all three nozzles.
For comparison, we also show the initial spectrum (black dashed line) on the same plot.This is normalized by taking the initial spectrum of one particle and multiplying it with the total number of particles at t = 50.Hence, we have a one to one comparison between the initial and final spectrum for the different nozzles.
Essentially, we see that compared to the initially injected spectrum, the spectrum for all three nozzles has a significantly larger contribution from higher energies.Some particles are accelerated to such high energies going even above the initially defined γ max = 10 8 .The shape of the spectrum is represented by a broken powerlaw for the three nozzles, featuring a major peak at γ 100.This peak is a relic from the initial γ min = 100.
For all three nozzles, the spectrum falls for energies below the initial minimum cut-off energy γ min .We emphasize that we do not inject particles at these energies and, thus, the presence of particles at these low energies is a result of physical cooling that is implemented in our treatment.Interestingly, from the peak at γ = 100 till γ 10 3.5 , we see the spectrum declining following a power-law with a slope of 6.This part of the spectrum is indeed similar to the initial spectrum, which follows a power-law with an index of 6 as well.We note that particle are continuously being injected from the jet nozzle.So, the particles found following the initial spectrum at later times are newly injected particles, that are not yet accelerated or cooled down.Note that the cooling time for those low energy particles is long, and, seemingly, longer than their current age since injection.
In general, beyond γ 10 3.5 the spectrum flattens, following an ankle into another power-law.The deviation from the injection spectrum is resulting from a number of reasons.On the the pure increase in number of particles that are accelerated or have cooled.However, we see that this deviation is also different for the different injection nozzles, thus, strongly indicating that the overall particle acceleration indeed depends on the jet injection mechanism.
Beyond γ 10 3.5 , the spectrum for setup prc flattens the most, following a power-law with index α 2.3 till γ 10 4.5 where it steepens, forming a knee, and another power-law with α 3. The subsequent power law again steepens after passing a second knee at γ 10 8.5 , again followed by rapid steepening of the spectrum.
The spectrum for setup std, on the other hand, after flattening at γ 10 3.5 , follows a power-law with index α 2.6 till γ 10 7.5 , followed by a subsequent decline in the spectra.
Finally, the setup var follows a power-law with α 3 between γ 10 3.5−6 .Beyond this, the spectrum flattens till a γ 10 8 , and a subsequent rapid steepening.
Of the three injection nozzles we consider in this study, the precession setup prc is clearly the most efficient in accelerating particles, with γ going up to values even beyond ∼ 10 9 , corresponding to a physical energy reaching few hundreds of TeV.This is followed by he steady injection std, while the variable jet setup var seems the least efficient of the three (apart from a narrow range of energies around γ 10 8 , where it is more efficient than the std setup).
After all, and after having in particular compared the energy evolution in various channels for the dynamical evolution in the previous section, the interesting question remains about what the energy that is deposited in the particles.We thus now compare the total energy of all the particles in the domain for different injection nozzles setups for a certain point in time, t = 50.
For this purpose, we define the energy of a single (macro-) particle i as E p,i = γmax,i γmin,i N (γ i )γ i dγ i , where γ max,i and γ min,i denote the upper and lower energy cutoffs for that particle, respectively, and γ i is the Lorentz factor of the electrons of the particle i.Then, the total energy of all the particles is E p = Σ i E p,i .We find that for the precessing nozzle the particle energy is largest, E p 5.9 × 10 6 (in code units), followed by the variable nozzle with E p 5.5 × 10 6 , and the steady nozzle with E p 4.1 × 10 5 .
Note at this point that the total particle energy is not only the energy gained by the particles through shocks but also contains a portion of the injected particle energy which the particle has to start with.In addition, after the particle is injected, the high-energy end of its spectrum decays rapidly as the cooling time for higher energies is very small as compared to the dynamical timescales.Cooling, of course, also happens to the particles that become shock accelerated, thus a part of the gained energy is lost as well.
Overall, the injected energy of the particles, as well as the gained energy by shocks, can be quickly lost as a result of cooling, except at the low-energy regime that have a cooling time larger than the age of the jet (i.e.t = 50).This can be seen e.g. in the total particle energy spectra shown in Figure 8, which are dominated by the injected energy spectrum for energies lower than γ 10 3 .
Hence, to quantify the level of total energy gained by the particles, we define an effective total particle energy E ef p similar to the total particle energy E p , but taking into account only the energy contribution above a threshold of γ = 10 3 .Thus, we neglect the injected particle energy and only consider the contribution from the energy regime where the particles were re-energized by shocks.With that find that the prc nozzle simulation leads to the highest effective particle energy E ef p 8.93×10 3 , while the nozzles var and std lead to effective particle energies E ef p 0.95 × 10 3 and 1.44 × 10 3 , respectively.
Note that these values are only snapshots in time, and a thorough energetics of the particle energy must include considerations of cooling and re-acceleration along the particle path, involving some time integration of the relevant processes.We will come back to this issue in our follow-up paper.

Spectra of Jet Components
Astrophysical jets appear as structured beams of material, however, we don't know yet how these structures form.We don't know either whether the observed jet sub-structures indicate some radiation pattern or actual Here, the particles in the jet and the entrained material are defined by having a tracer value > 0.8 and < 0.8, respectively.
material substructure.Sub-structure observed in motion could thus be pattern motion or material motion.One way to disentangle the structure formation process is to study the emission processes of the material.Here, we compare the particle spectra of the different outflow components.
In order to investigate the role of the different jet components for the acceleration of particles, we have therefore studied the spectrum from these regions.
Along the length of the jet, we are interested in the spectrum near the jet base, the jet head, and at some intermediary region in the jet.Note that comparing spatially different areas along the jet also corresponds to investigate a time evolution of these features, as, due to the jet propagation, these areas have evolved for different periods of time after injection.Of particular interest is steady internal shock feature we have observed in simulation run std, where we can as well study its effect on the particle acceleration and spectrum.
In addition to substructures along the jet, we want to compare substructures across the jet.We thus distinguish between the particles in the jet and in the entrained material (a mix of the ambient and jet material), and investigate particle acceleration in these areas.In this regard, we define the jet and entrained material as the regions where the value of tracer is ≥ 0.8 or < 0.8, respectively.We show the spectral energy distribution for the particles in the jet and entrained material for all three nozzles in Figure 9.
For the steady injection nozzle, the jet is dominant in accelerating particles as compared to the entrained material at all particle energies.This is also true for the nozzle with variable injection.However, we find that the material entrained by the steady nozzle setup std is much more efficient than that in the var nozzle setup in order to accelerate particles to energies beyond γ 10 4 .
In contrast, the entrained material in the precessing jet is a much more efficient accelerator.Here the spectra from particles located in the entrained material have a contribution comparable to those located in the jet for the energy range γ = 10 3.5 to 10 6 .At the highest particle energies, γ ≥ 8.5, for the setup with the precessing nozzle the acceleration in the entrained material surpasses even that in the corresponding jet.This is ultimately a consequence of the dynamics of the precessing nozzle.The gas in the entrained material in this setup is being constantly hit and perturbed by the jet material, generating many shocks, and subsequent enhanced shock acceleration in the entrained material.
To study the effect of strong steady internal shock feature in the steady nozzle setup std, we analyze the spectra from the area blocks of region above, below and enclosing the steady internal shock.These blocks are essentially cylindrical slices with radius spanning the whole domain 5 and a height of ∆z = 1, centered at z = 2 (below the steady internal shock region, at the jet base), z = 14 (in the steady internal shock region), and z = 20 (above the steady internal shock region).Hence, the average spectrum from a block is the sum of spectra from each of the particles inside the block normalized to the number of particles in the block.Additionally, we also analyze the spectrum from the jet head which we define as the region where z > 32.
We expect a stronger acceleration for particles in the block located at z = 14 due to the presence of the strong shock around that location.We expect the same for the particles in the jet head, owing to the presence of a termination shock there.This is precisely what we see in Figure 10 where we show the spectra from these four blocks.We clearly see a larger contribution to the spectrum for energies beyond γ 10 4 from the blocks at z = 14 and z > 32 as compared to the two other blocks.
Interestingly, the steady internal shock is capable of accelerating particles to energies higher than γ 10 8 5 Although the size of the blocks span the whole domain, it should be noted that the particles are concentrated just in the jet and the cocoon and there are no particles in the ambient medium.(corresponding to electrons with energies of few tens of TeV), surpassing the efficiency of the termination shock as well.Since this shock at z = 14.8 is (i) steady and (ii) filled with high energetic particles, we expect to see a stationary knot-like feature in the emission map (which we will explore in a subsequent paper) at this location.
As discussed above, comparing spatially different areas along the jet also corresponds to study the timeevolution of the jet spectrum.This is also evident in Figure 10.At z = 2, we are in the close vicinity of the jet base from where the new particles are being continuously injected.As a result, the spectrum at that location closely resemble the initially injected spectrum with a peak at γ 100, a consequence of the initial minimum cutoff γ = 100.Beyond this, the spectrum follows a power-law with slope 6, same as the initially injected power-law index α = 6.As we move downstream along the jet, the particles enter the steady internal shock and are accelerated to high energies, as discussed.Further downstream, at z = 20, the particles have cooled down after leaving the strong shock to lower energies, decreasing both the upper and lower energy cutoffs.At last, particles meet the termination shock and are accelerated again to high energies.

CONCLUSIONS
We have applied the PLUTO code (Mignone et al. 2007) to perform 3D relativistic MHD simulations in combination with a hybrid module (Vaidya et al. 2018) that is capable of following Lagrangian macro-particles with the fluid.The injected Lagrangian macro-particles constitute ensembles of electrons that are close in space.Their initially defined energy distribution changes as they evolve in time, considering diffusive shock acceleration and cooling due to synchrotron radiation, inverse Compton processes and adiabatic expansion of the jet.
We investigate how physically different jet injection nozzles affect the particle spectrum that is generated as a consequence of the underlying jet dynamics.We have considered three different nozzle setups, a standard steady jet injection into an unstratified ambient gas, a time-dependent jet injection, and a jet nozzle undergoing precession.We adopt jet injection profiles applying an equilibrium solution from Bodo et al. (2019).This is essential in order to avoid artificial instabilities that mimic turbulence and, hence, lead to unphysical particle acceleration.We adopt a Taub-Matthews equation of state, as we consider a relativistic gas embedded in a non-relativistic ambient medium.
Considering our resolution study, we apply a resolution of 25 cells per jet radius, more than ever applied in 3D MHD jet simulations.This resolution is needed for both, resolving dynamically interesting features, as well as capturing the structure and strength of the shocks that are involved in acceleration.
While our simulations are run in code units and are applicable to any relativistic jet source, from the prescribed injection Lorentz factor of 10, the corresponding astrophysical scaling of our simulations would be that of a pc-scale AGN jet, with a size of our simulation domain of 10 × 10 × 20 pc.Naturally, micro-quasars with a similar initial speed can also be described, however, the physical magnetic field strength we have to assume for particle cooling, is chosen for AGN-like conditions.
Our results are as follows.
(i) For all nozzle geometries investigated the interaction of injected jet material with the ambient gas and with the previously injected gas leads to the formation of strong shocks.Essentially, the location and the properties of these shocks vary for the different nozzles.
(ii) The steadily injected jet forms a well known jetcocoon structure.In this simulation we find a special shock feature -a strong and stationary shock resulting from head on collision of jet and a back-flow along the jet axis inside the jet.The location and structure of this feature converges only for sufficiently high numerical resolution, i.e. 25 cells per jet radius.
(iii) The jet with variable injection evolves into a narrower cocoon, resulting in an overall structure more collimated than the steady jet.This jet forms multiple bow shocks as a result of various episodes of time-dependent activity, leading to a skeletal appearance.
(iv) The dynamics of the precessing jet is more turbulent as compared to the other nozzles we consider.As the jet axis change over time, multiple shocks are observed at the positions where the (multiple) jet heads interact with the ambient medium and with the jet.
(v) We quantify the level of turbulence mainly by the fluctuations of the energy density.We compare the energy budget involved in turbulent motions with the mean fluid motion and other energy channels.We find that about 1% of the bulk energies is turned into the turbulent energy channels.
(vi) On encountering a shock, the injected Lagrangian macro-particles change their spectral energy distribution depending on the shock compression ratio they experience.We find that the variable jet injection nozzle produces the strongest shocks (with more of the highest compression ratios), followed by the precessing jet and the steady jet.On the other hand, the total number of shocks is largest for the precessing jet, followed by the steady and the variable jet.The overall particle acceleration depends on the interplay between the available shock strength and the number of shocks they experience.
(vii) The macro-particles cool down due to synchrotron, inverse Compton and expansion losses.While particle acceleration happens instantaneously when a macro-particle crosses a shock, cooling depends on magnetic field strength and particle energy.For the chosen magnetic field strength, the cooling time for the high energy tail of the particle energy spectrum is well below the dynamical time scale of our jet simulations.
(viii) The total energy spectral distribution induced by the different nozzles shows a significant increase in the number of high-energy particles over time.The precessing jet leads to the formation of comparatively more shocks that are most efficient in accelerating particles, followed by the steady jet, while the shocks in the variable jet are the least efficient of the three.In the precessing jet electron are accelerated up to a few hundreds TeV with a power-law index of α 2.3 to 3.1.For the steady and the variable jets electron energies are up to a few tens of TeV, with spectra of similar but somewhat steeper slopes.From this we find that particle acceleration is governed more by the total number of shocks available than by the strength of the shocks.
(ix) We have compared various regions of the jet concerning their particle energy distribution.Compared to the jet-entrained ambient material the jet plays a dominant role for the particle acceleration for the steady nozzle and the variable injection nozzle.For the variable nozzle, the particle acceleration in the entrained structure is highly suppressed for higher electron ener-gies.For the precessing jet, entrained material and jet are comparable in particle acceleration efficiency, except at the very high-energy end, where the entrained material is even more dominant than the jet for acceleration.Here, particles in the entrained material are being reenergized through shocks continuously -a result of the highly turbulent jet dynamics.Hence, we expect precessing AGN jets to be the best sites for highly efficient particle acceleration.
(x) The total energy of accelerated particles in the precessing and variable jet is higher than the steady jet as a result of more shocks available due to the turbulent motion induced by the time-dependent injection.The energy the particles have gained till the end of the simulation can be up to 5% of the turbulent energy, i.e. below 1% of the injected jet energy.
(xi) We have studied the effect of the strong, stationary shock feature that we have discovered, on the particle spectrum.Like the jet termination shock, this strong stationary shock is a particular site where electrons are accelerated to very high energies of few tens of TeV.As these highly energized particles cool rapidly, we expect that this dynamical feature can be observed as a stationary knot of enhanced intensity in emission.The similarity to the famous HST-1 knot seems intriguing.
To summarize we have demonstrated and quantified how different jet nozzles convert the injected jet power to different levels of turbulence, which in turn leads to the generation of shocks that can accelerate electrons to energies up to tens and even hundreds of TeV.The resulting particle energies are about 10% of the turbulent energy and below 1% of the available jet energy.
In a subsequent paper, we will apply our results presenting synthetic emission maps of these jets.This will require a specific choice of the physical scaling for densities and pressure, and can be applied for different relativistic jet sources, such as AGN or micro quasars, as well as for different viewing angles.

ACKNOWLEDGMENTS
This project is financed through grant DFG-FOR5195 of the German Science Foundation DFG.B.V. acknowledges funding by the Max Planck Partner Group located at IIT Indore.R.D. acknowledges travel funds by the International Max Planck Research School for Astronomy & Cosmic Physics at the University of Heidelberg.We thank Andrea Mignone for providing and sustaining the PLUTO code.We acknowledge enlightening discussions with Karl Mannheim and John Kirk.All simulations were performed at the MPCDF computing center of the Max Planck Society in Garching (utilizing the compute clusters Vera, Raven & Cobra).We are grateful to an unknown referee for a detailed and helpful report that considerably improved the clarity of the paper.

A. JET INJECTION PROFILES
In our study, we adopt the equilibrium solution derived by Bodo et al. (2019), who solved for the radial component of the momentum equation for constant gas pressure p.
In cylindrical coordinates (r, φ, z) they find where One may define a pitch angle of the magnetic field, δ = rB z /B φ , the Lorentz factor in z−direction, γ z (r) = 1 + γ c − 1 cosh (r/r j ) 6 , the azimuthal velocity, and the function Here, r is the radial coordinate in the cylindrical coordinate system, γ c is the Lorentz factor at r = 0, r j = 1.0 is the radius of the injection nozzle, a = 0.6 is commonly denoted as the magnetization radius at which the magnetization of the jet nozzle has its maximum value, and Ω c is the angular velocity.Choosing a magnetization radius a less than the jet radius r j , we make sure that we do not have an infinite Lorentz force or a strong shear acting at the jet edges.Note that a subscript c denotes the value of a parameter at r = 0.This gives the following injection profiles within the cylindrical injection nozzle (see Bodo et al. 2019 for detailed derivation) where The parameter ζ governs the strength of rotation and E is the error function.Note that not all combinations of input parameters lead to a physical solution.In general, it must be required that B 2 z and B 2 φ are positive everywhere.
Overall, with these equations we can define a rotating, magnetized, relativistic jet of density ρ j and pressure p, launched from the injection nozzle, which is propagating along the z−direction and is governed by the choice of the parameters ρ j , γ c , Ω c , B zc , δ c and p.Our choice of these parameters for our simulations is summarized in Table 2.We show the profiles of the z− and φ−components of the velocity and magnetic field injected through the steady jet nozzle setup as a function of radius R, following the aforementioned equations in

B. RESOLUTION STUDY
A resolution study is essential for numerical simulation work.It is important for our study in particular as we are not only interested in the gross dynamical behavior of the jet flow, but also in its small-scale structure that is eventually responsible for particle acceleration.Thus, it is necessary to avoid (artificial) numerical effects in our study and achieve a sufficient resolution of turbulent structures.
We have performed a number of simulation runs with different numerical resolution applying the setup with the steady jet nozzle.To this end, defining the resolution as the number of grid cells inside the jet radius r j , we run simulations with a resolution of 10, 15, 20, 25 and 30 pixels per r j .In Fig B .1 we show 2D slices of the density distribution in the x − z plane (at y = 0) at time t = 50 for various resolutions for the simulation setup std, in order to compare the effect of resolution on the dynamics of the jet.
We see that with a low resolution of 10 grid cells per jet radius, we miss a lot of small-scale structure of the jet which could be essential for particle acceleration.Also, we do not see the formation of the strong steady internal shock structure, which is a very efficient site of particle acceleration as we discussed.The formation of the strong steady internal shock can be detected as we move to a higher resolution of 15 grid cells per jet radius.Still the position of the strong, steady internal  shock still depends on resolution, which converges only even higher resolution.As we go on from a resolution of 15 to a resolution of 20 cells per r j , we see additional shocks near the the jet head seen as an over-dense region of filamentary structure.In addition we see an upstream shifting of the strong, steady internal shock.
These differences to the lower resolution study are even more visible as we move to a resolution of 25 cells per r j .Here, we see the formation of two finger-like (crab-like?) structures near the base of the jet, which is not present at lower resolutions.We also see the location of the strong steady internal shock further upstream, although now the difference in position is smaller than between the lower resolutions.We thus think that we have rule out simulation runs applying a resolution of less than 25 cells per r j .On the other hand, runs with a resolution of 25 and 30 cells per r j look really similar, suggesting a convergence with regards to the resolution.
We also investigated convergence beyond the dynamical jet structure.In particular, we are interested in studying how the particle acceleration is affected, for example by effects that remain invisible by a comparison of the pure dynamical structure.In Fig. B.2 we show the histograms of the shock compression ratio η the particles have experienced until t = 50 for the simulation setup std for the different resolutions of 20, 25, and 30 grid cells per jet radius .We clearly see that the difference in the statistics between the runs with resolution 25 and 30 grid cells per jet radius is lower than that between the runs with resolution of 20 and 25 grid cells per jet radius.
We conclude that for both, the dynamical evolution and the particle evolution, the simulations with resolution 25 grid cells per jet radius are sufficiently converged compared to the resolution of 30 grid cells per jet radius, which justifies the use of a resolution of 25 grid cells per jet radius.

Figure 1 .
Figure 1.Distribution of density (in log scale) in the domain at a time (left to right) t = 20, 30, 40 and 50 (in code units) in the x − z plane for the simulation run std (upper panel), var (middle panel) and prc (lower panel).Magnetic field lines projected onto the plane are also shown (black lines).

Figure 2 .
Figure 2. Cross section of the jet density distribution (in log scale) at t = 50 in the x − y plane at z = 5, 15, 25 (see panel title) for the simulation run std.Note the different sizes and different color bars of the first and the other three panels.

Figure 3 .
Figure 3. Top: Distribution of (left to right) pressure (in log scale), plasma β (in log scale), temperature Θ (in log scale), and proper sonic Mach number Ms (in log scale) in the domain at a time t = 50 (in code units) in the x − z plane for the simulation run std.Bottom: Distribution of (left to right) vz, the y−component of v φ , the y−component of B φ , and the passive tracer tr in the domain at a time t = 50 (in code units) in the x − z plane for the simulation run std.

Figure 4 .
Figure 4. Time-evolution of global outflow energetics.Shown are integrated values of the kinetic energy density E kin , the thermal energy density E th , and the electromagnetic energy density Eem (from left to right).We consider computational values beyond t = 1 and integrated from above z = 2.

Figure 5 .
Figure 5. Distribution of particles with different compression ratio (left) and shock speed (vz−component) (middle) at a time t = 50.Right panel shows example trajectories of five representative particles from injection till t = 50, with the background showing the density distribution for better visualization.The colored dots on the trajectories represent the locations where the particle encounters shock.

Figure 6 .
Figure 6.Histograms of number of particles N with different compression ratio η (left) and total number of shocks N shk (right) for the simulation run std (green), var (blue) and prc (red) at a time t = 50 (in code units).

Figure 7 .
Figure 7. Energy spectrum of a particle (labeled P3 in the right panel of Fig. 5) at different times (as mentioned in code units) for simulation run std.

Figure 8 .
Figure 8. Energy spectrum of all particles in the domain for the simulation run std, var and prc at t = 50.The injected spectrum for each particle normalized to the total number of particles at t = 50 is shown by the black dashed line.The blue and green dashed lines show the slope of spectra with power-law index α = 2.6 and 3.1, respectively, for comparison with the particle spectra.The magenta dashed line shows the asymptotic limit of α = 2.23 for ultra-relativistic shocks.

Figure 9 .
Figure9.Energy spectra considering all the particles in the jet structure (solid line) and the surrounding entrained material (dashed line) for the simulation runs std (blue), var (orange) and prc (green) at a time t = 50 (in code units).Here, the particles in the jet and the entrained material are defined by having a tracer value > 0.8 and < 0.8, respectively.

Figure 10 .
Figure 10.Average energy spectra of all the particles in the cylindrical slices of height 1 (in code units) centered at z = 6 (below the steady internal shock region), z = 14 (in the steady internal shock region), and z = 20 (above the steady internal shock region).

Figure A. 1 .
Figure A.1.Variation of the z− and φ−components of the velocity and magnetic field in the steady injection nozzle as a function of radius R.

Figure B. 1 .
Figure B.1.Distribution of density (in log scale) in the domain at a time t = 50 (in code units) in the x − z plane (for y = 0) for the simulation run std at a resolution of (from left to right) 10, 15, 20, 25 and 30 cells per jet radius.

Figure B. 2 .
Figure B.2. Histogram of compression ratios η.Shown is the number of particles N in a specific range of η, normalized by the total number of particles in the domain Ntot (in log scale) for different compression ratio for the simulation run std at a resolution of 20 (green), 25 (blue) and 30 (red) at time t = 50.

Table 2 .
Input and Initial Parameters

Table 3 .
Output Characteristic Values