A Simulation Study of Ultra-relativistic Jets - II. Structures and Dynamics of FR-II Jets

The conservation equations for special RHDs can be written as

(e.g., Landau & Lifshitz, 1959). The conserved quantities, $D$, $M$, and $E$, are the mass, momentum, and total energy densities in the laboratory frame, respectively, which are related to the primitive variables, the rest-mass density, $\rho$, the fluid three-velocity, $v$, and the isotropic gas pressure, $p$, as

Here, $c$ is the speed of light, $\Gamma=1/\sqrt{1-(v/c)^{2}}$ with $v=|\mbox{\boldmath$v$}|$ is the Lorentz factor, $h=(e+p)/\rho$ is the specific enthalpy, and $e$ is the sum of the internal and rest-mass energy densities.

We adopt the EOS proposed by Ryu et al. (2006), which was named as the RC EOS in Paper I:

where $\Theta\equiv p/\rho c^{2}$ is a temperature-like variable. It closely approximates the EOS of the single-component perfect gas in relativistic regime, called the RP EOS in Paper I. Hence, our RC EOS correctly describes the fluids of both non-relativistic temperature ($\Theta\lesssim 1$) and relativistic temperature ($\Theta\gtrsim 1$).

Although the adiabatic index, $\gamma$, does not explicitly appear in the RHD equations, below we present $\gamma$ which is estimated with $\rho$ and $p$ from

It ranges between $\gamma=5/3$ for $\Theta\ll 1$ and $4/3$ for $\Theta\gg 1$.

The domain of our simulations is a rectangular box in the 3D Cartesian coordinate system, whose bottom surface defined by $z=0$ contains a circular jet nozzle with the radius of $r_{j}=1$ kpc, centered at $(x,y)=(0,0)$. At the onset of simulations, the box is filled with a static uniform background medium with the density, $\rho_{b}$, and the pressure, $p_{b}$. Without including dissipative processes such as radiative losses, the RHD equations in (1)-(3) are scalable for arbitrary length, time, and density, so are our simulations. However, we adopt the typical parameters for the hot ICM, that is, $n_{\rm H,ICM}\approx 10^{-3}$ cm^-3 for the hydrogen number density and $T_{\rm ICM}\approx 5\times 10^{7}$ K for the temperature, since we are interested in radio galaxies in galaxy clusters. Then, $\rho_{b}\approx 2.34\times 10^{-24}{\rm g}\cdot n_{\rm H,ICM}$ and $p_{b}\approx 2.3\cdot n_{\rm H,ICM}k_{\rm B}T_{\rm ICM}$, where $k_{\rm B}$ is the Boltzmann constant. We do not consider the stratification of ICM halos, because the jets expand out only up to $\sim 50-60$ kpc in our simulations. Aiming to explore the dependence of jet structures and flow dynamics on jet parameters with a simple background configuration, our study focuses on the early evolution of radio galaxies on several tens of kpc scales in the ICM.

In the next sections, the results of our simulations are presented in units of $r_{0}=r_{j}=1$ kpc, $v_{0}=c$, and $\rho_{0}=\rho_{b}=2.34\times 10^{-27}{\rm g~{}cm^{-3}}$ for the length, velocity, and density, respectively; the time and pressure are expressed in units of $t_{0}=r_{j}/c=3.26\times 10^{3}~{}{\rm years}$ and $p_{0}=\rho_{0}c^{2}=2.1\times 10^{-6}~{}{\rm erg~{}cm^{-3}}$, respectively. Then, the pressure of the background medium is given as $p_{b}/p_{0}=7.64\times 10^{-6}$ in the dimensionless unit, and its adiabatic index is $\gamma\approx 5/3$ with $p_{b}/\rho_{b}c^{2}\ll 1$.

The jet is injected through the nozzle at $z=0$ with $\rho_{j}=\eta\rho_{b}$, $p_{j}=\zeta p_{b}$, and $v_{j}$ or $\Gamma_{j}=1/\sqrt{1-(v_{j}/c)^{2}}$. Then, the jet power, $Q_{j}$, the amount of the kinetic plus internal energy (excluding the mass energy) injected into the background medium per unit time, is given as

This is the primary parameter that governs the morphological and dynamical evolution of jets through the interactions with the ICM, as mentioned in the introduction. The density and pressure ratios, $\eta$ and $\zeta$, are the secondary parameters, which may be combined into the momentum injection rate, or the jet thrust,

In this study, we specify the three parameters, $Q_{j}$, $\eta$, and $\zeta$; then, $v_{j}$ (and $\Gamma_{j}$) and $\dot{M}_{j}$ are determined.

Table 1 shows the parameters of the jet models considered here. The first column lists the model name, which consists of three elements, the exponents of $Q_{j}$, $\eta$, and $\zeta$. The three fiducial models in the first group have the same $\eta=10^{-5}$ and $\zeta=1$, but different $Q_{j}$’s, Q45, Q46, and Q47. They are intended to represent the low-power, intermediate-power, and high-power jets of FR-II radio galaxies, respectively, and the Lorentz factor of jet flows ranges $\Gamma_{j}\approx 7.3-71$. Those attached with H denote the high-resolution models. High-resolution simulations have been run only for the Q46 and Q47 models, owing to the much longer computational time required for the Q45 model. The models in the second group include the jets of higher densities with $\eta=10^{-4}-10^{-3}$, while those in the third group include over-pressured jets with $\zeta=10$. In our jet models, for a fixed value of $Q_{j}$, the higher density (larger $\eta$) or higher pressure (larger $\zeta$) of the jet leads to smaller $v_{j}$ or $\Gamma_{j}$, whereas larger $\eta$ or smaller $\zeta$ results in larger $\dot{M}_{j}$. Note that for very high power jets with large $\Gamma_{j}$, $\dot{M}_{j}\approx Q_{j}/c$. The adiabatic index of injected jet material is fixed by the ratio of $\eta/\zeta$; for $\eta/\zeta=10^{-5}$, the temperature is relativistic with $p_{j}/\rho_{j}c^{2}=0.764$ and $\gamma=1.43$, while for $\eta/\zeta=10^{-3}$, $\gamma\approx 5/3$. We point that the temperature of shocked jet material in the hot spot is higher than that of injected material, and the adiabatic index there approaches $\gamma=4/3$ in some of our models.

For stable FR-II type jets, assuming a balance between the jet ram pressure and the background pressure, the advance speed of the jet head was estimated approximately as

(Martí et al., 1997). Here, $\eta_{r}=(\rho_{j}h_{j}\Gamma^{2}_{j})/(\rho_{b}h_{b})$ is the relativistic density contrast. In general, $\eta_{r}>\eta$, but $\eta_{r}$ approaches $\eta$ for non-relativistic jet speeds and internal energies. We define the jet crossing time as

which can be used as a characteristic timescale for the jet evolution. The values of $v_{\rm head}^{*}$ and $t_{\rm cross}$ are given in the 8th and 9th columns of Table 1. Both $v_{\rm head}^{*}$ and $t_{\rm cross}$ are fixed mainly by $Q_{j}$; for instance, $t_{\rm cross}$ for Q45 is $6-7$ times longer than $t_{\rm cross}$ for Q47. Although $v_{\rm head}^{*}$ and $t_{\rm cross}$ depend also on $\eta$ and $\zeta$ (and $\dot{M}_{j}$), the dependence is weak. Below the results of jet simulations are described in terms of $t_{\rm cross}$.

The jet flow is directed upward mostly along the $z$-axis with $v_{z}\approx v_{j}$. However, to break the rotational symmetry, a slow, small-angle precession with period $\tau_{\rm prec}=10~{}t_{\rm{cross}}$ and angle $\theta_{\rm prec}=0.5^{\circ}$ is applied to the jet velocity. Also to ensure a smooth launching of jets, preventing possible developments of unphysical structures in the start-up, the jet velocity is modified with a window function as

in the early stage of simulations. Here, $r=\sqrt{x^{2}+y^{2}}$ is the radial distance from the jet axis. The windowing is applied for one jet crossing time, $t=t_{\rm cross}$, and turned off afterward.

For the models listed in Table 1, simulations have been carried out using the newly developed 3D RHD code presented in Paper I. The version used for this study includes (1) a 5th-order accurate, finite-difference WENO scheme, WENO-Z (Borges et al., 2008), (2) a 4th-order accurate time-integration method, the strong stability preserving Runge–Kutta (SSPRK) (Spiteri & Ruuth, 2003), and (3) a 4th-order accurate averaging of fluxes along the transverse directions in smooth flow regions, which improves the performance in multi-dimensional problems involving complex flows (Buchmüller & Helzel, 2014). In addition, to suppress the carbuncle instability, which often appears at the bow shock surface, the code incorporates a modification of eigenvalues for the acoustic modes; the local sound speed is limited to ${c}^{\prime}_{s}=\min(\phi|{v}_{x}|,{c}_{s})$ in the calculation of eigenvalues (Fleischmann et al., 2020). For the Q46 and Q47 models, $\phi=10$ is used, while $\phi=5$ is used for the Q45 models. For the CFL condition, ${\rm CFL}=0.8$ is used. The code implements Message Passing Interface (MPI) for parallel computing. The details of numerical schemes and performance tests can be found in Paper I.

Simulations have been performed in the boxes elongated along the $z$-direction, which consist of either $(400)^{2}\times 600$ (for default models) or $(600)^{2}\times 1200$ (for high-resolution models) uniform grid zones. The jet radius $r_{j}=1$ kpc occupies either $N_{j}=12$ grid zones (for default models) or 18 grid zones (for high-resolution models). Hence, the size of grid zones is $\Delta x=1/12$ kpc or $\Delta x=1/18$ kpc, and the volume of the simulation domain is $(33.3)^{2}\times 50$ kpc³ or $(33.3)^{2}\times 66.7$ kpc³. The number of grid zones and the resolution in terms of the jet radius are listed in the 10th and 11th columns of Table 1.

It is known that the imposed condition at the bottom boundary of $z=0$ affects the properties of simulated jets (e.g., Donohoe & Smith, 2016). The commonly used condition is either “outflow” or “reflection”. With the outflow condition, some of the material along with the energy and momentum escapes from the simulation domain (see Figure 4 below), while it is conserved with the reflection condition. As a consequence, for instance, the jet morphology turns out to be more elongated with the outflow condition than with the reflection condition (e.g., Koessl & Mueller, 1988). We apply the continuous outflow condition to all the six faces of the simulation domain, including the $z=0$ plane except at the jet nozzle. Simulations stop when the bow shock reaches either the top $z$-boundary or the side $x$ and $y$-boundaries, and the end time of simulations, $t_{\rm end}$, is given in the last column of Table 1.

Previous numerical studies have shown that the characteristic morphology of light, relativistic, FR-II-type jets may include the following features: (1) a terminal shock (or “working surface”) at the head of the jet where the flow is halted and reversed, (2) a bow shock plowing through the denser background medium and representing the outer surface of the entire jet structures, (3) a board cocoon of the shocked jet material flowing backward, (4) the shock-heated background medium encompassing the cocoon, and (5) recollimation shocks formed in the jet spine flow (e.g., English et al., 2016; Li et al., 2018; Perucho et al., 2019). In reality, a stable distinct terminal shock may not appear, because the flows in the head become turbulent (e.g., Hardcastle & Krause, 2013; Li et al., 2018). Note that the cocoon filled with relativistic plasma is expected to be observed as radio lobe, so the two terms, cocoon and lobe, are often used interchangeably in describing jet structures.

In our simulated jets, we classify the structures bounded by the bow shock into three regions: (1) the highly under-dense, jet spine flow with $\rho<0.1~{}\rho_{b}$ and $v_{z}>0.1~{}c$, which is injected from the nozzle and keeps focused by recollimation shocks, (2) the low density, backflowing jet material with $\rho<0.1~{}\rho_{b}$ and $v_{z}<0.1~{}c$, which is halted and reversed at the jet head, and (3) the higher density, shocked ICM with $\rho>0.1\rho_{b}$ behind the bow shock. The left panel of Figure 1 illustrates those regions in one of our jet models, Q46-$\eta 5$-$\zeta 0$-H, moving from the inside to the outside, the jet spine flow, the backflow, the shocked ICM, and the bow shock, which is the outermost surface (see also the top panels of Figure 2).

Figure 2 shows the two-dimensional (2D) slice images of $\log\rho$ and the adiabatic index, $\gamma$, for the three fiducial models, Q45-$\eta 5$-$\zeta 0$, Q46-$\eta 5$-$\zeta 0$-H, and Q47-$\eta 5$-$\zeta 0$-H. The density distributions in the top panels demonstrate how the morphology of jets depends on the jet power $Q_{j}$. Comparing the three fiducial models, one can see that in the models with higher $Q_{j}$ (and larger $\Gamma_{j}$), the jet advances faster in terms of $t_{\rm cross}$, resulting in more elongated structures. Note that the jet crossing time, $t_{\rm cross}$, is even shorter in the models with higher $Q_{j}$ (see Table 1). On the other hand, more extended cocoons of backflow, filled with shocks and turbulence, develop in the models with lower $Q_{j}$. The backflow mixes with the shocked ICM, as also can be seen in the distributions of $\gamma$ in the bottom panels. At the jet head, the adiabatic index is close to $\gamma=4/3$ (blue); in the backflow, it increases due to the turbulent mixing as well as the adiabatic expansion, and smoothly merges to the value of the shocked ICM, $\gamma=5/3$ (red).

In the right half ($x>0$) of each panel in Figure 3, $\log\rho$ is plotted for jets with different $\eta$ and $\zeta$. Overall, while the influence of $\eta$ and $\zeta$ on the jet morphology is less significant than that of $Q_{j}$, as already known from previous studies (see the introduction), the figure shows that those secondary parameters affect some properties of the jets, such as the advance speed, in our simulations. In the high-power Q47 models, the jet advance speed is not very sensitive to $\eta$ and $\zeta$, while the jet of Q47-$\eta 3$-$\zeta 0$ propagates a bit faster where the momentum injection rate, $\dot{M}_{j}$, is slightly larger (see Table 1). In the lower power models (Q46 and Q45), the jet expansion rate differs somewhat for different $\eta$ and $\zeta$. For instance, among the four Q46 models, the jet of Q46-$\eta 3$-$\zeta 0$ with the largest $\dot{M}_{j}$ advances the fastest. Although the other three Q46 jets have similar $\dot{M}_{j}$, they propagate with somewhat different speeds. We find that the amount of the $z$-momentum contained in the jet spine flow, $M_{z,\rm{jet}}$, which pushes the jet head forward, turns out to be different. Figure 4 shows the accumulated $M_{z,\rm{jet}}$ as a function of $t/t_{\rm cross}$ for the four Q46 models. Even though $\dot{M}_{j}$ is similar among the three models, $M_{z,\rm{jet}}$ may develop differently owing to different dynamical evolutions. For instance, in the Q46-$\eta 4$-$\zeta 1$ model where the jet pressure is higher, the amount of the material that escapes from the simulation domain through the $z=0$ boundary is less than in the Q46-$\eta 5$-$\zeta 0$ and Q46-$\eta 5$-$\zeta 0$ models, because the jet expands more laterally (see the middle panels of Figure 3). The $z$-momentum in the whole jet structure as well as $M_{z,\rm{jet}}$ in the jet spine are smaller in the Q46-$\eta 4$-$\zeta 1$ model, since less material with negative $z$-momentum escapes from the system. Hence, the difference in the advance speed should be attributed at least partly to the outflow condition we impose at the $z=0$ boundary. A similar trend is found for the Q45 models.

In the left half ($x<0$) of each panel in Figure 3, the vertical velocity, $v_{z}$, is plotted. The upward-moving jet spine flow is shown in red, while the downward-moving backflow is shown in blue. The shocked ICM surrounding the cocoon has relatively small vertical velocities ($v_{z}\sim 0$, white), while it expands mostly in the lateral direction behind the bow shock. The interface between the jet spine flow and the backflow and also the interface between the backflow and the shocked ICM become turbulent via the Kelvin-Helmholtz (KH) instability due to strong velocity shear. Moreover, the working surface is not very apparent owing to turbulent flows and the small precession in the injected flow at the jet nozzle.

A notable point in the $\log\rho$ images of Figures 2 and 3 is that the bow shock surface includes kink-like structures in the Q46 and Q47 jets, while it is relatively smooth in the low-power Q45 jets. In addition, structures resembling herringbone patterns appear in the density of the shocked ICM behind the bow shock for the Q46 and Q47 jets. Such structures were observed in previous simulations of high-power jets using high-accurate codes (e.g., see Perucho et al., 2019). Below, we will see similar patterns also emerge in the vorticity distribution. We expect that the development of these structures would be the consequence of the interaction between the bow shock and the turbulent flows in the cocoon.

As the jet propagates into the background medium, while the whole jet-induced structure expands, the head where the jet spine flow stops advances. After the initial adjustment, the advance speed of the jet head is expected to approach the analytic estimate in Equation (11), roughly $v_{\rm head}^{*}\approx 0.05c$, $0.1c$, and $0.3c$ for the Q45, Q46, and Q47 models, respectively (see Table 1). The left panel of Figure 5 shows the time evolution of the advance speed, $v_{\rm head}$, which is determined with the actual position of the jet head in simulations, for all the models considered here. The values of $v_{\rm head}$ fluctuate around $v_{\rm head}^{*}$, but after $t/t_{\rm cross}\sim 10-20$, they tends to approach asymptotic values. We find that the asymptotic values are roughly $v_{\rm head}\approx 0.025c$, $0.1c$, and $0.45c$ for the Q45, Q46, and Q47 models, respectively. That is, the numerically estimated asymptotic values are somewhat larger than $v_{\rm head}^{*}$ in the high-power Q47 models, while they are smaller than $v_{\rm head}^{*}$ in the low-power Q45 models. In the Q46 models, the asymptotic values are close to $v_{\rm head}^{*}$.

The shape of the lobe (cocoon) may be quantified with the axial ratio, $\mathcal{A}\equiv\mathfrak{L}/\mathcal{W}$, where $\mathfrak{L}$ is the vertical length of the cocoon and $\mathcal{W}$ is the lateral width at its midpoint (Hardcastle & Krause, 2013). The middle panel of Figure 5 shows $\mathcal{A}$, which is measured from the simulation results. The axial ratio undergoes variations due to the competition between the longitudinal advancement and the lateral expansion. Overall, $\mathcal{A}$ is larger if the jet is more powerful. In the high-power Q47 models, the shape of the cocoon is highly elongated with $\mathcal{A}\sim 6$ or even larger, which still increases at $t_{\rm end}$ in our simulations.

In the Q45 and Q46 models, on the other hand, $\mathcal{A}$ on average increases up to $t/t_{\rm cross}\sim 20$ and then approaches asymptotic values. In these models, while the jet advances slowly, the over-pressured cocoon expands laterally, resulting in smaller asymptotic axial ratios; roughly $\mathcal{A}\sim 2$ and $\sim 4$ for the Q45 and Q46 models, respectively. As noted above, $v_{\rm head}$, the increment speed of $\mathfrak{L}$, is about four times larger in the Q46 models than in the Q45 models. The lateral expansion speed, $v_{\rm lateral}$, could be estimated from $\mathcal{W}$ shown in the right panel of Figure 5, where the slope gives $v_{\rm lateral}\times t_{\rm cross}$. We find that $v_{\rm lateral}$ is about twice larger in the Q46 models. Since $v_{\rm head}/v_{\rm lateral}$ is about twice larger, we get $\mathcal{A}$ that is about twice larger in the Q46 models than in the Q45 models. In general, both $v_{\rm head}$ and $v_{\rm lateral}$ tend to increase with the jet power.

Comparing the thin and thick green lines for Q46-$\eta 5$-$\zeta 0$ and the thin and thick blue lines for Q47-$\eta 5$-$\zeta 0$ in Figure 5, we see that the overall morphology of the jets is well converged in the default and high-resolution simulations (see also Figures 2 and 3). Yet, we expect that finer structures and more flow motions develop in smaller scales in higher resolutions (see the next section).

As noted above, the jet is injected to the uniform background in our simulations, assuming that the jet propagation is confined within the cluster core region. In the case of giant radio galaxies that expand out to several 100 kpc into stratified halos, the jet head is expected to be decelerated due to the lateral expansion; then, the evolutionary trend of $\mathcal{A}$ may differ somewhat from that shown in Figure 5.

In short, the morphology of the jet including the shape of the cocoon is primarily governed by the jet power, $Q_{j}$, and less dependent on $\eta$ and $\zeta$. Hence, we present mainly the fiducial models with different $Q_{j}$ in describing jet flow dynamics in the next section. The other models will be used to examine the parameter dependence of the problem.

Before we describe jet flow dynamics, we examine the kinetic and internal energies contained in the different regions of the jet. We point that while the RHDs is formulated based on the energy-momentum conservation, the strict decomposition of the energy into the kinetic and internal energies is not plausible (e.g., Landau & Lifshitz, 1959). Yet, following previous works (e.g., English et al., 2016), we divide the energy in Equation (6) as follows:

where the three terms in the right hand side could be regarded as the kinetic, internal, and mass energies, respectively. In terms of the adiabatic index in Equation (8), the internal energy can be written as $\gamma\Gamma^{2}p/(\gamma-1)-p$. In the non-relativistic limit, the kinetic and internal energies reduce to the usual forms, $(1/2)\rho v^{2}$ and $p/(\gamma-1)$, respectively.

Panel (a) of Figure 6 shows the kinetic plus internal energy, $E_{\rm KE+IE}$ (excluding the mass energy), contained in the region encompassed by the bow shock surface, relative to the accumulated energy injected through the jet nozzle (labeled as theory), for three fiducial models. Note that the initial internal energy of the background material should be counted for the exact estimation of the expected energy, but it is small in our jet models, $\sim 5~{}\%$, $\sim 1~{}\%$, and $\sim 0.2~{}\%$ of $E_{\rm KE+IE}$ around $t_{\rm end}$ for the Q45, Q46, and Q47 models, respectively. We find that the energy inside the bow shock surface is smaller than the injected energy, because a part of the energy escapes through the outflow boundary at $z=0$ (see, also English et al., 2016). The fraction of the escaped energy is greater for lower $Q_{j}$, since broader, more turbulent cocoons develop.

Panels (b)-(d) of Figure 6 show the kinetic and internal energies in the different regions of the jet for the three fiducial models. The energy partitioning differs for different models. Roughly, inside the bow shock surface, the kinetic energy is several times smaller than the internal energy. The kinetic plus internal energy in the cocoon (including both the jet spine flow and backflow) is somewhat smaller than that in the shocked ICM. The kinetic energy in the cocoon, which is the manifested quantity in the flow dynamics described below, is always a small fraction of the jet energy; it is estimated to be $\sim 5~{}\%$, $\sim 3~{}\%$, and $\sim 2~{}\%$ of $E_{\rm KE+IE}$ inside the bow shock surface around $t_{\rm end}$ for the Q45, Q46, and Q47 models, respectively.

The jet-induced structure contains two types of shock surfaces in our simulations: (1) distinct, connected surfaces such as the bow shock and recollimation shocks, and (2) less prominent, somewhat disordered shock surfaces in turbulent flows such as the jet spine flow, backflow, and the shocked ICM (see Figure 8).

Each shock surface is composed of many “shock zones” (numerical grid elements), which are identified in a post-processing step by applying a widely used algorithm (e.g., Ryu et al., 2003). Along each coordinate axis, grid zones are tagged as “shocked”, if they satisfy the following conditions: (1) $\mbox{\boldmath$\nabla$}\cdot\mbox{\boldmath$v$}<0$, i.e., the locally converging flow, and (2) $\Delta p/p>\epsilon_{p}$, i.e., the pressure jump in the adjacent zones larger than $\epsilon_{p}p$. Here, $\epsilon_{p}$ is a free parameter that should depend on the minimum Mach of the shocks to be identified. In our simulations, each shock transition spreads typically over two to three “shocked” zones, so the zone with the minimum value of $\mbox{\boldmath$\nabla$}\cdot\mbox{\boldmath$v$}$ is identified as a “shock zone”. The preshock and postshock states are then estimated across the shock transition. With the density and pressure in the preshock and postshock regions, the Mach number along each coordinate axis is calculated (see below); the Mach number of the shock zone is obtained as $M_{s}=\max(M_{s,x},M_{s,y},M_{s,z})$.

We here find shocks with $M_{s}\geq 1.01$, expecting that shocks with $M_{s}<1.01$ would be dynamically unimportant (see the next subsection). For non-relativistic shocks of $M_{s}=1.01$, $\Delta p/p\approx 0.025$ across the shock transition. So we set $\epsilon_{p}=0.005$, since the shock transition spreads over grid zones.

In Newtonian HDs, $M_{s}$ can be estimated with the pressure ratio across the shock, $p_{2}/p_{1}$, using $p_{2}/p_{1}=(2\gamma M_{s}^{2}-\gamma+1)/(\gamma+1)$, where the adiabatic index $\gamma=5/3$ for thermally non-relativistic, monatomic gas. Hereafter, the subscripts 1 and 2 denote the preshock and postshock states, respectively. In RHDs, however, $M_{s}$ cannot be determined by $p_{2}/p_{1}$ alone. In addition, $M_{s}$ is a quantity that depends on the frame where it is obtained.

From the density, momentum, and energy conservations across the shock, the shock speed can be expressed as

in the shock-rest frame. On the other hand, for instance, if the shock moves with $v_{t}$ along the direction transverse to the shock normal, the shock speed is modified to $v_{s}/{\Gamma_{t}}$, where $\Gamma_{t}=1/\sqrt{1-(v_{t}/c)^{2}}$. It is technically non-trivial to find the shock-rest frames for all the shock zones in our simulations, and hence it is practically difficult to estimate the velocity with which each shock zone moves in the computational frame. Thus, we calculate the Mach numbers of shock zones in the shock-rest frame, rather than in the computational frame. This would introduce uncertainties in the characteristic properties of identified shocks especially in the jet spine flow, as some of those have $v_{s}$ close to $c$. Such effects should not be substantial for shocks in other parts of the jet (see Figure 9).

For identified shock zones, the estimation of $v_{s}$ with Equation (15) using numerical values of $h_{1}/h_{2}$ and $\rho_{1}/\rho_{2}$ would not be robust, since those two ratios are close to unity and not very sensitive to $M_{s}$, especially at weak shocks. Instead, we find that it is more reliable to use the ratio $p_{2}/p_{1}$ in estimating $M_{s}$. Again, from the conservations across the shock, we can get

Then, for given values of $p_{2}/p_{1}$ and $\Theta_{1}\equiv p_{1}/\rho_{1}c^{2}$, the ratios, $h_{1}/h_{2}$ and $\rho_{1}/\rho_{2}$, can be obtained using Equations (7) and (16), and the shock speed can be calculated from Equation (15). With the RC EOS, the sound speed of the preshock gas is given as

Hence, $M_{s}\equiv v_{s}/c_{s,1}$ can be calculated for given values of $p_{2}/p_{1}$ and $\Theta_{1}$. In practice, we have built a numerical table for $M_{s}$ as a function of $p_{2}/p_{1}$ and $\Theta_{1}$. The left panel of Figure 7 shows $M_{s}$ in the 2D parameter space of $(p_{2}/p_{1},\Theta_{1})$. We use this table for the estimation of $M_{s}$ with the numerical values of $p_{2}/p_{1}$ and $\Theta_{1}$ obtained for shock zones in simulated jets.

We then estimate the kinetic energy dissipation rate at the shock, $\Delta\Phi_{\rm shock}$, from the difference between the entering and leaving kinetic energy fluxes across the shock. It can be written in the shock-rest frame as

where $v_{1}=v_{s}$, and $\Gamma_{1}$ and $\Gamma_{2}$ are the Lorentz factors of the preshock and postshock flow speeds, $v_{1}$ and $v_{2}$, respectively. It reduces to $(1/2)\rho_{1}v_{s}(v_{s}^{2}-v_{2}^{2})$ in the non-relativistic limit. Below, we use $\Delta\Phi_{\rm shock}$ to estimate the total amount of the jet energy dissipated at shocks in the jet-induced flows. This quantity is also frame-dependent. For instance, if the shock moves with a transverse velocity of $v_{t}$, it is given as $\Delta\Phi_{\rm shock}\Gamma_{t}$. We here employ $\Delta\Phi_{\rm shock}$ in the shock-rest frame, rather than the one in the computational frame; hence our estimation of the shock dissipation should be considered only approximate.

As in the case of $M_{s}$, we have built a numerical table for $\Delta\Phi_{\rm shock}/(\rho_{1}c_{s,1}^{3})$ as a function of $p_{2}/p_{1}$ and $\Theta_{1}$, which is plotted in the right panel of Figure 7. We use this table for the estimation of $\Delta\Phi_{\rm shock}$ with the numerical values of $p_{2}/p_{1}$,$\Theta_{1}$, and $\rho_{1}$ obtained for shock zones in simulated jets.

Panel (b) of Figure 1 illustrates the 3D spatial distribution of shock zones formed in the different regions of the Q46-$\eta 5$-$\zeta 0$-H jet. Figure 8 shows the 2D distributions of shock zones for three fiducial models. The bow shock surface surrounds the entire jet-induced structure and separates the shocked ICM from the background medium. While the so-called recollimation shocks appear along the jet spine, the first one is the most distinctive and its surface has a conical shape that stretches further upward from the jet nozzle in the models with higher $Q_{j}$. Although the surfaces of the bow shock and the first recollimation shock are composed of many shock zones with different $M_{s}$, we find that those shocks can be characterized with typical values of $M_{s}$ (see below). On the other hand, the surfaces of shocks formed in turbulent flows such as the backflow and shocked ICM are more chaotic and less distinct. While there are many of them, the sizes of connected shock surfaces are much smaller than that of the bow shock surface. We note that the termination shock does not clearly appear in our model jets due to the turbulence in the head and backflow. Below, we quantify the properties of shock zones to examine the characteristics of the shocks formed in the jet-induced flows.

Panels (a)-(c) of Figure 9 plot the probability distribution functions (PDFs) of the shock Mach number, $N(M_{s})N_{j}/N_{t}$, for shock zones in the different regions of the jet for five fiducial models. Here, $N(M_{s})$ is the number of shock zones with $M_{s}$ in the range of [$\log M_{s}$,$\log M_{s}+d\log M_{s}$]. $N_{t}$ is the total number of grid zones in the volume encompassed by the bow shock surface, while $N_{j}=r_{j}/\Delta x$ is the number of grids across the jet radius. Note that $N(M_{s})$ is proportional to the area of shock surfaces, while $N_{j}$ and $N_{t}$ are proportional to the jet radius and the volume of the jet-induced structure; hence, $N(M_{s})N_{j}/N_{t}$ is effectively a dimensionless quantity.

The Mach number is the highest for the bow shock (cyan) and the next highest for shocks in the jet spine flow (red), and it is relatively low for shocks in the backflow (blue) and shocked ICM (green). But this is not necessarily the order of the shock speed (see below), since $v_{s}=M_{s}c_{s,1}$ depends on the preshock sound speed, $c_{s,1}$, as well. While $c_{s,1}\ll c$ in the background ICM, it is much larger in the jet spine flow and backflow; in particular, $c_{s,1}$ is close to $c/\sqrt{3}$ in the regions where the adiabatic index is close to $\gamma=4/3$ (see Figure 2).

Panels (a)-(c) of Figure 9 manifest two distinct populations of shocks: (1) the population for peaked PDFs with characteristic $M_{s}$’s, which consists of shock zones of the bow shock and recollimation shock surfaces, and (2) the population for power-law-like PDFs, which is associated with turbulent flows (e.g., Park & Ryu, 2019). The PDF for shock zones associated with the bow shock, peaks at $M_{s}\simeq 3.0$, 6.5, and 12.6 in the Q45, Q46, and Q47 models, respectively; the strength of the bow shock increases with increasing $Q_{j}$, as expected. The characteristic Mach number of shock zones associated with the first recollimation shock is $M_{s}\simeq 2.9$ and 3.0 in the Q46 and Q47 models. In the case of Q45-$\eta 5$-$\zeta 0$, the peak due to the recollimation shock is not noticeable, as it is expected to occur at smaller $M_{s}$. Again, the strength of the first recollimation shock increases with $Q_{j}$, which can be also guessed with the adiabatic index in the bottom panels of Figure 2. On the other hand, shocks in turbulent parts of the jet-induced structure are relatively weak with $M_{s}\lesssim$ a few. Shock zones associated with these disordered shocks follow fairly steep power-law-like distributions, $N(M_{s})\propto M_{s}^{-q}$; the slope ranges $q\sim 5-10$ for shocks of low $M_{s}$ in the jet spine flow, $q\sim 12-13$ in the backflow, and $q\sim 15-20$ in the shocked ICM. The averaged $M_{s}$ of these shocks is less dependent on $Q_{j}$.

The integration of the PDF for all the shocks (black lines) gives the total number of shock zones with $M_{s}\geq 1.01$, $N_{\rm shock}$, and the ratio, $N_{t}/(N_{\rm shock}N_{j})$, provides a measure for the mean distance between shock surfaces (in units of $r_{j}$) over the whole jet-induced structure. We point that $N_{\rm shock}$ and $N_{t}$ increase with time, but their ratio does not change much, once the jet-induced structures have more-or-less fully developed, for instance, at $t\gtrsim 20~{}t_{\rm cross}$ (see Figure 5). The mean distance measured at $t=30~{}t_{\rm cross}$ increases with $Q_{j}$ as $N_{t}/(N_{\rm shock}N_{j})=0.43$, 0.46, and 0.57 for Q45-$\eta 5$-$\zeta 0$, Q46-$\eta 5$-$\zeta 0$, and Q47-$\eta 5$-$\zeta 0$, respectively. This indicates that shocks are more frequent for smaller $Q_{j}$, consistent with the fact that more extended cocoons of backflow develop for lower power jets. The convergence of the PDFs in the default (thin lines) and high-resolution (thick lines) simulations for the Q46 and Q47 models, shown in panels (b) and (c), looks good. Yet, there are more shocks in higher resolution jets, $N_{t}/(N_{\rm shock}N_{j})=0.32$ and 0.38 for the Q46-$\eta 5$-$\zeta 0$-H and Q47-$\eta 5$-$\zeta 0$-H models, respectively. This confirms the previous statement that finer structures and more flow motions develop at smaller scales in higher resolution simulations.

Although the time evolution of the PDF of $M_{s}$ is not shown here, we find that shocks are somewhat stronger at earlier times, in particular, at $t\lesssim 10~{}t_{\rm cross}$. On the other hand, after $t\sim 20~{}t_{\rm cross}$, the overall distribution of $M_{s}$ does not change significantly over the time period of our simulations, whereas the jet propagation speed fluctuates somewhat (see Figure 5).

Panels (d)-(f) of Figure 9 show the PDFs of the shock speed, $N(\beta_{s})N_{j}/N_{t}$, for shock zones in the different regions of the jet, where $\beta_{s}=v_{s}/c$. The speed is the largest for shocks in the jet spine flow and the next largest for shocks in the backflow; shocks in the jet spine flow are relativistic with $\beta_{s}\sim 0.2-1.0$, while those in the backflow are mildly or sub-relativistic with $\beta_{s}\sim 0.01-0.4$. On the other hand, the bow shock and shocks in the shocked ICM are non-relativistic with characteristic values of $\beta_{s}\lesssim 0.005-0.05$, which increase with increasing $Q_{j}$.

To measure the importance of shocks in the jet flow dynamics, we quantify the shock dissipation by estimating the rate of the jet kinetic energy converted into heat at shocks. Panels (g)-(i) of Figure 9 show the energy dissipation rate at shock zones with the Mach number in the range of [$\log M_{s}$,$\log M_{s}+d\log M_{s}$],

which is normalized to the energy injection rate of the jet, $Q_{j}$. Here, $A_{\rm shock}=1.19(\Delta x)^{2}$ is the mean projected area of a shock zone within 3D space for random shock normal orientations (e.g., Ryu et al., 2003). The plots demonstrate that shocks in the jet spine flow and backflow are dominant in the shock dissipation, while the bow shock, even with high $M_{s}$’s, is relatively unimportant. Despite much lower preshock densities, shocks in the jet spine flow and backflow has higher $\Delta\Phi_{\rm shock}$ than those in the shocked ICM and the bow shock, because they have much higher $v_{s}$.

In Figure 10, we examine the integrated energy dissipation rate at shocks,

as a function of time. Here, the minimum value of $M_{s}$ is $M_{\rm min}=1.01$. Note that $\mathcal{L}_{\rm shock}/Q_{j}$ is the fraction of the jet-injected energy, which is dissipated into heat at shocks. Panel (a) shows the total dissipation rate, $\mathcal{L}_{\rm shock,tot}/Q_{j}$ due to all the shocks in the jet-induced flows for all the models considered here, while panels (b)-(d) show $\mathcal{L}_{\rm shock}/Q_{j}$ due to shocks in the different regions for the Q45, Q46, and Q47 models, respectively. As the jet penetrates into the background medium, a cocoon filled with shocks and turbulence develops. This leads to the initial increase of the shock dissipation. After the jet-induced structures have more-or-less fully developed, for instance, at $t\gtrsim 20~{}t_{\rm cross}$, $\mathcal{L}_{\rm shock}/Q_{j}$ gradually approaches asymptotic values, although it seems to increase somewhat even close to $t_{\rm end}$ in our simulations.

The following points are noticeable in Figure 10. (1) In panel (a), $\mathcal{L}_{\rm shock,tot}/Q_{j}$ is larger for smaller $Q_{j}$, as shocks are more frequent in lower power jets. However, still $\mathcal{L}_{\rm shock,tot}$ itself is larger for larger $Q_{j}$, since shocks on average have higher $v_{s}$ and dissipate larger amounts of energy in higher power jets. (2) Given the same $Q_{j}$, $\mathcal{L}_{\rm shock,tot}$ is somewhat smaller for the models with higher $\rho_{j}$ (smaller $\eta$). The increase of the jet pressure by an order of magnitude ($\zeta=10$) has only marginal effects on $\mathcal{L}_{\rm shock,tot}$. (3) In panels (b)-(d), the shock dissipation occurs mostly at shocks in the backflow (blue) and the jet spine flow (red), while the dissipation at shocks in the shocked ICM (green) and the bow shock (cyan) is much less. Thus, regarding the shock dissipation, the backflow is most important in the Q45 and Q46 models, while the backflow and the jet spine flow are about equally important in the Q47 models. (4) The comparison of $\mathcal{L}_{\rm shock}/Q_{j}$ for high-resolution models (bold-solid lines) with that for the corresponding default-resolution models (solid lines) demonstrates that the shock dissipation fraction is higher owing to more frequent shocks at higher resolution simulations.

As noted in Section 3.2, the kinetic energy is not strictly defined in RHDs. We assume the form in Equation (14), which reduces to the Newtonian kinetic energy in the non-relativistic limit. In addition, the kinetic energy dissipation rate at shocks is a frame-dependent quantity, and we employ $\Delta\Phi_{\rm shock}$ in Equation (18), which is valid in the shock-rest frame. Thus, our estimations of $\mathcal{L}_{\rm shock,tot}$ should be considered to be only approximate. Nevertheless, $\mathcal{L}_{\rm shock,tot}/Q_{j}\sim 0.5-1$, $0.45-0.8$, and $0.1-0.15$ at $t=t_{\rm end}$ for the Q45, Q46, and Q47 models, respectively. These results indicate that a substantial fraction of the jet energy is dissipated at shocks in the jet-induced flows; shocks in the backflow is important for the shock dissipation, even though their speeds are mildly relativistic.

In diffusive shock acceleration (DSA) theory, cosmic rays (CRs) are scattered off MHD waves in the shock converging flows and gain energy by crossing repeatedly back and forth across the shock front (Bell, 1978). In non-relativistic shocks, it predicts a power-law distribution of CRs, $f(p)\propto p^{-q}$, with slope, $q=4M_{s}^{2}/(M_{s}^{2}-1)$, in the test-particle regime (Drury, 1983). The slope becomes $q=4$ for strong non-relativistic shocks with large $M_{s}$. In relativistic shocks, on the other hand, a power-low spectrum with a steeper slope, $q\sim 4.2-4.3$, is predicted (e.g., Achterberg et al., 2001), and the acceleration becomes less efficient at the ultra-relativistic limit due to the anisotropic particle distribution and the limited residence time at both the upstream and downstream regions (e.g., Bykov et al., 2012; Sironi et al., 2015). Hence, among shocks in the jet-induced flows, we expect that mildly relativistic shocks with $v_{s}/c\lesssim 0.1$ in the backflow could play important roles in producing UHECRs, as previously suggested (Bell et al., 2018; Matthews et al., 2019).

On the other hand, as shown above, shocks are quite frequent with the mean distance between shock surfaces, $\langle l_{\rm sh}\rangle\sim 0.3-0.6~{}r_{j}$, for shocks with $M_{s}\geq 1.01$ in the jet-induced flows. Considering $r_{j}=1$ kpc in our models, UHECRs with $E_{\rm CR}\gtrsim 1$ EeV may have, $r_{g}\approx 1.1~{}{\rm kpc}(E_{\rm CR}/10^{18}~{}{\rm eV})(B/1\mu{\rm G})^{-1}>\langle l_{\rm sh}\rangle$, that is, the gyroradius larger than the mean distance between shock surfaces. Hence, UHECRs may encounter more than one shocks in a scattering length and get accelerated differently, compared to a single episode of DSA. This rather complex scenario involving multiple shocks needs to be further examined.

Although not shown as plots, the comparison of the higher $\rho_{j}$ models ($\eta 3$ and $\eta 4$) with the fiducial models ($\eta 5$) indicates that for the same $Q_{j}$, $N(M_{s})$ for the first recollimation shock including the peak shifts to higher $M_{s}$ with higher $\rho_{j}$, since $p_{j}/\rho_{j}$ is lower and hence $c_{s,1}$ is smaller in the injected jet material. In addition, $N(\beta_{s})$ for shocks in the jet spine flow and backflow shifts to lower $\beta_{s}$ with higher $\rho_{j}$ (smaller $v_{j}$ and $\Gamma_{j}$). Comparing the $\eta 4$-$\zeta 1$ models with the $\eta 4$-$\zeta 0$ models, in the over-pressured jets, $N(M_{s})$ for the first recollimation shock shifts to lower $M_{s}$, owing to higher $c_{s,1}$ in the injected material. Otherwise, other characteristics of the PDFs seem to be less affected by different $\eta$ and $\zeta$.

As shown in the left half of each panel in Figure 3, a relativistic velocity shear with $\Delta v_{z}\sim c$ develops at the interface between the upward-moving jet spine flow (red) and the downward-moving backflow (blue). In addition, a non-relativistic velocity shear with $\Delta v_{z}<0.1c$ develops at the interface between the backflow and the shocked ICM (white). These shear interfaces are unstable against the KH instability, so turbulence emerges in the jet spine flow and backflow as well as in the shocked ICM.

Assuming the presence of turbulent magnetic fields, UHECRs may be accelerated through elastic scatterings off turbulent waves moving with the underlying shear flow, especially around the jet-backflow boundary (see, e.g., Rieger, 2019, for a review). In the so-called discrete shear acceleration, which operates when the particle scattering length is larger than the width of the velocity shear layer, the mean energy gain is given as $\Delta E/E\sim\Gamma_{j}-1$ for each crossing of the jet-backflow boundary, if the CR distribution is nearly isotropic around the boundary (e.g., Ostrowski, 1998). On the other hand, when the particle scattering length is smaller than the velocity shear scale, the so-called gradual shear acceleration operates and the energy gain is $\Delta E/E\propto(\bar{v}/c)^{2}$ (Fermi-II), where $\bar{v}=(\partial v_{z}/\partial r)\lambda$ is the effective velocity difference that the particles with the mean free path, $\lambda$, experience in the shear flow with $\mbox{\boldmath$v$}=v_{z}(r)\hat{\mbox{\boldmath$z$}}$ (Rieger, 2019). With the relativistic gradual shear of the jet, the acceleration timescale is inversely proportional to the relativistic shear coefficient defined as,

where $\Gamma_{z}=1/\sqrt{1-(v_{z}/c)^{2}}$ (Webb et al., 2018).

Figure 11 shows the 2D slice images of the magnitude of the velocity shear, $\Omega_{\rm shear}\equiv|{\partial v_{z}}/{\partial r}|$ (in the left half of each panel, $x<0$), and the relativistic shear coefficient, $\mathcal{S}_{r}$ (in the right half of each panel, $x>0$), for three fiducial models. Figure 12 plots their PDFs, $N(\Omega_{\rm shear})$ (top panels) and $N(\mathcal{S}_{r})$ (bottom panels), in the different regions for five fiducial models. The value of $\mathcal{S}_{r}$ in the shocked ICM is very small, so its PDF is not shown. While the velocity shear is somewhat larger for higher power jets, the dependence on $Q_{j}$ is not strong. The peaks of $N(\Omega_{\rm shear})$ lie at $\Omega_{\rm shear}(r_{j}/c)\gtrsim 1$ for the jet spine flow, $\sim 0.1-1$ for the backflow, and $\lesssim 0.01$ for the shocked ICM, respectively. On the other hand, $\mathcal{S}_{r}$ in the jet spine flow is noticeably larger for higher power jets, as can be seen with the purple tone shade in Figure 11, owing to the weighting factor of $\Gamma_{z}^{4}$. Figure 12 shows that it extends up to $\mathcal{S}_{r}(r_{j}/c)^{2}\sim 10^{3}-10^{5}$ inside the jet spine flow (red dashed lines), depending on $Q_{j}$. In the backflow, $\mathcal{S}_{r}$ depends only weakly on $Q_{j}$, and $N(\mathcal{S}_{r})$ peaks at $\mathcal{S}_{r}(r_{j}/c)^{2}\sim 10^{-3}-10^{-2}$. Although not shown here, $\Omega_{\rm shear}$ and $\mathcal{S}_{r}$ are larger in the earlier stage, but their PDFs converge as the jet-induced structures approach more or less steady states at $t\gtrsim 20~{}t_{\rm cross}$. Again although not shown here, the models with the same $Q_{j}$ but different $\eta$ and $\zeta$ have similar $N(\Omega_{\rm shear})$ and $N(\mathcal{S}_{r})$, except that $\mathcal{S}_{r}$ is smaller in the models with higher $\rho_{j}$ (smaller $\Gamma_{j}$).

While both $|{\partial v_{z}}/{\partial r}|$ and $\mathcal{S}_{r}$ are the largest in the jet spine flow, the cocoon of the backflow occupies a much larger volume than the jet spine. Hence, we expect that lower energy CRs would be accelerated via the gradual shear acceleration mostly in the backflow (e.g., Rieger & Duffy, 2004). In comparison, higher energy CRs with the gyroradius $\gtrsim 1$ kpc may diffuse across the jet-backflow boundary and gain energy via the discrete shear acceleration (e.g., Kimura et al., 2018). However, the exact process of the shear acceleration will depend on the details such as the jet geometry, the origin of seed particles, the magnetic field strength and configuration, the magnetic fluctuations that scatter CRs, and etc.

As discussed above, turbulence develops in the jet-induced flows as a consequence of the KH instability in the shear boundaries and also the ensuing complex flow dynamics. Then, along with the shear acceleration, the stochastic Fermi-II acceleration due to elastic scatterings off turbulent magnetic fluctuations should operate, especially inside the cocoon (e.g., O’Sullivan et al., 2009; Hardcastle, 2010). As a measure of the turbulence, we define the vorticity excluding the shear considered in the previous subsection as

where $\mbox{\boldmath$\Omega$}_{t}=\mbox{\boldmath$\nabla$}\times\mbox{\boldmath$v$}$ is the total vorticity.

Figure 13 shows the 2D slice images of $\Omega_{t}=|\mbox{\boldmath$\Omega$}_{t}|$ (in the left half of each panel, $x<0$) and $\Omega_{-}=|\mbox{\boldmath$\Omega$}_{-}|$ (in the right half of each panel, $x>0$) for three fiducial models. The spatial distributions of $\Omega_{t}$ and $\Omega_{-}$ are similar to those of the velocity shear, $\Omega_{\rm shear}$, shown in Figure 11. They look more or less homogeneous in the cocoon, as expected in fully developed turbulence, but patterns appear in the shocked ICM. Some of the patterns follow the shocks shown in Figure 8, and others match the herringbone patterns in the density distributions in Figure 3. While those patterns look interesting, overall the flow dynamics in the shocked ICM is not very important in terms of the energy dissipations, as shown in Figure 10 and also in Figure 16 below.

Figure 14 plots the PDFs of $\Omega_{t}$ and $\Omega_{-}$, $N(\Omega_{t})$ (top panels) and $N(\Omega_{-})$ (bottom panels), for five fiducial models. They behave similarly as $N(\Omega_{\rm shear})$ shown in Figure 12; the mean values of $\Omega$’s are the greatest in the jet spine flow (red lines) and the smallest in the shocked ICM (green lines). Moreover, they show only weak dependence on the jet power, $Q_{j}$, although the vorticities increase slightly with increasing $Q_{j}$. We note that on average $\Omega_{-}$ is comparable to, or somewhat larger than, $\Omega_{\rm shear}$, and $N(\Omega_{-})$ peaks at slightly higher values than $N(\Omega_{\rm shear})$. This indicates that turbulence develops fully in the shear boundaries, and hence the radial and vertical components of $\mbox{\boldmath$\Omega$}_{t}$ is comparable to $\Omega_{\rm shear}$. As in the case of $\Omega_{\rm shear}$, both $\Omega_{t}$ and $\Omega_{-}$ are larger in the earlier stage, but their PDFs converge at $t\gtrsim 20~{}t_{\rm cross}$. Furthermore, the models with the same $Q_{j}$ but different $\eta$ and $\zeta$ have similar PDFs.

Our finding, $\Omega_{-}\gtrsim\Omega_{\rm shear}$, implies that the stochastic turbulent acceleration of UHECRs could be as important as the gradual shear acceleration for CRs with the diffusion length smaller than the width of the velocity shear layers of the jet. However, the real situation could be quite complex, since the two acceleration processes may operate simultaneously and their interplay would occur. Then, they may need to be examined as the turbulent shear acceleration (e.g., Ohira, 2013).

The turbulence developed in the jet-induced flows cascades down to smaller scales and is eventually converted into heat. Here, we examine this process, known as the turbulent dissipation, to measure the importance of turbulence in the jet flow dynamics and also to assess the relative importance of the turbulent versus shock dissipations in relativistic jets.

Despite the relevance to many astrophysical phenomena such as relativistic jets and accretion disks, turbulence in relativistic regime has been less explored than its non-relativistic counterpart. Recently, however, systematic studies of RHD and RMHD turbulence using high-resolution simulations have been performed (e.g., Zrake & MacFadyen, 2012, 2013; Radice & Rezzolla, 2013). In particular, the RHD studies showed that the velocity power spectrum follows the $k^{-5/3}$ power-law and the kinetic energy dissipation scaled as $v_{\ell}^{3}/\ell$ is independent of $\ell$, in good agreement with the predictions of the classical Kolmogorov theory, at least, for mildly relativistic, hydrodynamic flows with the mean Lorentz factor of $\lesssim$ a few (see, e.g., Zrake & MacFadyen, 2013; Radice & Rezzolla, 2013). Here, $k$ is the wavenumber and $v_{\ell}$ is the characteristic speed of turbulent flows at the length scale, $\ell$.

Our simulated jets contain the bulk motions of the upward-moving jet spine flow and the downward-moving backflow. Hence, the flow motions associated with turbulence need to be separated from the large-scale bulk motions. The extraction of turbulent flow velocity could be done by applying a “filtering”, and we employ the scheme used in Vazza et al. (2017). The mean velocity is computed as

where the summation runs over the cubic box of the size $\Lambda$ around the position $r$. Here, $w_{i}$ is a weighting function, and we set $w_{i}=1$. Then, the turbulent velocity is estimated as $\mbox{\boldmath$v$}_{\rm turb}(\mbox{\boldmath$r$})=\mbox{\boldmath$v$}(\mbox{\boldmath$r$})-\left<\mbox{\boldmath$v$}(\mbox{\boldmath$r$})\right>_{\Lambda}$. The estimated $\mbox{\boldmath$v$}_{\rm turb}$ should depend on the filtering size, $\Lambda$. Vazza et al. (2017) suggested that the mean size of eddies in turbulent flow motions would be a good choice for $\Lambda$. Here, we set $\Lambda=r_{j}$, since the jet radius represents a characteristic scale in the jet-induced flows.

Figure 15 shows the power spectrum, $P_{v}$, of $\mbox{\boldmath$v$}_{\rm turb}$ in the different regions for three fiducial models. A few points are evident. (1) The peaks of $P_{v}$ occur around $k\sim 1$, corresponding to the wavelength $\lambda\sim r_{j}$, due to the filtering scale we impose. (2) At $k\gtrsim 1$, i.e., at smaller scales of $\lambda<r_{j}$, $P_{v}$ exhibits the Kolmogorov power-law of $k^{-5/3}$, despite the fact that there are numerous shocks in the flows. If the velocity spectrum is dominated by shocks in non-relativistic regime, $P_{v}$ follows $k^{-2}$, the Burger’s power spectrum (e.g., Kim & Ryu, 2005). The Kolmogorov spectrum here should reflect the characteristics of RHD turbulence (Zrake & MacFadyen, 2013; Radice & Rezzolla, 2013). In addition, the contribution of shocks to $P_{v}$ may not be substantial, because most shocks are weak with $M_{s}\lesssim$ a few. (3) $P_{v}$ is the largest in the jet spine flow, the next largest in the backflow, and the smallest in the shocked ICM, revealing the rms magnitude of turbulent flow motions. In addition, although not shown here, the solenoidal mode, $P_{\rm sol}$, dominates over the compressive mode, $P_{\rm comp}$, in the jet spine flow and backflow; $P_{\rm sol}$ is about an order of magnitude larger than $P_{\rm comp}$ there. On the other hand, in the shocked ICM, while $P_{\rm sol}$ is still larger than $P_{\rm comp}$, the difference of the two is small.

We then estimate the energy dissipation rate due to turbulent flow motions, $\mathcal{L}_{\rm turb}$, using the following approximate quantity:

where $\varphi$ is a free parameter of order of unity and to be set as $\varphi=1$. Here, $E_{\rm turb}\approx\Gamma_{\rm turb}\left(\Gamma_{\rm turb}-1\right)\rho c^{2}$ is the kinetic energy of turbulent flows, where $\Gamma_{\rm turb}=1/\sqrt{1-(v_{\rm turb}/c)^{2}}$. And $\tau_{\rm cascade}\approx\lambda_{\rm peak}/v_{\rm rms}$ is the timescale of turbulent cascade, where $v_{\rm rms}$ is the rms of $\mbox{\boldmath$v$}_{\rm turb}$ and $\lambda_{\rm peak}\sim r_{j}$ is the peak scale of $P_{v}$. Then, $\mathcal{L}_{\rm turb}/Q_{j}$ can be regarded as the fraction of the jet-injected energy, dissipated into heat through turbulent cascade.