© The Authors 2023

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 jet-driven 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 non-thermal 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 X-ray 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_cavt_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 X-ray, as well as the volume of the cavity, respectively. Bîrzan et al. (2008, 2004) found that the mechanical power estimated from the X-ray 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 non-radiating 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 jet-ICM 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 non-axisymmetric modes, namely Kelvin–Helmholtz (Bodo et al. 1994), Rayleigh–Taylor (Matsumoto & Masada 2013), and current-driven kink modes (Mizuno et al. 2009; Mignone et al. 2010; Porth & Komissarov 2015). Thus, the non-linear evolution of these instabilities plays a significant role in the jet dynamics and their large-scale 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 two-temperature 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². 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 non-axisymmetric 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 two-temperature 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 two-temperature three-dimensional (3D) MHD simulations of semi-relativistic jets in galaxy clusters. In Sect. 2, basic equations and sub-grid 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 shock-finding algorithm.

2. Numerical method

2.1. Basic equations

Multiple methods have been developed to incorporate electron thermodynamics into single-fluid simulations self-consistently 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 single-temperature MHD code CANS+ (Matsumoto et al. 2019) to a two-temperature 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:

(1)

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:

(2)

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_en_e + m_pn_p ≈ m_pn. The total energy, E, and total pressure, p_T, are:

(3)

(4)

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:

(5)

The source terms are

(6)

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:

(7)

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:

(8)

(9)

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 trans-relativistic regime for electrons, and use the following approximate entropy formula derived by Sadowski et al. (2017) in Appendix A:

(10)

where θ_e ≡ kT_e/m_ec² is the dimensionless temperature. In this approximate formula, we can easily obtain the electron temperature for given s_e and ρ_e as:

(11)

The adiabatic index for electrons is calculated as follows:

(12)

In contrast, protons are non-relativistic in this simulation. Thus, we use the non-relativistic entropy formula for protons:

(13)

where γ_p = 5/3. The thermal energies of protons, electrons, and gas are as follows:

(14)

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):

(15)

2.2. Sub-grid models of electron heating

We considered two sub-grid 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 shock-finding 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.,

(16)

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 particle-in-cell (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:

(17)

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:

(18)

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 two-temperature 3D MHD simulations in Cartesian coordinates with the z-axis 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 xz-plane at z = 0. Other boundaries are imposed on the free-boundary condition.

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 X-ray data (Wilson et al. 2006; Smith et al. 2002). The density profile of ICM is given by:

(19)

where , n₀, 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.05 cm⁻³. 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 z-axis, 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 sub-relativistic flows to be consistent with the outburst energy of Cygnus A jets, 0.6 − 0.8 × 10⁴⁶ erg s⁻¹ (Bîrzan et al. 2004; Snios et al. 2018). The injected jet power can be written as follows:

(20)

We set the kinetic power as L_kin = 5.5 × 10⁴⁵ erg s⁻¹, and the thermal power as L_th = 4.4 × 10⁴⁴ erg s⁻¹.

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¹⁰ 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 small-amplitude (1%) random pressure perturbation for the injection flow is adapted to model non-axisymmetric features. We list common parameters for ICM and jets in Table 2.

The jets have a purely toroidal magnetic field B_ϕ = B_jetsin⁴(π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 matter-dominated 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 yz-plane 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 shocked-ICM, compressed by the bow shock, and the low-density cocoon are formed. We observe that the number density of the cocoons is very low, ∼10⁻⁴ cm⁻³, 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 non-axisymmetric 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 z-component of velocity vz = 0.5vjet (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),

(21)

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 two-dimensional 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 non-axisymmetric mode. We compute the beam barycenter as follows:

(22)

Fig. 3. Two-dimensional distribution maps of beam barycenter RfG as a function of time for model A (top), B (middle), and C (bottom), respectively. Horizontal yellow solid lines depict the distance, lkink, for the development of the external kink model, and vertical yellow solid lines are the time, τkink, at which the jets propagate to distance, lkink.

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 non-axisymmetric 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 xy-plane 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 low-magnetized 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 non-rotating 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 cross-section, ∼10 kpc, of positive velocity field (v_z > 0: from blue to yellow region), and finger-like structures. Owing to the high-mixing 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 z-direction component of velocity for models A, B, and C.

Figure 5 shows the slices of the x-direction 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 x-direction 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 small-scale turbulence is excited. To discuss the length-scale of the magnetic filament quantitatively, we evaluate the length-scales parallel to the magnetic field as (Schekochihin et al. 2004; Bodo et al. 2011; Mukherjee et al. 2020):

(23)

Figure 6 shows the volume-weighted probability distribution function (PDF) of the L_∥ in the jet cocoon where the electron temperature is higher than 10⁸ 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 108 K and z > 40 kpc.

4.2. Temperature distribution

Figure 7 shows the distribution of the electron temperature at the yz-plane 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 sub-grid 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 two-dimensional 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⁹ − 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 108 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

(24)

where we adopt γ_p = γ_e = 5/3. For models A and B, a quasi-steady 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: Te/Tp − β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 Ue − Umag histogram. Dashed line plots Ue = Umag.

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. Small-scale 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 non-thermal 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², 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 volume-renders 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 FR-II type jets. When model A enters the non-linear 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 small-scale 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 FR-I 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 two-temperature 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 X-ray 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:

(25)

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 power-law distribution for non-thermal electrons, dN/dγ = N₀γ^−s, where s, N₀, and γ are the electron energy index, number of non-thermal electrons, and Lorentz factor of non-thermal electrons. Meanwhile, the two-temperature MHD simulations account only for the evolution of thermal electrons. Because previous studies suggest that the non-thermal electron energy is proportional to the thermal gas energy or magnetic energy as a first-order 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 two-type model for approximation of non-thermal electrons as cases 1T and 2T. Case 1T is N₀ = C₀u_p, and case 2T is N₀ = C₀u_e. Here, , and η = 0.2 is a parameter, which is the ratio of the non-thermal electron energy density to the electron thermal energy density. Furthermore, we set γ_min = 100 and γ_max → ∞. Case 1T is corresponds to the single-temperature 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 z-axis, is 80°, and the observed frequency is 144 MHz.

We list the calculated radio powers in Table 3. The radio power of model A-1T, which is most prominent model, is two orders of magnitude higher than that of model C-2T. The 1T models have same order of non-thermal electron energy, because the proton temperatures have roughly the same value in lobes (see left panel of Fig. 9). Notably, 1T models assume that non-thermal 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 two-temperature models is sensitive to the magnetic field energy.

5.2.2. X-ray 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⁻¹. 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 X-ray observations.

(26)

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⁻¹⁰ erg cm⁻³ 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⁸ K. We confirm that the region whose electron temperature exceeds 10⁸ K has a lower density than ICM.

First, we mention the morphological properties of the X-ray 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 semi-minor 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 (kinetic-dominated) 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 low-density jet (Ohmura et al. 2021) and/or a low-power jet (Mukherjee et al. 2020). Otherwise, a long-periodic 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 large-scale 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 semi-minor axis of the cavity, b. Black circles show the radio-filled 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_gasV = 2.4 × 10⁵⁹, 2.6 × 10⁵⁹, and 3.7 × 10⁵⁹ 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⁻¹ (corresponding to ICM temperature T = 5 KeV). Therefore, we estimate the mechanical power P_cav = 1.5 × 10⁴⁴, 1.6 × 10⁴⁴, and 2.3 × 10⁴⁴ erg s⁻¹ 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 over-pressured 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 non-radiating proton thermal energy. Naively, protons can be energetically dominant over radiatively inefficient lobes (P_cav ≫ P_radio), as the pressure of non-radio emitting protons supports the expansion of the cocoon. Our jets are active during the simulation time. Therefore, we compare our results with radio-filled cavities, except for radio-ghost cavities. Although we plot the radio-filled 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 Pradio versus jet mechanical power estimated from the X-ray cavity system: (adopted from Laura Bîrzan’s PhD thesis Rafferty 2007). Filled black symbols show radio-filled 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 two-temperatures 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 protons-electrons 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 radiation-efficient lobe, an acceleration mechanism of η ≫ 10 is needed, where the energy of non-thermal 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 electron-proton 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 pair-plasma 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 electron-positron-proton mixture jets likewise achieved consistent results with observed FR-II 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 two-temperature simulations is necessary.

We discuss only the radio-fill cavities, as our jet is always active during the simulation time. However, several radio and X-ray observations imply multi-episodic 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 trans-relativistic 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 high-quality 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 inverse-Compton (CMB-IC) spectrum. If the energy of thermal electrons exceeds GeV, the thermal CMB-IC spectrum would be observed in X-ray 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 two-temperature 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 sub-grid scale of the radio lobes and, thus, two sub-grid electron heating mechanisms were considered. Because the jets have both the turbulence and the strong shock waves, we used the sub-grid 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 non-axisymmetric, current-driven kink mode. Meanwhile, weakly magnetized jets are decelerated by the high-mixing 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)

  Small-scale 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, shocked-electrons 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, high-temperature and high-magnetization multiple hotspots are formed in the beam.

• (vi)

  The low-density cavity of our model is narrow compared with observed radio-filled 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 X-ray 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 one-temperature approximation. Two-temperature 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 non-thermal 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 X-ray observations (e.g., Snios et al. 2018). We confirm that shocked-ICM plasma occurs in two-temperature states (see Fig. 8), hence, we aim to report thermodynamics and X-ray properties of shocked-ICM 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 post-K 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/Oakforest-PACS provided by JCAHPC through the HPCI System Research project (project id: hp190125, hp200092).

References

Appendix A: Numerical integration

In this simulation, we numerically solve Eq. 1. The steps of numerical interaction are as follows:

1. Calculate conservation variables from principle variables at the end of the previous time step, and calculate the effective adiabatic index γ_gas via Eq. 15.

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

3. 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.

4. 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 right-hand sides of Eqs. 8 and 9 are set to zero. We then solve these equations by the finite-difference method. We let x_i be the cell center of a uniform grid, Δx is the cell width. Time describes tⁿ = nΔt, where Δt is the time step. We adopted third order TVD-Runge–kutta schemes, such that finite-difference equations are given as follows:

(A.1)

(A.2)

(A.3)

where (1) and (2) denote each number of sub-time steps of the TVD-Runge–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.1-A.3 as follows:

(A.4)

(A.5)

(A.6)

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 finite-volume methods in Eqs. A.4-A.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.4-A.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_ekT_e/(γ_e(θ_e)−1), is a function of the density and entropy:

(A.7)

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 non-relativistic in our simulation (see Eq. 13). From the above calculation, we can obtain the thermal energy of electron and proton at each Runge-Kutta substeps as

(A.8)

(A.9)

(A.10)

We recall that through above process, the effect of non-adiabatic 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:

(A.11)

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:

(A.12)

(A.13)

The source term, namely energy transfer via Coulomb and radiative cooling, is updated implicitly by adopting the Newton–Raphson iteration at last sub-step of the TVD-Runge–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)

(A.14)

where α_f and σ_T denote the fine-structure constant and Thomson cross-section, respectively. In the above equation, F_ei, and F_ee are the dimensionless radiation rates due to proton-electron and electron-electron collisions, respectively. The approximation formulas of F_ei, and F_ee are respectively

(A.15)

(A.16)

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);

(A.17)

where θ_m = θ_pθ_e/(θ_p + θ_e) and θ_p ≡ kT_p/m_pc². The parameters σ_T, and c denote the Thomson scattering cross-section and the speed of light, respectively. lnΛ is the Coulomb logarithm and approximated as:

(A.18)

for T > 4 × 10⁵ K (Spitzer 1962). Functions K₀, K₁, and K₂ are respectively modified Bessel functions of the second kind of orders 0, 1, and 2.

Appendix B: Shock-finding algorithm

To identify whether the MHD grid is inside the shock zone, we implemented a shock-finder 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 shock-finding. 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

(B.1)

(B.2)

(B.3)

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 Rankine-Hugoniot condition across shocks, which is given by (Pfrommer et al. 2017)

(B.4)

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:

(B.5)

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 electron-to-proton heating ratio of H10 is written by

(C.1)

where c₁ = 0.92,  c₂ = 1.6T_e/T_p, and c₃ = 18 + 5log₁₀(T_p/T_e) for T_p > T_e, while c₂ = 1.2T_e/T_p, and c₃ = 18 for T_p < T_e. H10 has similar behaviour with K19 for β_p < 1. For K19, the proton-to-electron 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 shocked-ICM. 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 proton-to-electron 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 sub-grid 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 high-beta 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 dentist-drill effect (Scheuer 1974). This trend also found in previous three-dimensional 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. Te/Tp − β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 volume-renders with the strength of J ⋅ E.

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

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 z-component of velocity vz = 0.5vjet (0.15c).
In the text

	Fig. 3. Two-dimensional distribution maps of beam barycenter RfG as a function of time for model A (top), B (middle), and C (bottom), respectively. Horizontal yellow solid lines depict the distance, lkink, for the development of the external kink model, and vertical yellow solid lines are the time, τkink, at which the jets propagate to distance, lkink.
In the text

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

	Fig. 5. Slices (in the y − z plane) of the x-direction 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 108 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 108 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: Te/Tp − β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 Ue − Umag histogram. Dashed line plots Ue = Umag.

In the text

	Fig. 11. Projected distance from core to cavity center, R, versus Projected semi-minor axis of the cavity, b. Black circles show the radio-filled 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 Pradio versus jet mechanical power estimated from the X-ray cavity system: (adopted from Laura Bîrzan’s PhD thesis Rafferty 2007). Filled black symbols show radio-filled 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. Te/Tp − β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 volume-renders 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