Issue 
A&A
Volume 679, November 2023



Article Number  A160  
Number of page(s)  19  
Section  Extragalactic astronomy  
DOI  https://doi.org/10.1051/00046361/202244690  
Published online  30 November 2023 
Simulations of twotemperature jets in galaxy clusters
I. Effect of jet magnetization on dynamics and electron heating
^{1}
Institute for Cosmic Ray Research, The University of Tokyo, 515 Kashiwanoha, Kashiwa, Chiba 2778582, Japan
email: tohmura@icrr.utokyo.ac.jp
^{2}
Department of Physics, Faculty of Sciences, Kyushu University, 744 Motooka, Nishiku, Fukuoka 8190395, Japan
^{3}
Division of Science, National Astronomical Observatory of Japan, 2211 Osawa, Mitaka, Tokyo 1818588, Japan
email: mami.machida@nao.ac.jp
^{4}
Astronomical Science Program, The Graduate University for Advanced Studies, SOKENDAI, 2211 Osawa, Mitaka, Tokyo 1818588, Japan
^{5}
Department of Astrofusion Plasma Physics (AFP), Headquarters for CoCreation Strategy, National Institutes of Natural Sciences, Minatoku, Tokyo 1050001, Japan
Received:
5
August
2022
Accepted:
9
September
2023
Context. Nonradiating protons in the radio lobes play an essential role in shaping the jet morphology, as demonstrated in recent radio and Xray observations. However, since protons and electrons are not always in energy equilibrium due to weak Coulomb coupling, it is difficult to estimate the energy contribution of protons for the inflation of radio lobes.
Aims. The focus of this study is to examine the effect of the variable model for electron heating by turbulence and shock waves on the thermal energy distribution of electrons and protons.
Methods. We performed twotemperature threedimensional magnetohydrodynamic (3D MHD) simulations of subrelativistic jets in the galaxy cluster, while varying the jet magnetization parameters. Because the energy partition rate between electrons and protons in shock and turbulence is determined by plasma kinetic scale physics, our global simulations include electron instantaneous heating subgrid models for shock waves and turbulence.
Results. We find that most of the bulk kinetic energy of the jet is converted into the thermal energy of protons through both shocks and turbulence. Thus, protons are energetically dominant. Meanwhile, thermal electrons stored in the lobe evolve toward energy equipartition with magnetic energy through turbulent dissipation. We further estimated the radio power and the mechanical jet power of radio lobes following the same method used for radio and Xray observations, then we compared these powers with that of the observed radio jets. The twotemperature model quantitatively explains the radiatively inefficient radio cavities, but it cannot reproduce the radiatively efficient cavity, even for strongly magnetized jets. This implies that a significant population of pairplasma is needed to explain radiatively efficient radio cavities.
Key words: galaxies: jets / magnetohydrodynamics (MHD) / radio continuum: galaxies
© The Authors 2023
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
This article is published in open access under the Subscribe to Open model. Subscribe to A&A to support open access publication.
1. Introduction
Jets that are driven by the active galactic nuclei play a significant role in galaxy and cluster evolution. They propagate beyond the spatial scale of its host galaxy and transport its own kinetic energy into the surrounding intracluster medium (ICM). This heating energy prevents catastrophic cooling within the ICM and keeps it from falling into the host galaxy (Fabian 1994, 2012; McNamara & Nulsen 2007). To better understand and discuss the heating energy of ICM in a quantitative manner, estimating the total kinetic energy of the jet is of great importance. Relativistic electrons and magnetic fields contained in the jetdriven lobes can be constrained by radio observations using the equipartition energy condition for the energy of cosmic ray and magnetic fields (Perley et al. 1984; Beck & Krause 2005). It is, however, difficult to estimate the total kinetic power of the jet, as most of the radiation stems from nonthermal electron origin (Hardcastle & Croston 2020), that is, the energy of thermal electrons and protons cannot arise directly from the source signals.
In the context of jets in galaxy clusters, the jet kinetic energy is constrained by observation of the Xray cavity, which is produced as a result of the radio lobe inflation (Fabian et al. 2000; McNamara et al. 2001). The minimum energy needed to form a cavity has been estimated as E = P_{cav}t_{age} = 4 pV, where P_{cav}, t_{age}, p, and V are the cavity power, source age, and pressure of the ICM observed by thermal Xray, as well as the volume of the cavity, respectively. Bîrzan et al. (2008, 2004) found that the mechanical power estimated from the Xray cavity seem to be correlated with the radio luminosity (sum of core and lobes) and that the cavity increases with radio luminosity, , where 0.35 ≤ α ≤ 0.70. The median ratio of the mechanical power to radio luminosity (radiative efficiency) is P_{cav}/P_{Radio} ∼ 100. These results imply that the energy contribution of nonradiating protons is needed for cavity formation. However, there is a large scatter in this relation. For example, Cygnus A is the most radiatively efficient system, with P_{cav}/P_{Radio} ∼ 1. Although there are several physical factors of these scatters, such as electron cooling, estimation of ages, plasma composition, and the variable activity of active galactic nuclei (AGNs), their contribution in creating this scatter is not yet understood.
Fluid simulations can follow highly complex jet flows and are a useful approach to examine the energy transport in the jetICM system. To approach the realistic conditions of radio lobes, some numerical models implement key physics, such as magnetic fields (Massaglia et al. 2019), special relativity (Aloy et al. 1999), difference plasma compositions (Scheck et al. 2002; Perucho et al. 2014), cosmic ray electrons (Jones et al. 1999; Mendygral et al. 2012; Mukherjee et al. 2021), and cosmic ray protons (Mathews & Guo 2010; Weinberger et al. 2017). Several studies have pointed out that magnetohydrodynamic (MHD) instabilities develop nonaxisymmetric modes, namely Kelvin–Helmholtz (Bodo et al. 1994), Rayleigh–Taylor (Matsumoto & Masada 2013), and currentdriven kink modes (Mizuno et al. 2009; Mignone et al. 2010; Porth & Komissarov 2015). Thus, the nonlinear evolution of these instabilities plays a significant role in the jet dynamics and their largescale morphology (e.g., Tchekhovskoy & Bromberg 2016). In particular, the development of instabilities that cause the jet deceleration and/or jet disruption, could directly link to the physical reasons of the Fanaroff–Riley (FR) distinction (Fanaroff & Riley 1974).
Our previous studies (Ohmura et al. 2019, 2020) have focused on the twotemperature plasma, where the electrons and protons in jets are not in thermal equilibrium, since Coulomb coupling is inefficient in the tenuous plasma (Braginskii 1965; Stepney & Guilbert 1983). Here, we discuss the results of axisymmetric simulations with a constant fraction model for electron and proton heating. Because thermal protons heated up at the internal shocks, and as the Coulomb coupling did not work effectively in the jets, the proton temperature was several times higher than electron temperature. Thus, we found that thermal protons support the expansion of the cocoon, rather than thermal electrons.
Recent theoretical studies clarified the physical picture of electron heating in the collisionless turbulence (Howes 2010; Kawazura et al. 2019, 2020). These provide the variable model that describes the partition of heating energy between protons and electrons for the dissipation at the plasma kinetic scale. The heating fraction between electrons and protons in this model is not constant and represents an increasing function of the proton plasma, β, β_{p} ≡ 8πnkT_{p}/B^{2}. However, in our previous simulation, we did not focus on the electron heating in turbulent conditions. Our previous results indicated that the magnetic fields accumulate and the proton plasmaβ decreased in the cocoon. Furthermore, the magnetic field can be amplified locally by nonaxisymmetric motion. Therefore, the turbulence is the dominant heating source for thermal electrons, compared with shock heating.
Jet magnetization parameters, namely, jet Alfvén Mach number ℳ_{A} and plasmaββ_{gas}, are important for both electron heating and dynamical evolution. However, none of the studies considered twotemperature plasma, thus they’ve failed to explore the effect of different jet magnetization on the electron heating with the development of MHD instabilities. The purpose of this study is to examine the effect of the variable model for electron heating under turbulence on the distribution of electron temperature while varying jet magnetization parameters.
In this study, we present the results of twotemperature threedimensional (3D) MHD simulations of semirelativistic jets in galaxy clusters. In Sect. 2, basic equations and subgrid models for electron heating are presented. We describe the setup of our simulation in Sect. 3 and the results are presented in Sect. 4. In Sect. 5, we discuss the observational implications of this study. A brief summary of key findings is provided in Sect. 6. The appendix describes detailed numerical methods for solving entropy equations and the shockfinding algorithm.
2. Numerical method
2.1. Basic equations
Multiple methods have been developed to incorporate electron thermodynamics into singlefluid simulations selfconsistently in the context of the general relativistic MHD simulations for hot accretion flow (Ressler et al. 2015; Sadowski et al. 2017). In this work, we extend the singletemperature MHD code CANS+ (Matsumoto et al. 2019) to a twotemperature framework as follow the method in Sadowski et al. (2017). The total gas (summed electrons and protons) evolves by the MHD equations in the following conservation form:
where U, F, and S are the vector of conserved quantities, the vectors of flux, and the vector of source term, respectively. The conserved quantities are expressed as:
where ρ, v and B are the mass density, the bulk velocity, and the magnetic field, respectively. We assume the gas to be hydrogen, such that n_{e} ∼ n_{p}, where n_{e, p} represents the number density of protons and electrons, respectively. Then, ρ = m_{e}n_{e} + m_{p}n_{p} ≈ m_{p}n. The total energy, E, and total pressure, p_{T}, are:
where p_{gas} = p_{p} + p_{e} is the gas pressure, which sums the proton pressure, p_{p}, and electron pressure, p_{e}. The fluxes are then:
The source terms are
q_{rad} is the radiative energy loss rate. In this study, we assume the radiation process as bremsstrahlung emission. Notably, an adiabatic index of gas γ_{gas} and the radiative energy loss rate by bremsstrahlung radiation, q_{rad}, are functions of the electron and proton temperature. The primitive variables of the above system of equations are:
In addition to solving MHD equations, we solve the entropy equations of the two species to obtain each temperature. The entropy equations of electrons and protons can be expressed as:
where q_{ie}, f_{e}, and q_{heat} are the energy transfer ratio via Coulomb coupling, the fraction of electron heating, and the dissipation heating rate, respectively. The detailed procedures of numerical integration are described in Appendix A.
We address the transrelativistic regime for electrons, and use the following approximate entropy formula derived by Sadowski et al. (2017) in Appendix A:
where θ_{e} ≡ kT_{e}/m_{e}c^{2} is the dimensionless temperature. In this approximate formula, we can easily obtain the electron temperature for given s_{e} and ρ_{e} as:
The adiabatic index for electrons is calculated as follows:
In contrast, protons are nonrelativistic in this simulation. Thus, we use the nonrelativistic entropy formula for protons:
where γ_{p} = 5/3. The thermal energies of protons, electrons, and gas are as follows:
From the relationship between the gas pressure and gas thermal energy, the effective adiabatic index for the gas can be calculated as (Ressler et al. 2015):
2.2. Subgrid models of electron heating
We considered two subgrid models for the fraction of electron heating, f_{e}. One model represents turbulence heating, f_{e,turb}, and another model represents the shock heating, f_{e,shock}. Therefore, f_{e} is determined by plasma properties at each simulation grid. First, we identify a shock zone by shockfinding method based on Ryu et al. (2003) and Schaal & Springel (2015; see detail in Appendix B). The fraction of shock heating is adopted only in the shock zone, and the other region is adopted by the fraction of turbulence heating, i.e.,
We note that part of the dissipation energy is small in the region of laminar flow and, hence, the heating fraction of turbulence heating spontaneously works in the turbulence zone.
For the shock zone, we modeled a constant electron heating fraction, f_{e,shock} = 0.05. This value is justified by the observation data in the solar system and supernova remnants shocks (Vink et al. 2015). Furthermore, some theoretical simulations, based on particleincell (PIC) simulation, indicate that electrons irreversibly heat up, while protons are primarily heated during collisionless shocks (Matsukiyo 2010; Guo et al. 2018; Crumley et al. 2019; Tran & Sironi 2020). The validity of this parameter is previously discussed in Sect. 4.1 of Ohmura et al. (2020).
For electron to proton heating rates of MHD turbulence, there are two models, which were proposed by Howes (2010, hereafter H10) and Kawazura et al. (2019, hereafter K19). The model comparison between H10 and K19 is shown in Appendix C. Both the models are based on the gyrokinetics approach to damping weakly collisional MHD turbulence. H10 provided the heating model derived from the linear theory for the first time, while K19 treated the nonliner evolution of the turbulence by numerical simulations. Thus, K19 would be favored. We use the K19 for the unshocked zone:
where P_{compr} and P_{AW} are the compressive energy injection and Alfvénic energy injection, respectively. Therefore, the fraction of electron heating, f_{e}, is defined by:
It is a difficult task to estimate the ratio of the compressive energy injection and the Alfvénic energy injection in fluid simulations. Therefore, we assume pure Alfvénic turbulence (i.e., P_{compr}/P_{AW} → 0). Notably, this assumption leads to an overestimation of the amount of electron heating. This heating model represents a weak dependence of the temperature ratio T_{e}/T_{p}, but a strong dependence of proton plasma beta, β_{p}. In the case of β_{p} < 1, most of the dissipation energy is absorbed by electrons and vice versa.
3. Simulation setup
We carried out the twotemperature 3D MHD simulations in Cartesian coordinates with the zaxis pointing along the jet direction. The computational domain is x ∈ ( − L_{x}/2, L_{x}/2), y ∈ ( − L_{y}/2, L_{y}/2), and z ∈ (0, L_{z}), where L_{x}, L_{y}, and L_{z} denote the length of the computational domain. We use a uniform mesh of (N_{x}, N_{y}, N_{z}) with size Δ_{x} = Δ_{y} = Δ_{z} = 0.1 kpc. The grid number and length of the computational domain are given in Table 1. We permit the backflow to escape from the boundary at z = 0. Therefore, the absorbing boundary condition applies to the xzplane at z = 0. Other boundaries are imposed on the freeboundary condition.
Numerical models.
3.1. Initial condition
To study the interaction between jets and the ICM, we initialize the surrounding ICM in the form of β profile (King 1962). Our cluster model is roughly consistent with the environment of Cygnus A from Chandra Xray data (Wilson et al. 2006; Smith et al. 2002). The density profile of ICM is given by:
where , n_{0}, r_{c} and β′ are the radius, core density, core radius, and ratio of the specific energy in galaxies to the specific thermal energy in the ICM, respectively. We set β′ = 0.5, r_{c} = 20 kpc, and n_{0} = 0.05 cm^{−3}. We also assume that our atmosphere is initially isothermal with the temperature kT_{p} = kT_{e} = 5 keV. We employ the uniform magnetic field B_{ICM} that is parallel to the zaxis, and B_{ICM} = 0.44 μG. The blue lines of Fig. 1 show the initial density (top panel) and pressure (bottom panel), respectively.
Fig. 1. Number density (top) and gas pressure (bottom) profiles of initial ICM as a function of radius. Blue dots represent the jets number density and the jets gas pressure, respectively. The initial gas temperature of ICM is 5 KeV over the entire simulation domain. 
3.2. Jet model
Focusing on the effect of jet magnetization on electron heating and jet stability, we carried out simulations with different magnetic field strengths. We modeled magnetized supersonic subrelativistic flows to be consistent with the outburst energy of Cygnus A jets, 0.6 − 0.8 × 10^{46} erg s^{−1} (Bîrzan et al. 2004; Snios et al. 2018). The injected jet power can be written as follows:
We set the kinetic power as L_{kin} = 5.5 × 10^{45} erg s^{−1}, and the thermal power as L_{th} = 4.4 × 10^{44} erg s^{−1}.
To generate the jet beams, we injected supersonic and magnetized flows inside a constant cylindrical nozzle at the origin. The radius and length of the nozzle are 1 kpc and 1.2 kpc, respectively. Although it is difficult to determine the real flow radius of jets by observation. Radio observations indicate that the beam radius of Cyguns A is 0.1 to 1.0 kpc at 1 kpc from the central engine (Nakahara et al. 2019). We assume that the jet temperature and the velocity are T_{p} = T_{e} = 10^{10} K and v_{jet} = 0.3c, respectively. Thus, the internal sonic Mach number is ℳ = 6.2. Our jet models satisfy the condition that the thermal pressure ratio between the jet and the ICM in the launching region is unity. In Fig. 1, the blue dots for each panel denote the density and pressure in the jet injection region, respectively. A smallamplitude (1%) random pressure perturbation for the injection flow is adapted to model nonaxisymmetric features. We list common parameters for ICM and jets in Table 2.
Jets and ICM common setup parameters.
The jets have a purely toroidal magnetic field B_{ϕ} = B_{jet}sin^{4}(πr/r_{jet}). As shown in Table 1, Models A, B, and C have different values of gas plasma β_{gas} equal to 1, 5, 100, respectively. In these case, the amplitudes of the injected magnetic field are B_{jet} = 138, 62, 14 μG for models A, B, and C, respectively. Our models are matterdominated jets, namely, therein, the kinetic energy of jets exceeds the Poynting flux energy.
4. Results
We conducted simulations with various magnetic energies to investigate the effect of jet dynamics and electron heating. First, we focused on the effect of the magnetic field strength on the jet dynamics. The strength of the magnetic field affects the development of instabilities, such as kink, Kelvin–Helmholtz, or Rayleigh–Taylor modes. The electron heating model for turbulence is the function of plasmaβ_{p}, hence the magnetization strongly affects the electron temperature distributions. Here, we report the time evolution for thermal electrons and protons in the jet lobe.
4.1. Overall morphology and beam stability
Figure 2 shows the density slice in the yzplane of x = 0 kpc at the end of the simulation (t = 9.52, 9.94, and 13.02 Myr) for models A, B, and C. The shockedICM, compressed by the bow shock, and the lowdensity cocoon are formed. We observe that the number density of the cocoons is very low, ∼10^{−4} cm^{−3}, such that Coulomb coupling is inefficient. The jet beam reaches jet tip despite suffering MHD instabilities (see red contours in Fig. 2), and a terminate shock is formed at the end of the jet for all models. We find that the nonaxisymmetric mode is developed in model A and B, as the shapes of the bow shock are affected by the bending motion of the jets for both models. Strong pressure waves are generated at the termination of the beam, which push up the bow shock.
Fig. 2. Slices (in the y − z plane) of number density distribution for model A, B, and C at t = 9.52, 9.94, and 13.02 Myr, respectively. The white lines represent contours of the zcomponent of velocity v_{z} = 0.5v_{jet} (0.15c). 
For models A and B, the jet develops a pronounced helical shape that is characteristic for kink instabilities. Although our jet has a purely toroidal magnetic field at the launching region, a helical component is generated during jet propagation. The timescale of the development of an external kink mode corresponds to the Alfvén crossing time in the beam (Moll et al. 2008; Mizuno et al. 2009),
where v_{A, ϕ} is the azimuthal Alfvén velocity. A fluid element in the beam has roughly a constant velocity, denoted v_{jet}. We verify that the bulk beam velocity changes slightly, but roughly maintains an injection velocity. Therefore, the kink mode develops after the jets propagate to the distance as l_{kink} ∼ v_{jet}τ_{kink}. For models A and B, the distances of l_{kink} ∼ 40, 70 kpc are, respectively, within the simulation domain. Meanwhile, l_{kink} ∼ 300 kpc is larger than the length of the simulation domain, hence model C jet is not expected to develop the kink mode in simulation time.
Figure 3 shows the twodimensional distribution maps of the beam barycenter R_{fG}, which is described as the distance from the origin, as function of time for each model. The large value of the barycenter indicates development of a nonaxisymmetric mode. We compute the beam barycenter as follows:
Fig. 3. Twodimensional distribution maps of beam barycenter R_{fG} as a function of time for model A (top), B (middle), and C (bottom), respectively. Horizontal yellow solid lines depict the distance, l_{kink}, for the development of the external kink model, and vertical yellow solid lines are the time, τ_{kink}, at which the jets propagate to distance, l_{kink}. 
A similar analysis was performed by Mignone et al. (2013). We see that l_{kink} is a good indicator of the kink instability for both models A and B. After the jet propagates to the distance l_{kink}, the barycenter is larger than 3 kpc for both models. For model B, the time at which the jet propagates to the distance l_{kink} is 7.5 Myr, which corresponds to the time when it starts to decelerate. Thus, the nonaxisymmetric mode developed by the kink instability induces deceleration by increasing the size of the jet head. Meanwhile, because the jet does not suffer the kink mode, the barycenter for model C is within 1 kpc in simulation time, namely, the jet propagates straight.
Magnetic fields also play an important role in the suppression of the Rayleigh–Taylor instability. Figure 4 shows the v_{z} distribution on the xyplane at z = 70 kpc for model A, B, and C. We observe that kink instability for models A and B bend the jet away from its initial launch axis (x = y = 0). Furthermore, for models A and B, the jets clearly separate between the beam flow (yellow region) and cocoon gas (red region), namely, the low mixing ratio between the beam and cocoon gas. Meanwhile, the lowmagnetized jet of model C is not in the development of kink instability, but in that of the Rayleigh–Taylor and Kelvin–Helmholtz instabilities. Similar results have been presented for relativistic MHD jet simulations in Mignone et al. (2010) and Mukherjee et al. (2020). In addition to this, our results are supported by the linear stability analysis for a relativistic nonrotating jet, which indicates a strong magnetic field can suppress a growth rate of the Kelvin–Helmholtz instability (Bodo et al. 2013). The onset condition of the Rayleigh–Taylor instability is given analytically by ρ_{jet} > ρ_{cocoon} in the hydrodynamic case. Notably, in their analysis, Komissarov et al. (2019) found that jets are stable for the Rayleigh–Taylor instability mode when ℳ_{A} < 40. Thus, the jet is in an unstable mode of the Rayleigh–Taylor instability for model C (see Fig. 2 and Table 1). The right panel of Fig. 4 shows the Rayleigh–Taylor and the Kelvin–Helmholtz mode forming a large crosssection, ∼10 kpc, of positive velocity field (v_{z} > 0: from blue to yellow region), and fingerlike structures. Owing to the highmixing ratio between the beam and cocoon gas, the jet of model C is decelerated.
Fig. 4. Slices (in the x − y plane at z = 80 kpc) of the distribution of the zdirection component of velocity for models A, B, and C. 
Figure 5 shows the slices of the xdirection component of the magnetic field for models A, B, and C, respectively. While the jet beam is injected with a pure toroidal magnetic field, the inverse field is randomly distributed in the cocoon. The toroidal magnetic field in the cocoon is weaker, that is, closer than z < 40 kpc. Further, reversing fields are dissipated in the cocoons of models A and B at z < 40 kpc. We describe the field structures in the beams in more detail in Sect. 5.1.
Fig. 5. Slices (in the y − z plane) of the xdirection component distribution of the magnetic field for models A, B, and C, respectively. 
A magnetic filament develops in the cocoon at z > 40 kpc in model A because the strong initial toroidal magnetic field that flows down with the countercurrent creates sufficient magnetic tension to suppress turbulent motion. The typical length of the filament is several kpc, longer than the filament in model B. Due to shock compression, filaments are formed around the jet head, and have stronger magnetic fields than the injected ones. In contrast, for the cocoon of model C, the smallscale turbulence is excited. To discuss the lengthscale of the magnetic filament quantitatively, we evaluate the lengthscales parallel to the magnetic field as (Schekochihin et al. 2004; Bodo et al. 2011; Mukherjee et al. 2020):
Figure 6 shows the volumeweighted probability distribution function (PDF) of the L_{∥} in the jet cocoon where the electron temperature is higher than 10^{8} K and z > 40 kpc. It can be seen from the PDF that the typical value of the length scale increases with stronger magnetization. The cocoon of model C is filled in smaller vortices, whose typical length scale is about 0.3 kpc. The volume occupations of the length scales that are longer than 1 kpc are 6.5, 16.5, and 28.0% for models A, B, and C, respectively. These trends are consistent with the results of previous MHD studies (Mukherjee et al. 2020).
Fig. 6. Probability distribution functions of the characteristic length scale parallel to the magnetic field in the cocoon where z > 40 kpc for model A (red dotted), B (blue solid), and C (green dashed) at the end of the simulations. We define the cocoon as grids with the electron temperature higher than 10^{8} K and z > 40 kpc. 
4.2. Temperature distribution
Figure 7 shows the distribution of the electron temperature at the yzplane for models A, B, and C. At first glance, the electron temperature in the jet is proportional to the strength of injected magnetic field. This implies that the subgrid for turbulence heating plays an important role in the evolution of the electron temperature (see also Sect. 4.3). Here, we recall that the electron heating fraction for turbulence is proportional to inverse plasma beta .
Fig. 7. Slices (in the y − z plane) of the electron temperature distribution for models A, B, and C, respectively. 
Thermal electrons propagating through the beam are not subject to the heating energy of the internal shock, but the electrons are heated in the jet termination region. Subsequently, hot electrons are stored in the cocoon. Although this physical image is similar to the result of the twodimensional case (see Fig. 1 in Ohmura et al. 2020), the difference between them is that the turbulence heating also works in the beam. Thus, the electrons are heated locally in the beam for models A and B.
In contrast to electrons, protons receive most of the shock heating, as we set f_{e} = 0.05. Thus, proton temperatures are several tens of times higher than electron temperatures in the cocoons for all models (Fig. 8). Electrons are in the relativistic temperature in range of T_{e} ∼ 10^{9} − 10^{10} K. Hotter electrons are located along magnetized filamentary structures formed by shock compression for models A and B. In particular, electron temperatures are higher than the protons’ temperature in some filaments of model A. Meanwhile, for model C, the distribution of the ratio of the proton to electron temperature is monochromatic, kT_{p}/kT_{e} ∼ 40, in the cocoon.
Fig. 8. Slices (in the y − z plane) of the ratio of proton to electron temperature distributions for models A, B, and C, respectively. 
4.3. Lobe energetics
The left panel of Fig. 9 displays the time evolution of different energy components of the cocoon for all models. The proton thermal energy is the dominant energy component of the cocoon (U_{p} ≫ U_{e}), that is, cocoons are supported by the proton pressure (see red and blue lines in the left panel of Fig. 9). We confirm that the kinetic energy is comparable to the proton thermal energy and that 20 percent of total injected energy is converted to the ICM at t = 10 Myr. Because the electron heating fraction of turbulence is an increasing function of , there is a positive correlation between the electron thermal energy and field strength (see dashed lines, solid lines, and dotted lines for magnetic energy and electron thermal energy in the left panel of Fig. 9).
Fig. 9. Energies in the cocoon as a function of time for all models. Left: Time evolution of different energy components of the cocoon for model A (dotted lines), B (solid lines), and C (dashed lines), respectively. We define the cocoon as grids with an electron temperature higher than 10^{8} K. Right: Time evolution of the ratio between the magnetic field and the electron energy for model A (dotted lines), B (solid lines), and C (dashed lines). 
In the right panel of Fig. 9, we plot the time evolution of the ratio between the magnetic and electron thermal energies in the cocoon. The energy ratio for models A and B saturates at ∼0.7. The electron pressure is described by p_{e} = (γ_{e} − 1)u_{e} = 2u_{e}/3 for γ_{e} = 5/3 (from Eq. (14)); hence, it is in pressure equilibrium with the magnetic field for model A and B. Meanwhile, in the case of model C, the electron pressure is larger than the magnetic pressure.
Thermal electrons evolve while their thermal energy is added by a large amount of dissipated energy due to shocks and turbulence. If the gas reaches a turbulence equilibrium discussed in Sadowski et al. (2017), the final temperature ratio is determined by the electron heating model, which is given by
where we adopt γ_{p} = γ_{e} = 5/3. For models A and B, a quasisteady state of MHD turbulence is observed in the T_{e}/T_{p} − β_{p} histogram (left panel of Figs. 10a,b). This histogram plots the gas stored in the cocoons. We observe that the gas distribution follows along the dashed line, which is plotted as Eq. (24). This indicates that energy components (U_{p}, U_{e}, and U_{mag}) evolve following the electron heating model of MHD turbulence. Therefore, the electron heating model of turbulence plays a significant role in the determination of the gas thermal evolution in our models. The heating ratio of protons to electrons, Q_{p}/Q_{e} in the turbulence heating model saturates at ∼30 for β_{i} > 10, and therefore the minimum temperature ratio is located at (T_{e}/T_{p}) = 1/30 (see Eq. (17)). Another view on the turbulence equilibrium is the relationship between U_{e} and U_{mag}, shown in the right panel of Fig. 10. The gas distributes along the line, U_{e} = U_{mag}, in this histogram. We confirmed that the volume occupations of the gas where 0.666 < U_{e}/U_{mag} < 3 are 0.82 and 0.80 for model A and B. Therefore, the electrons evolve toward energy equipartition (pressure equilibrium) with the magnetic energy when β_{p} < 10 in the cocoon for models A and B. Note that The fraction of electron heating become a constant, f_{e} ∼ 1/30, when the proton plasmaβ_{p} is higher than 10, namely, the fraction of electron heating does not depend on the magnetic fields in this regime. Thus, a part of gas does not distribute along the line that U_{e} = U_{mag} in a U_{e} − U_{mag} histogram.
Fig. 10. Relationship of three energy components: electron thermal energy, proton thermal energy, and magnetic field energy. Left: T_{e}/T_{p} − β_{p} histogram for regions in the cocoon for model A (a), model B (b), and model (c) at t = 9.52 − 9.94 and 13.02 Myr, respectively. The dashed line depicts the electron to proton temperature ratio corresponding to the equilibrium state for plasma β_{p}, as implied by the turbulence heating in Eq. (24). Right: Same as left panel, but displaying U_{e} − U_{mag} histogram. Dashed line plots U_{e} = U_{mag}. 
For model C, proton plasmaβ_{p} is higher than 10 across the whole region of the cocoon. Thus, U_{e}/U_{p} ∼ 1/30 in the cocoon for model C after 5 Myr and there is no relation between the electron thermal energy and the magnetic field energy (see Fig. 10c). The volume occupations of the gas where 0.666 < U_{e}/U_{mag} < 3 for model C are 0.35. Furthermore, the magnetic energy is subdominant compared with the electron thermal energy. Notably, U_{e}/U_{mag} saturates at ∼0.4 for model C in the right panel of Fig. 10. However, the saturation value depends on the proton plasma beta when β_{p} ≫ 10 in a cocoon.
5. Discussion
5.1. Smallscale dissipation in jet beam
As we report in Sect. 3, beams suffer MHD instabilities. The growth of instabilities leads to the formation of current sheets, where magnetic reconnection takes place. Magnetic reconnection is a dissipation mechanism that can energize nonthermal particles. Notably, although we do not explicitly deal with resistivity, magnetic reconnection arises due to numerical dissipation. The magnetic energy dissipation rate is given by ηj^{2}, where η is the resistivity, and j is the current density. However, it is difficult to measure the numerical resistivity in ideal MHD simulations. Thus, following Zhang et al. (2017), we quantify dissipation to be proportional to the strength of j ⋅ E, where E = −v × B is the electric field.
Figure D.1 displays the volumerenders of a physical quantity j ⋅ E at times when jets reach 60 kpc. Notably, the color bar scale of the panel (c) is 0.2 times narrower than that of panels (a) and (b). Peaks of dissipation take place around the jet head for all models due to the shock compression at the termination shocks. This feature can be regarded as typical for powerful FRII type jets. When model A enters the nonlinear phase for the kink mode, the beam core is fragmented, which is shown in the beam at 30 < z < 50 kpc in the left panel of Fig. 5. This fragmented structure identifies as a dissipation region, which shows a smallscale periodic structure at 30 < z < 50 kpc for model A. Meanwhile, the value of j ⋅ E is high at the shear layer between the beam and cocoon for model B. The growth of the Kelvin–Helmholtz mode drive gas mixing between the beam and cocoon in the shear layer (Mignone et al. 2010; Mukherjee et al. 2020). Therefore, the beam radius of model C in Fig. D.1 looks larger than that of models A and B.
After the jets propagate at 95 kpc, the dissipative structures of the three models are significantly different (Fig. D.1). We observe the fragmentation structures at 30 < z < 60 kpc for model A at this time. The magnetic energy is dissipated in this region, hence, the flow is laminar downstream of it (60 < z < 70 kpc). Dissipative spots are formed by shock, in particular, and the shock is induced by magnetic pinching at z = 70 kpc. For model B, we observe the abrupt change in the flow direction by the development of a beam kink at z = 60, 70, 80 kpc. Hence, these localized spots, where the flow is shocked and bent, could be reconnection layers, where efficient particle acceleration would take place. Such beam disruption can explain the formation mechanism of double hotspots in the western lobe of Cygnus A (Carilli & Barthel 1996) and of multiple knots observed 3C273 (Uchiyama et al. 2006). Meanwhile, the model C jet does not yield the formation of multiple hotspots. Although the model C jet has the dissipative spot at the jet head, the dissipation ratio of the magnetic energy is uniformly distributed. Therefore, the jet at the later phase of model C has a feature of FRI type jets such as M87, which has an extended diffusive radio lobe without a hotspot (Laing et al. 2011). Jet deceleration by the development of Rayleigh–Taylor instabilities is a possible explanation of FR dichotomy, consistently with previous simulations (Massaglia et al. 2016; Rossi et al. 2020).
5.2. Comparison with observations
We discuss that our twotemperature MHD models connect to the observational results. In particular, we focus on the relationship between the jet mechanical power and radio power. Observational data sets of Xray and radio properties are adopted from the PhD thesis by Laura Bîrzan (Rafferty 2007) and Rafferty et al. (2006).
5.2.1. Radio power
The radio power is obtained by the sum of the radio emissivity for the synchrotron emission in the optical thin limit, as follows:
where j_{ν} and ν are the synchrotron emissivity and the observed frequency, respectively. We used the formula of the synchrtoron emissitity provided in Pacholczyk (1970).
We assumed a single powerlaw distribution for nonthermal electrons, dN/dγ = N_{0}γ^{−s}, where s, N_{0}, and γ are the electron energy index, number of nonthermal electrons, and Lorentz factor of nonthermal electrons. Meanwhile, the twotemperature MHD simulations account only for the evolution of thermal electrons. Because previous studies suggest that the nonthermal electron energy is proportional to the thermal gas energy or magnetic energy as a firstorder approximation, this approximation is adopted to estimate the radio power from numerical simulations (e.g., Gomez et al. 1995). We follow this approximation in the current study. We adopt the twotype model for approximation of nonthermal electrons as cases 1T and 2T. Case 1T is N_{0} = C_{0}u_{p}, and case 2T is N_{0} = C_{0}u_{e}. Here, , and η = 0.2 is a parameter, which is the ratio of the nonthermal electron energy density to the electron thermal energy density. Furthermore, we set γ_{min} = 100 and γ_{max} → ∞. Case 1T is corresponds to the singletemperature model, namely, T_{e} = T_{p}. This model is motivated by the prior works demonstrating the radio emission from AGN jets on assuming thermal equilibrium between protons and electrons (e.g., Gomez et al. 1995). We assume that an electron energy index, s, is 2.05. For all calculations, the viewing angle, which is in respect to the zaxis, is 80°, and the observed frequency is 144 MHz.
We list the calculated radio powers in Table 3. The radio power of model A1T, which is most prominent model, is two orders of magnitude higher than that of model C2T. The 1T models have same order of nonthermal electron energy, because the proton temperatures have roughly the same value in lobes (see left panel of Fig. 9). Notably, 1T models assume that nonthermal electron energy is proportional to the proton temperature. Meanwhile, electron temperatures vary for the three simulation models, they are proportional to the inverse the proton plasmaβ. Thus, there are a large scatter in the radio powers of 2T models, greater than that seen for 1T models. This indicates that the radio power of twotemperature models is sensitive to the magnetic field energy.
Radio power and the amount of PdV work for the simulation results.
5.2.2. Xray cavity and jet mechanical power
Because our simulation assumes a constant energy input corresponding to the active phase of the jet, we can calculate the true mechanical power in units of erg s^{−1}. At the same time, we obtain snapshot quantities and the gas pressure of surrounding ICM through the actual observations. Therefore, the same observation method is adopted in this work. Meanwhile, mechanical power is measured using PdV work, as follows in Xray observations.
where t_{age} is the outburst age of jets. Lobes and ICM are approximately in the pressure equilibrium state. Thus, we used the initial ICM pressure around middle of lobe p_{gas} ∼ 2.0 × 10^{−10} erg cm^{−3} to calculate the cavity power (see Fig. 1). Three estimations are commonly used for the outburst age: the buoyancy time, t_{buoy}, refill time, t_{r}, and sound crossing time, t_{c}, generally t_{c} < t_{buoy} < t_{r} (McNamara & Nulsen 2007). It is appropriate to use the sound crossing time in our model, because jets are in an active phase. The cavity volume, V, is calculated by integrating numerical grids with an electron temperature exceeding 10^{8} K. We confirm that the region whose electron temperature exceeds 10^{8} K has a lower density than ICM.
First, we mention the morphological properties of the Xray cavity. In Fig. 11, we compare our simulation of model B with observations about the relationship between the projected distance from the core to the cavity center, R, and the projected semiminor axis of the cavity, b. Because R and b for model A and model C have similar values as for model B, we only show the result of model B in Fig. 11. The axis of the cavity of model B is roughly consistent with the observation, though the cavity is narrow compared with observed ones. This long and narrow structure is characteristic of a powerful (kineticdominated) jet, observed in a number of previous studies (e.g., Massaglia et al. 2016; Perucho et al. 2019; Mukherjee et al. 2020). To create a broad cavity, the jet must propagate slowly to have sufficient time to expand. Thus, a simple solution to forming a broad cavity is to model a lowdensity jet (Ohmura et al. 2021) and/or a lowpower jet (Mukherjee et al. 2020). Otherwise, a longperiodic precession may play an important role in forming observed broadened cavities (Horton et al. 2020). The jet for model B decelerates and has precession due to the development of largescale kink modes. However, the speed of the lateral expansion is also decelerated to about the sound speed of ICM. Therefore, we cannot expect the axis ratio (R/b) to decreases at later times (Mukherjee et al. 2020).
Fig. 11. Projected distance from core to cavity center, R, versus Projected semiminor axis of the cavity, b. Black circles show the radiofilled cavities taken from Rafferty et al. (2006). Blue square shows our result for model B at t = 9.94 Myr. Because R and b values for models A and C have similar values to model B, we do not show them. 
Next, we discuss the energy and age by the observational method. The cavity energies are calculated as E_{cav} = 4p_{gas}V = 2.4 × 10^{59}, 2.6 × 10^{59}, and 3.7 × 10^{59} erg for models A, B, and C, respectively. Further, the sound crossing time t_{c} = R/c_{s, ICM} is 51.6 Myr, where we adopt the following assumption: R = 45 kpc and c_{s} = 830 km s^{−1} (corresponding to ICM temperature T = 5 KeV). Therefore, we estimate the mechanical power P_{cav} = 1.5 × 10^{44}, 1.6 × 10^{44}, and 2.3 × 10^{44} erg s^{−1} for model A, B, and C, respectively. Hence, the actual age from our simulations is ∼10 Myr for all models, and the sound crossing time underestimates the mechanical power by a factor of ∼5. Even if we estimate the mechanical power using the simulation time, the injection energy of the jet is about ten times higher. The reason is attributed to the conversion of the jet energy to the thermal energy of ICM through shocks and sound waves. Furthermore, the cocoon of model B is still overpressured with respect to the ICM while the mechanical power is measured under the assumption that the radio lobe and ICM have reached in the pressure equilibrium state.
5.2.3. Relationship between radio power and mechanical power
The plot of the jet mechanical power P_{cav} versus the synchrotron radio power P_{radio} is shown in Fig. 12. It provides physical insights for jet energetics, including nonradiating proton thermal energy. Naively, protons can be energetically dominant over radiatively inefficient lobes (P_{cav} ≫ P_{radio}), as the pressure of nonradio emitting protons supports the expansion of the cocoon. Our jets are active during the simulation time. Therefore, we compare our results with radiofilled cavities, except for radioghost cavities. Although we plot the radiofilled cavities including the intermediate cases in Fig. 12, the samples have a large scatter in the P_{cav}/P_{radio} ∼ 1 − 1000 relation.
Fig. 12. Radio synchrotron power P_{radio} versus jet mechanical power estimated from the Xray cavity system: (adopted from Laura Bîrzan’s PhD thesis Rafferty 2007). Filled black symbols show radiofilled cavities (which include intermediate cases). Symbols and wide error bars denote the values of the mechanical power calculated using the sound speed. Open squares and filled squares show our results of cases 1T and 2T for models A (red), B (blue), and C (green), respectively. The diagonal dotted lines (dashed lines) represent ratios of constant mechanical power to radio luminosity. Blue dashed lines depict the total injection energy of our jet model, given by Eq. (20). 
We find that our twotemperatures model is able to explain radiatively inefficient lobes. In case 2T of all models, the lobes tend to be more radiatively inefficient than those of case 1T. Because the electrons lack thermal energy compared with case 1T, the radio powers are weak. Meanwhile, protons have a large contribution for the cavity power, P_{cav}. Thus, radiative efficiencies, P_{cav}/P_{radio}, for case 2T are 10–30 higher than the results for the 1T case. The ratio of the radio power between models A and C for case 2T is higher than that for case 1T. This difference is attributed to electron heating being coupled with the strength of magnetic fields (see Sect. 4.3).
This result indicates that the pure protonselectrons jet has difficulty to form radiatively efficient lobes, such as Cygnus A, which is located at P_{cav}/P_{radio} = 1 in Fig. 12, without exotic processes. One of the exotic processes, herein, is the efficient acceleration for electrons. To create a radiationefficient lobe, an acceleration mechanism of η ≫ 10 is needed, where the energy of nonthermal electrons is at least an order of magnitude larger than that of the thermal ones. Alternatively, thermal electrons are dominant over protons, as they are efficiently heated by shocks and turbulence. Nevertheless, we have not observed this image from several PIC simulations of shock and turbulence (e.g., Crumley et al. 2019; Zhdankin et al. 2019, 2021). Another possibility to form radiatively efficient lobes is a strong magnetic field of jets, β_{p} ≪ 1. However, this possibility is not favored in some observations (Croston et al. 2005; Isobe & Koyama 2015).
The existence of a large number of positrons could likewise explain the radiatively efficient lobes. Our simulations model the pure electronproton jet. Thus, to achieve P_{cav}/P_{radio} = 1 in our simulations of models A and B (case 2T), the number density of leptons must be at least a hundred times larger than that of protons, because P_{radio} is proportional to their number density. Notably, the plasma momentum would be represented by protons under this assumption, because the protons mass is a thousand magnitude greater than the lepton mass. Thus, it is expected that there is no difference in the electron heating process. Meanwhile, in the case of model C, significant population of the pairplasma is needed, because P_{cav}/P_{rad} ∼ 1000. In this condition, the electron (and positron) heating process would change, and, furthermore, the MHD approximation would not be guaranteed. Finally, analytical models of electronpositronproton mixture jets likewise achieved consistent results with observed FRII radio lobes in studies by Ito et al. (2008), Kawakatu et al. (2016) and Kino et al. (2012). However, these models do not consider the electron heating process at turbulence and shocks. Therefore, construction of the new model for the mixture jet based on twotemperature simulations is necessary.
We discuss only the radiofill cavities, as our jet is always active during the simulation time. However, several radio and Xray observations imply multiepisodic jet activity (e.g., Wise et al. 2007; Maccagni et al. 2020). These images make it difficult to model the radio lobes, leading to a large scatters in the P_{cav} − P_{radio} relation. To understand the physical condition of the radio lobes, it is important that we find radio lobes in active. To this end, the future radio survey by the Square Kilometer Array is of great value for studying the radio lobes.
5.3. Observational implications
We showed that electrons remain in a transrelativistic temperature, namely, the electron energy ranges over a few MeV. Although it is quite difficult to obtain the observational signals from these electrons, there are several possibilities. One possibility is that the highquality Sunyaev–Zel’dovich (SZ) radio observations are able to obtain information on the pressure of thermal electrons (Pfrommer et al. 2005). However, our results indicate that the thermal SZ signal from the radio lobe would not be able to detect this, because the electron pressure is lower than the pressures of ICM protons and radio lobe electrons in our model. Another possibility is the Cosmic Micro Background inverseCompton (CMBIC) spectrum. If the energy of thermal electrons exceeds GeV, the thermal CMBIC spectrum would be observed in Xray observations. However, these components are not yet detected in the contest of radio lobes. In contrast, thermal electrons whose energies range between a few MeV scatter photons in the infrared and optical range (Enßlin & Sunyaev 2002).
6. Summary
We carried out twotemperature MHD simulations for three models whose jets have different magnetic fields. The jets propagated along 90 kpc in the cluster center, whose environment is roughly consistent with the Cygnus cluster. We studied the dynamics and electron heating in the subgrid scale of the radio lobes and, thus, two subgrid electron heating mechanisms were considered. Because the jets have both the turbulence and the strong shock waves, we used the subgrid models at shock waves and at the turbulence in a hybrid manner. The main findings achieved in this study are as follows:

(i)
Strongly magnetized jets suffer from the development of a nonaxisymmetric, currentdriven kink mode. Meanwhile, weakly magnetized jets are decelerated by the highmixing ratio between the jet beam and cocoon gas, induced by a Rayleigh–Taylor and Kelvin–Helmholtz instability. We show that the Alfvén crossing time, τ_{kink}, is a good indicator for the time scale of the development of the kink mode, even if the jets have purely toroidal fields at the injected region.

(ii)
Electrons heat up at the jet termination region and hot electrons are stored in the cocoon. The electron heating fraction for turbulence is proportional to the inverse plasma beta, . Therefore, a jet with a strong magnetic field has a higher electron temperature in the cocoon.

(iii)
Smallscale turbulence develops in the weakly magnetized jets. In contrast, the strongly magnetized jets have magnetized filaments, as the magnetic tension suppresses the turbulence motion.

(iv)
The protons are energetically dominant over the electrons in the cocoon. First, most of the bulk kinetic energy of the jet is converted into thermal energy of protons through shocks. Second, while magnetic fields are relatively strong, shockedelectrons stored in the cocoon evolve toward energy equipartition with magnetic energy through turbulent dissipation.

(v)
The strong current is induced by the kink instability. Therefore, hightemperature and highmagnetization multiple hotspots are formed in the beam.

(vi)
The lowdensity cavity of our model is narrow compared with observed radiofilled cavities. This suggests that the density of jets is slightly lower than that determined by our model. The propagation velocity is faster than the sound speeds of ICM and some fractions of the jet energy are converted into the thermal energy ICM. Thus, the jet mechanical energy estimated by Xray observation is ten times lower than the injected kinetic energy in our simulations.

(vii)
Radio powers estimated from the electron thermal energy are ten times lower than those estimated from the proton thermal energy, which correspond to the onetemperature approximation. Twotemperature models quantitatively explain the radiatively inefficient lobes (P_{cav}/P_{radio} ∼ 100). Furthermore, our results indicate that it is difficult to explain radiatively efficient lobes, such as Cygnus A, unless nonthermal energy is more than one order of magnitude larger than the thermal electron energy.
In the current study, we only focused on the results of the property of the jetted plasma. However, the weak shocks around radio lobes are frequently observed in Xray observations (e.g., Snios et al. 2018). We confirm that shockedICM plasma occurs in twotemperature states (see Fig. 8), hence, we aim to report thermodynamics and Xray properties of shockedICM in Paper II (Ohmura et al. 2023).
Acknowledgments
We thank the anonymous referee for the useful comments that greatly improved the presentation of the paper. This work was supported by JSPS KAKENHI Grant Numbers JP22K14032 (T.O.), 19K03916, and 22H01272 (M.M.). Our numerical computations were carried out on the Cray XC50 at the Center for Computational Astrophysics of the National Astronomical Observatory of Japan. The computation was carried out using the computer resource by Research Institute for Information Technology, Kyushu University. This work was also supported in part by MEXT as a priority issue (Elucidation of the fundamental laws and evolution of the universe) to be tackled by using postK Computer and JICFuS and by MEXT as “Program for Promoting Researches on the Supercomputer Fugaku” (Toward a unified view of the universe: from large scale structures to planets). This research used computational resources of the K computer/OakforestPACS provided by JCAHPC through the HPCI System Research project (project id: hp190125, hp200092).
References
 Aloy, M. A., Ibáñez, J. M., Martí, J. M., Gómez, J. L., & Müller, E. 1999, ApJ, 523, L125 [CrossRef] [Google Scholar]
 Beck, R., & Krause, M. 2005, Astron. Nachr., 326, 414 [Google Scholar]
 Bîrzan, L., Rafferty, D. A., McNamara, B. R., Wise, M. W., & Nulsen, P. E. J. 2004, ApJ, 607, 800 [Google Scholar]
 Bîrzan, L., McNamara, B. R., Nulsen, P. E. J., Carilli, C. L., & Wise, M. W. 2008, ApJ, 686, 859 [Google Scholar]
 Bodo, G., Massaglia, S., Ferrari, A., & Trussoni, E. 1994, A&A, 283, 655 [NASA ADS] [Google Scholar]
 Bodo, G., Cattaneo, F., Ferrari, A., Mignone, A., & Rossi, P. 2011, ApJ, 739, 82 [NASA ADS] [CrossRef] [Google Scholar]
 Bodo, G., Mamatsashvili, G., Rossi, P., & Mignone, A. 2013, MNRAS, 434, 3030 [Google Scholar]
 Braginskii, S. I. 1965, Rev. Plasma Phys., 1, 205 [Google Scholar]
 Carilli, C. L., & Barthel, P. D. 1996, A&ARv, 7, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Croston, J. H., Hardcastle, M. J., Harris, D. E., et al. 2005, ApJ, 626, 733 [NASA ADS] [CrossRef] [Google Scholar]
 Crumley, P., Caprioli, D., Markoff, S., & Spitkovsky, A. 2019, MNRAS, 485, 5105 [NASA ADS] [CrossRef] [Google Scholar]
 Dermer, C. D., Liang, E. P., & Canfield, E. 1991, ApJ, 369, 410 [Google Scholar]
 Enßlin, T. A., & Sunyaev, R. A. 2002, A&A, 383, 423 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Fabian, A. C. 1994, ARA&A, 32, 277 [Google Scholar]
 Fabian, A. C. 2012, ARA&A, 50, 455 [Google Scholar]
 Fabian, A. C., Sanders, J. S., Ettori, S., et al. 2000, MNRAS, 318, L65 [Google Scholar]
 Fanaroff, B. L., & Riley, J. M. 1974, MNRAS, 167, 31P [Google Scholar]
 Gomez, J. L., Marti, J. M. A., Marscher, A. P., Ibanez, J. M. A., & Marcaide, J. M. 1995, ApJ, 449, L19 [Google Scholar]
 Guo, X., Sironi, L., & Narayan, R. 2018, ApJ, 858, 95 [NASA ADS] [CrossRef] [Google Scholar]
 Hardcastle, M. J., & Croston, J. H. 2020, New Astron. Rev., 88, 101539 [Google Scholar]
 Horton, M. A., Krause, M. G. H., & Hardcastle, M. J. 2020, MNRAS, 499, 5765 [NASA ADS] [CrossRef] [Google Scholar]
 Howes, G. G. 2010, MNRAS, 409, L104 [NASA ADS] [Google Scholar]
 Isobe, N., & Koyama, S. 2015, PASJ, 67, 77 [NASA ADS] [CrossRef] [Google Scholar]
 Ito, H., Kino, M., Kawakatu, N., Isobe, N., & Yamada, S. 2008, ApJ, 685, 828 [CrossRef] [Google Scholar]
 Jones, T. W., Ryu, D., & Engel, A. 1999, ApJ, 512, 105 [NASA ADS] [CrossRef] [Google Scholar]
 Kawakatu, N., Kino, M., & Takahara, F. 2016, MNRAS, 457, 1124 [NASA ADS] [CrossRef] [Google Scholar]
 Kawazura, Y., Barnes, M., & Schekochihin, A. A. 2019, Proc. Natl. Acad. Sci., 116, 771 [NASA ADS] [CrossRef] [Google Scholar]
 Kawazura, Y., Schekochihin, A. A., Barnes, M., et al. 2020, Phys. Rev. X, 10, 041050 [NASA ADS] [Google Scholar]
 King, I. 1962, AJ, 67, 471 [Google Scholar]
 Kino, M., Kawakatu, N., & Takahara, F. 2012, ApJ, 751, 101 [CrossRef] [Google Scholar]
 Komissarov, S. S., Gourgouliatos, K. N., & Matsumoto, J. 2019, MNRAS, 488, 4061 [NASA ADS] [CrossRef] [Google Scholar]
 Laing, R. A., Guidetti, D., Bridle, A. H., Parma, P., & Bondi, M. 2011, MNRAS, 417, 2789 [NASA ADS] [CrossRef] [Google Scholar]
 Maccagni, F. M., Murgia, M., Serra, P., et al. 2020, A&A, 634, A9 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Massaglia, S., Bodo, G., Rossi, P., Capetti, S., & Mignone, A. 2016, A&A, 596, A12 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Massaglia, S., Bodo, G., Rossi, P., Capetti, S., & Mignone, A. 2019, A&A, 621, A132 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mathews, W. G., & Guo, F. 2010, ApJ, 725, 1440 [NASA ADS] [CrossRef] [Google Scholar]
 Matsukiyo, S. 2010, Phys. Plasmas, 17, 042901 [NASA ADS] [CrossRef] [Google Scholar]
 Matsumoto, J., & Masada, Y. 2013, ApJ, 772, L1 [NASA ADS] [CrossRef] [Google Scholar]
 Matsumoto, Y., Asahina, Y., Kudoh, Y., et al. 2019, PASJ, 71, 83 [NASA ADS] [CrossRef] [Google Scholar]
 McNamara, B. R., & Nulsen, P. E. J. 2007, ARA&A, 45, 117 [NASA ADS] [CrossRef] [Google Scholar]
 McNamara, B. R., Wise, M. W., Nulsen, P. E. J., et al. 2001, ApJ, 562, L149 [NASA ADS] [CrossRef] [Google Scholar]
 Mendygral, P. J., Jones, T. W., & Dolag, K. 2012, ApJ, 750, 166 [NASA ADS] [CrossRef] [Google Scholar]
 Mignone, A., Rossi, P., Bodo, G., Ferrari, A., & Massaglia, S. 2010, MNRAS, 402, 7 [NASA ADS] [CrossRef] [Google Scholar]
 Mignone, A., Striani, E., Tavani, M., & Ferrari, A. 2013, MNRAS, 436, 1102 [NASA ADS] [CrossRef] [Google Scholar]
 Mizuno, Y., Lyubarsky, Y., Nishikawa, K.I., & Hardee, P. E. 2009, ApJ, 700, 684 [NASA ADS] [CrossRef] [Google Scholar]
 Moll, R., Spruit, H. C., & Obergaulinger, M. 2008, A&A, 492, 621 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mukherjee, D., Bodo, G., Mignone, A., Rossi, P., & Vaidya, B. 2020, MNRAS, 499, 681 [Google Scholar]
 Mukherjee, D., Bodo, G., Rossi, P., Mignone, A., & Vaidya, B. 2021, MNRAS, 505, 2267 [NASA ADS] [CrossRef] [Google Scholar]
 Nakahara, S., Doi, A., Murata, Y., et al. 2019, ApJ, 878, 61 [NASA ADS] [CrossRef] [Google Scholar]
 Ohmura, T., Machida, M., Nakamura, K., et al. 2019, Galaxies, 7, 14 [NASA ADS] [CrossRef] [Google Scholar]
 Ohmura, T., Machida, M., Nakamura, K., Kudoh, Y., & Matsumoto, R. 2020, MNRAS, 493, 5761 [CrossRef] [Google Scholar]
 Ohmura, T., Ono, K., Sakemi, H., et al. 2021, ApJ, 910, 149 [NASA ADS] [CrossRef] [Google Scholar]
 Ohmura, T., Machida, M., & Akamatsu, H. 2023, A&A, 679, A161 (Paper II) [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Pacholczyk, A. G. 1970, Radio Astrophysics. Nonthermal Processes in Galactic and Extragalactic Sources (San Francisco: Freeman) [Google Scholar]
 Perley, R. A., Dreher, J. W., & Cowan, J. J. 1984, ApJ, 285, L35 [NASA ADS] [CrossRef] [Google Scholar]
 Perucho, M., Martí, J.M., Quilis, V., & Ricciardelli, E. 2014, MNRAS, 445, 1462 [NASA ADS] [CrossRef] [Google Scholar]
 Perucho, M., Martí, J.M., & Quilis, V. 2019, MNRAS, 482, 3718 [NASA ADS] [CrossRef] [Google Scholar]
 Pfrommer, C., Enßlin, T. A., & Sarazin, C. L. 2005, A&A, 430, 799 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Pfrommer, C., Pakmor, R., Schaal, K., Simpson, C. M., & Springel, V. 2017, MNRAS, 465, 4500 [NASA ADS] [CrossRef] [Google Scholar]
 Porth, O., & Komissarov, S. S. 2015, MNRAS, 452, 1089 [Google Scholar]
 Rafferty, F. 2007, Ph.D. thesis, Ohio University, USA [Google Scholar]
 Rafferty, D. A., McNamara, B. R., Nulsen, P. E. J., & Wise, M. W. 2006, ApJ, 652, 216 [NASA ADS] [CrossRef] [Google Scholar]
 Ressler, S. M., Tchekhovskoy, A., Quataert, E., Chandra, M., & Gammie, C. F. 2015, MNRAS, 454, 1848 [Google Scholar]
 Rossi, P., Bodo, G., Massaglia, S., & Capetti, A. 2020, A&A, 642, A69 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ryu, D., Kang, H., Hallman, E., & Jones, T. W. 2003, ApJ, 593, 599 [NASA ADS] [CrossRef] [Google Scholar]
 Sadowski, A., Wielgus, M., Narayan, R., et al. 2017, MNRAS, 466, 705 [NASA ADS] [CrossRef] [Google Scholar]
 Schaal, K., & Springel, V. 2015, MNRAS, 446, 3992 [NASA ADS] [CrossRef] [Google Scholar]
 Scheck, L., Aloy, M. A., Martí, J. M., Gómez, J. L., & Müller, E. 2002, MNRAS, 331, 615 [NASA ADS] [CrossRef] [Google Scholar]
 Schekochihin, A. A., Cowley, S. C., Taylor, S. F., Maron, J. L., & McWilliams, J. C. 2004, ApJ, 612, 276 [NASA ADS] [CrossRef] [Google Scholar]
 Scheuer, P. A. G. 1974, MNRAS, 166, 513 [NASA ADS] [CrossRef] [Google Scholar]
 Smith, D. A., Wilson, A. S., Arnaud, K. A., Terashima, Y., & Young, A. J. 2002, ApJ, 565, 195 [CrossRef] [Google Scholar]
 Snios, B., Nulsen, P. E. J., Wise, M. W., et al. 2018, ApJ, 855, 71 [NASA ADS] [CrossRef] [Google Scholar]
 Spitzer, L. 1962, Physics of Fully Ionized Gases (New York: Interscience) [Google Scholar]
 Stepney, S., & Guilbert, P. W. 1983, MNRAS, 204, 1269 [NASA ADS] [Google Scholar]
 Svensson, R. 1982, ApJ, 258, 335 [Google Scholar]
 Tchekhovskoy, A., & Bromberg, O. 2016, MNRAS, 461, L46 [NASA ADS] [CrossRef] [Google Scholar]
 Tran, A., & Sironi, L. 2020, ApJ, 900, L36 [NASA ADS] [CrossRef] [Google Scholar]
 Uchiyama, Y., Urry, C. M., Cheung, C. C., et al. 2006, ApJ, 648, 910 [NASA ADS] [CrossRef] [Google Scholar]
 Vink, J., Broersen, S., Bykov, A., & Gabici, S. 2015, A&A, 579, A13 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Weinberger, R., Ehlert, K., Pfrommer, C., Pakmor, R., & Springel, V. 2017, MNRAS, 470, 4530 [CrossRef] [Google Scholar]
 Wilson, A. S., Smith, D. A., & Young, A. J. 2006, ApJ, 644, L9 [NASA ADS] [CrossRef] [Google Scholar]
 Wise, M. W., McNamara, B. R., Nulsen, P. E. J., Houck, J. C., & David, L. P. 2007, ApJ, 659, 1153 [NASA ADS] [CrossRef] [Google Scholar]
 Zhang, H., Li, H., Guo, F., & Taylor, G. 2017, ApJ, 835, 125 [NASA ADS] [CrossRef] [Google Scholar]
 Zhdankin, V., Uzdensky, D. A., Werner, G. R., & Begelman, M. C. 2019, Phys. Rev. Lett., 122, 055101 [Google Scholar]
 Zhdankin, V., Uzdensky, D. A., & Kunz, M. W. 2021, ApJ, 908, 71 [NASA ADS] [CrossRef] [Google Scholar]
Appendix A: Numerical integration
In this simulation, we numerically solve Eq. 1. The steps of numerical interaction are as follows:

Calculate conservation variables from principle variables at the end of the previous time step, and calculate the effective adiabatic index γ_{gas} via Eq. 15.

Adopt an operator split method (S = 0 in Eq. 1) and solve the conservative equation (see Matsumoto et al. (2019) for details).

Recalculate conservation variables, except for the gas thermal energy u_{gas}, from principle variables. Note: updated electron and proton temperatures are unknown, thus we cannot calculate u_{gas} and γ_{gas} at this step.

Solve the entropy equations for electrons and protons (Eqs. 8 and 9). The detailed procedure is described in the following section. Then, calculate u_{gas} and γ_{gas} using updated electron and proton temperatures.
A.1. Solving entropy equations
We describe the procedure of numerical integration for the entropy equations of electrons and protons. For simplicity, we describe the method for the 1D coordinate x. First, we need to estimate the energy dissipation at every grid for each time step. The dissipation energy is the difference between the thermal energy of the total gas, u_{gas}, which as obtained from the energy equation and the sum of the thermal energies that evolved by purely adiabatic evolution for electrons and protons (u_{e, ad} and u_{p, ad}). To compute adiabatic evolution, the righthand sides of Eqs. 8 and 9 are set to zero. We then solve these equations by the finitedifference method. We let x_{i} be the cell center of a uniform grid, Δx is the cell width. Time describes t^{n} = nΔt, where Δt is the time step. We adopted third order TVDRunge–kutta schemes, such that finitedifference equations are given as follows:
where (1) and (2) denote each number of subtime steps of the TVDRunge–kutta scheme. We note that the entropy formulas for electrons and protons are of the same form and we do not distinguish between them.
According to Sadowski et al. (2017), we arrange Eqs. A.1A.3 as follows:
where f and f′ are the fractions of the final state represented by the three contributing grids of gas. When two individual gases are mixed in a constant volume, the total energy is the sum of the initial energies of the gases. In contrast, the total entropy is not to be the sum of the initial entropy of the gases. Hence, the finitevolume methods in Eqs. A.4A.6 are incorrect. However, to overcome this problem, we must treat the dynamics of each gas individually. In this study, according to Sadowski et al. (2017), we solve Eqs. A.4A.6 by replacing the entropy with the thermal energy. For electrons, the relationship between the entropy and the dimensionless temperature is known in Eq. 11. Then, given the dimensionless temperature, we can calculate the adiabatic index by Eq. 14. Therefore, the electron thermal energy, u_{e} = n_{e}kT_{e}/(γ_{e}(θ_{e})−1), is a function of the density and entropy:
This process corresponds to assuming gas mixing at a constant gas density. To be safe, we used upwind values of the entropy, and . We also calculated the proton thermal energy following the same procedure described above, but protons are nonrelativistic in our simulation (see Eq. 13). From the above calculation, we can obtain the thermal energy of electron and proton at each RungeKutta substeps as
We recall that through above process, the effect of nonadiabatic process is ignored. Therefore, this is the internal energy that evolved by purely adiabatic evolution, hereafter denoted as u_{ad}.
Then, we calculate the dissipation heating at each grid as:
where dQ_{heat}/dt = q_{heat}. Further, if necessary, the fraction of the electron heating, f_{e}, is calculated using MHD quantities. Dividing the dissipation heat into electrons and protons, their thermal energies are updated as follows:
The source term, namely energy transfer via Coulomb and radiative cooling, is updated implicitly by adopting the Newton–Raphson iteration at last substep of the TVDRunge–Kutta scheme. Finally, we can easily recover the entropy s^{(1),(2),n + 1} ( and the temperature T^{(1),(2),n + 1} ) for each species using ρ^{(1),(2),n + 1} and u^{(1),(2),n + 1}.
A.2. Bremsstrahlung radiation cooling
The bremsstrahlung cooling rate per unit volume for relativistic plasma is given by(Svensson 1982)
where α_{f} and σ_{T} denote the finestructure constant and Thomson crosssection, respectively. In the above equation, F_{ei}, and F_{ee} are the dimensionless radiation rates due to protonelectron and electronelectron collisions, respectively. The approximation formulas of F_{ei}, and F_{ee} are respectively
Here, η_{E} = exp( − γ_{E}) and γ_{E} ≈ 0.5772 is Euler’s number.
A.3. Coulomb coupling
The rate of energy transfer, q^{ie}, from ions to electrons per unit volume through Coulomb collisions is determined as follows (Stepney & Guilbert 1983; Dermer et al. 1991);
where θ_{m} = θ_{p}θ_{e}/(θ_{p} + θ_{e}) and θ_{p} ≡ kT_{p}/m_{p}c^{2}. The parameters σ_{T}, and c denote the Thomson scattering crosssection and the speed of light, respectively. lnΛ is the Coulomb logarithm and approximated as:
for T > 4 × 10^{5} K (Spitzer 1962). Functions K_{0}, K_{1}, and K_{2} are respectively modified Bessel functions of the second kind of orders 0, 1, and 2.
Appendix B: Shockfinding algorithm
To identify whether the MHD grid is inside the shock zone, we implemented a shockfinder in MHD code CANS+ and adopted an approach similar to that followed by Ryu et al. (2003) and Schaal & Springel (2015). Although our method is based on the theory for hydrodynamic shock, the influence of omitting the magnetic field is insignificant for shockfinding. The inclusion of magnetic fields complicate the system to add two types of compressible shocks, and the measurement of the Mach number for a MHD shock is difficult task.
Here, we used the Cartesian coordinate (x, y, z), and subscript i ∈ (x, y, z) to denote the direction of each coordinates. We divide each grid in or out of the shock zone to employ the following criteria
where and ℳ_{min} are the estimated Mach number and a minimum Mach number. In the simulations, the divergence operator replaces the central differences, and we adopt ℳ_{min} = 1.3. The Mach number of each grid is estimated from the RankineHugoniot condition across shocks, which is given by (Pfrommer et al. 2017)
where y ≡ p_{gas,2}/p_{gas,1} and . Up and downstream quantities are denoted by subscripts 1 and 2, respectively. We determined the direction of shock propagation, d_{s}, in each grid using the temperature gradient:
Appendix C: Comparison of the electron heating models
We performed the axisymmetric MHD simulations to compare the results by using different electron heating models for MHD turbulence (H10 and K19, Howes 2010; Kawazura et al. 2019). The electrontoproton heating ratio of H10 is written by
where c_{1} = 0.92, c_{2} = 1.6T_{e}/T_{p}, and c_{3} = 18 + 5log_{10}(T_{p}/T_{e}) for T_{p} > T_{e}, while c_{2} = 1.2T_{e}/T_{p}, and c_{3} = 18 for T_{p} < T_{e}. H10 has similar behaviour with K19 for β_{p} < 1. For K19, the protontoelectron heating ratio saturate at ∼35 at highβ_{p}, while it monotonically increase with inverse β_{p} for H10. The simulation setups are the same as the model B, except for numerical resolution and the coordinate. We use the cylindrical coordinate (r, ϕ, z), and the grid size is Δ_{r} = Δ_{z} = 0.05 kpc, which is half the size of our model B.
In Fig. C.1, we show the distributions of the electron temperature using the different electron heating models (H10 and K19). We find that the distributions of the electron temperature do not differ in the shockedICM. Meanwhile, the electron temperature of the model using H10 is slightly higher than that of the model using K19 in the cocoon, except for in the sheath of the beam. Since the plasmaβ_{p} in the sheath region is much higher than 10, the protontoelectron heating ratio of H10 is significant low (see also Fig. C.2). The heating ratio for both models is sensitive for plasmaβ_{p}, and hence the electron is locationally heated up around the jet head. The ratio of the thermal electron energies, U_{e, H10}/U_{e, K19}, in the cocoon is 1.28 so that the choice of the subgrid models would not have a significant affect for the main results in the case of the model A and B in Section 5. On the other hand, in the case of model C, the radio power using H10 would be much lower than that using K19 since the cocoon is consist of the highbeta plasma. A comparison with the 3D simulation of the model B, the propagation velocity in a 3D simulation is faster than that in axisymmetirc simulation due to the dentistdrill effect (Scheuer 1974). This trend also found in previous threedimensional simulation (Perucho et al. 2019). Because the kink mode do not develop in the axisymetric condition, the electrons are only heated up at the internal shocks. Meanwhile, in the 3D case, the dissipation due to the development of the kink mode also is dominant heating source for the electrons (see detail in section 5.1).
Fig. C.1. Distribution of the electron temperature using H10 (left side of the panel) and K19 (right side of the panel) at t = 16.86 Myr. 
Fig. C.2. T_{e}/T_{p} − β_{p} histogram for regions in the cocoon using H10 (left) and K19 (right) at t = 16.86 Myr, respectively. The dashed line depicts the electron to proton temperature ratio corresponding to the equilibrium state for plasma β_{p}, as implied by the turbulence heating in Eq. 24 for H10 and K19. 
Appendix D: Magnetic dissipative structures
Figures. D.1 and D.2 shown magnetic dissipative structures for all models.
Fig. D.1. Three panels representing magnetic dissipative structures at early stages for all models. In each panel, we show volumerenders with the strength of J ⋅ E. 
All Tables
All Figures
Fig. 1. Number density (top) and gas pressure (bottom) profiles of initial ICM as a function of radius. Blue dots represent the jets number density and the jets gas pressure, respectively. The initial gas temperature of ICM is 5 KeV over the entire simulation domain. 

In the text 
Fig. 2. Slices (in the y − z plane) of number density distribution for model A, B, and C at t = 9.52, 9.94, and 13.02 Myr, respectively. The white lines represent contours of the zcomponent of velocity v_{z} = 0.5v_{jet} (0.15c). 

In the text 
Fig. 3. Twodimensional distribution maps of beam barycenter R_{fG} as a function of time for model A (top), B (middle), and C (bottom), respectively. Horizontal yellow solid lines depict the distance, l_{kink}, for the development of the external kink model, and vertical yellow solid lines are the time, τ_{kink}, at which the jets propagate to distance, l_{kink}. 

In the text 
Fig. 4. Slices (in the x − y plane at z = 80 kpc) of the distribution of the zdirection component of velocity for models A, B, and C. 

In the text 
Fig. 5. Slices (in the y − z plane) of the xdirection component distribution of the magnetic field for models A, B, and C, respectively. 

In the text 
Fig. 6. Probability distribution functions of the characteristic length scale parallel to the magnetic field in the cocoon where z > 40 kpc for model A (red dotted), B (blue solid), and C (green dashed) at the end of the simulations. We define the cocoon as grids with the electron temperature higher than 10^{8} K and z > 40 kpc. 

In the text 
Fig. 7. Slices (in the y − z plane) of the electron temperature distribution for models A, B, and C, respectively. 

In the text 
Fig. 8. Slices (in the y − z plane) of the ratio of proton to electron temperature distributions for models A, B, and C, respectively. 

In the text 
Fig. 9. Energies in the cocoon as a function of time for all models. Left: Time evolution of different energy components of the cocoon for model A (dotted lines), B (solid lines), and C (dashed lines), respectively. We define the cocoon as grids with an electron temperature higher than 10^{8} K. Right: Time evolution of the ratio between the magnetic field and the electron energy for model A (dotted lines), B (solid lines), and C (dashed lines). 

In the text 
Fig. 10. Relationship of three energy components: electron thermal energy, proton thermal energy, and magnetic field energy. Left: T_{e}/T_{p} − β_{p} histogram for regions in the cocoon for model A (a), model B (b), and model (c) at t = 9.52 − 9.94 and 13.02 Myr, respectively. The dashed line depicts the electron to proton temperature ratio corresponding to the equilibrium state for plasma β_{p}, as implied by the turbulence heating in Eq. (24). Right: Same as left panel, but displaying U_{e} − U_{mag} histogram. Dashed line plots U_{e} = U_{mag}. 

In the text 
Fig. 11. Projected distance from core to cavity center, R, versus Projected semiminor axis of the cavity, b. Black circles show the radiofilled cavities taken from Rafferty et al. (2006). Blue square shows our result for model B at t = 9.94 Myr. Because R and b values for models A and C have similar values to model B, we do not show them. 

In the text 
Fig. 12. Radio synchrotron power P_{radio} versus jet mechanical power estimated from the Xray cavity system: (adopted from Laura Bîrzan’s PhD thesis Rafferty 2007). Filled black symbols show radiofilled cavities (which include intermediate cases). Symbols and wide error bars denote the values of the mechanical power calculated using the sound speed. Open squares and filled squares show our results of cases 1T and 2T for models A (red), B (blue), and C (green), respectively. The diagonal dotted lines (dashed lines) represent ratios of constant mechanical power to radio luminosity. Blue dashed lines depict the total injection energy of our jet model, given by Eq. (20). 

In the text 
Fig. C.1. Distribution of the electron temperature using H10 (left side of the panel) and K19 (right side of the panel) at t = 16.86 Myr. 

In the text 
Fig. C.2. T_{e}/T_{p} − β_{p} histogram for regions in the cocoon using H10 (left) and K19 (right) at t = 16.86 Myr, respectively. The dashed line depicts the electron to proton temperature ratio corresponding to the equilibrium state for plasma β_{p}, as implied by the turbulence heating in Eq. 24 for H10 and K19. 

In the text 
Fig. D.1. Three panels representing magnetic dissipative structures at early stages for all models. In each panel, we show volumerenders with the strength of J ⋅ E. 

In the text 
Fig. D.2. Same as Fig. D.1, but at a later phase. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.