Simulations of two-temperature jets in galaxy clusters I. Effect of jet magnetization on dynamics and electron heating

,


Introduction
Jets driven by the active galactic nuclei play a significant role in the 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 ICM from catastrophic cooling and falling into the host galaxy.(Fabian 1994;McNamara & Nulsen 2007;Fabian 2012).To better understand and discuss the heating energy of ICM in a quantitative manner, estimating the total kinetic energy of the jet is of high 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), i.e., 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 cav t age = 4pV, where P cav , t age , p, and V are the cavity power, source age, pressure of ICM observed by thermal X-ray, and volume of cavity, respectively.Bîrzan et al. (2008Bîrzan et al. ( , 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 cavity increases with radio luminosity, P cav ∝ P α radio , 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, P cav /P Radio ∼ 1.Although there are several physical factors of these scatters, such as electron cooling, estimation of ages, plasma composition, and variable activity of AGNs, their contribution in creating this scatter is not 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 pointed out that magnetohydrodynamics (MHD) instabilities develop non-axisymmetric modes, namely Kelvin-Helmholtz modes (Bodo et al. 1994), Rayleigh-Taylor modes (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(Ohmura et al. , 2020) ) focused on the two-temperature plasma, where the electrons and protons in jets are not in thermal equilibrium, because Coulomb coupling is inefficient in the tenuous plasma (Braginskii 1965;Stepney & Guilbert 1983).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. 2019Kawazura et al. , 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 M A and plasma-β β gas , are important for both electron heating and dynamical evolution.However, none of the studies considered two-temperature plasma, and thus 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 three-dimensional, two-temperature MHD simulations of semi-relativistic jets in galaxy clusters.In section 2, basic equations and sub-grid models for electron heating are presented.We describe the setup of our simulation in section 3, and the results are presented in section 4. In section 5, we discuss the observational implications of this study.A brief summary of key findings is provided in section 6.The appendix describes detailed numerical methods for solving entropy equations and the shock-finding algorithm.

Basic equations
Methods have been developed to incorporate electron thermodynamics into single-fluid simulations self-consistently in the con-text 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 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 read where ρ, v and B are the mass density, the bulk velocity, and the magnetic field, respectively.We assume the gas to hydrogen, such that n e ∼ n p , where n e,p are respectively the number density of proton and electron.Then, ρ = m e n e + m p n p ≈ m p n.The total energy E and total pressure p T are respectively 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 trans-relativistic 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 non-relativistic in this simulation.Thus, we use the non-relativistic entropy formula for protons: where γ p = 5/3.The thermal energies of protons, electrons, and gas are as follows: , u e = p e γ e (T e ) − 1 , u gas = u p +u e = p gas γ gas (T p , T e ) − 1 . ( 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) γ gas (T p , T e ) = 1+(γ e −1)(γ p −1) 1 + T p /T e (γ p − 1) + (γ e − 1)T p /T e .(15)

Sub-grid models of electron heating
We consider 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., f e (x, y, z) = f e,shock (for shock zone) f e,turb (for otherwise) (16) 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 model 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 Section 4.1 of Ohmura et al. (2020).
For electron to proton heating rates of MHD turbulence, there are two models proposed by (H10 Howes 2010) and (K19 Kawazura et al. 2019).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.

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.

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 where r = x 2 + y 2 + z 2 , 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 z-axis, and B ICM = 0.44 µG.The blue lines of figure 1 show the initial density (top panel) and pressure (bottom panel), respectively.Gas Pressure [erg/cc] 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.

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 observation indicates 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 M = 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 figure 1, the blue dots for each panel denote the density and pressure in the jet injection region, respectively.A small-amplitude (1 percent) 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 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 matter-dominated jets, i.e., the kinetic energy of jets exceeds the Poynting flux energy.

Results
We conducted simulations with various magnetic energies to investigate the effect of jet dynamics and electron heating.First, we focus on the effect of the magnetic field strength on the jet dynamics.The strength of the magnetic field affects the development of instability such as kink, Kelvin-Helmholtz, and Rayleigh-Taylor modes.The electron heating model for turbulence is the function of plasma-β p , and hence magnetization strongly affects the electron temperature distributions.Next, we report the time evolution for thermal electrons and protons in the jet lobe.

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 −4 cm −3 , such that Coulomb coupling is inefficient.The jet beam reaches jet tip despite suffering MHD instabilities (see red contours in figure 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.
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, 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 l kink ∼ v jet τ kink .For models A and B, the distances 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, and hence the 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 G , 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: R G (t, z) = x y rv z (x, y, z, t)dxdydz x y v z (x, y, z, t)dxdydz Similar analysis is 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.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).Further, for models A and B, the jets clearly separate between the beam flow (yellow region) and cocoon gas (red region), i.e., 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, that 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, Komissarov et al. (2019) found analytically that jets are stable for the Rayleigh-Taylor instability mode when M A < 40.Thus, the jet is in an unstable mode of the Rayleigh-Taylor instability for model C (see figure 2 and table 1).The right panel of figure 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.
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, i.e., 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 details in section 5.1.
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, 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) 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 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).

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 section 4.3).Here, we recall that the electron heating fraction for turbulence is proportional to inverse plasma beta β p −1 .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 Figure 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 ten times higher than electron temperatures in the cocoons for all models (figure 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.

Lobe energetics
The left panel of figure 9 displays the time evolution of different energy components of the cocoon for all models.Proton thermal energy is the dominant energy component of the cocoon (U p ≫ U e ), i.e., cocoons are supported by the proton pressure (see red and blue lines in the left panel of figure 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 β p −1 , 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 figure 9).
In the right panel of figure 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 ( 14)), and 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 where we adopt γ p = γ e = 5/3.For model Aand B, a quasisteady state of MHD turbulence is observed in the T e /T p −β p histogram (left panel of figure 10ab).This histogram plots the gas stored in the cocoons.We observe that the gas distribution follows along the dashed line, which is plotted as equation 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 equation 17).Another view on the turbulence equilibrium is the relationship between U e and U mag , shown in the right panel of figure 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, i.e., 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.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 be-tween the electron thermal energy and the magnetic field energy (see figure 10c).The volume occupations of the gas where 0.666 < U e /U mag < 3 for model C are 0.35.Further, 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 figure 10.However, the saturation value depends on the proton plasma beta when β p ≫ 10 in a cocoon.

Small-scale dissipation in jet beam
As we reported in section 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 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 in Appendix D displays the volume-renders of a physical quantity j • E at times that 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 figure 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 figure 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 (figure 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, and 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).

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

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: 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.
where j ν and ν are the synchrotron emissivity and the observed frequency, respectively.We use the formula of the synchrtoron emissitity provided in Pacholczyk (1970).
We assume a single power-law distribution for non-thermal electrons, dN/dγ = N 0 γ −s , where s, N 0 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 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 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, i.e., 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 degrees, and the observed frequency is 144 MHz.
We list the calculated radio powers in table 3 (see also figure 12).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 figure 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, larger than that of 1T models.This indicates that the radio power of two-temperature models is sensitive to the magnetic field energy.

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 observation.
where t age is the outburst age of jets.Lobes and ICM are approximately in the pressure equilibrium state.Thus, we use the initial ICM pressure around middle of lobe p gas ∼ 2.0×10 −10 erg cm −3 to calculate the cavity power (see figure 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 X-ray cavity.In figure 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 figure 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 numerous previous studies (e.g., Massaglia et al. 2016;Perucho et al. 2019;Mukherjee et al. 2020).To create a T. Ohmura and M. Machida: Simulations of two-temperature jets in galaxy clusters Fig. 8. Slices (in the y − z plane) of the ratio of proton to electron temperature distributions for models A, B, and C, respectively.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).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).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 model A, B, and C, respectively.Further, the sound crossing time t c = R/c s,ICM is 51.6 Myr, where we adopt that 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 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 ∼ 5.Even if we estimate the mechanical power using the simulation time, the injection energy of the jet is ∼ 10 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 un-der the assumption that the radio lobe and ICM have reached in the pressure equilibrium state.

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 figure 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 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 figure 12, the samples have a large scatter in the relationship, P cav /P radio ∼ 1 − 1000.
We find that our two-temperatures model explains 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 those 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 section 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 figure 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. 2019Zhdankin et al. , 2020)).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 twotemperature 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, and lead 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.

Observational implications
We showed that electrons remain in a trans-relativistic temperature, i.e., 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, the thermal electrons whose energies range with a few MeV scatter photons in the infrared and optical range (Enßlin & Sunyaev 2002).

Summary
We carried out two-temperature MHD simulations for the 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 study 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 use the sub-grid models at shock waves and at the turbulence in a hybrid manner.12. Radio synchrotron power P radio versus jet mechanical power estimated from X-ray cavity system P cav = 4pVt −1 c (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 equation 20.
The main findings achieved in this study are listed 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, which were induced by 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 β p −1 .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-fill 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 10 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).Further, 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 X-ray observations (e.g., Snios et al. 2018).We confirm that shocked-ICM plasma occurs in two-temperature states (see figure 8), and hence we aim to report thermodynamics and X-ray properties of shocked-ICM in paper II.
where (1) and (2) denote each number of sub-time steps of the TVD-Runge-kutta scheme.Note that the entropy formulas for electrons and protons are of the same form, and we do not distinguish between them.
According ot Sadowski et al. (2017), we arrange equations A.1 -A.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 finite-volume methods in equations 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 equations 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 equation 11.Then, given the dimensionless temperature, we can calculate the adiabatic index by equation 14.Therefore, the electron thermal energy, u e = n e kT e /(γ e (θ e ) − 1), is a function of the density and entropy: u n i,i±1/2 = u(s n i,i±1/2 , ρ (1) i ), u (1) i,i±1/2 = u(s (1) i,i±1/2 , ρ (2) i ), u (2) i,i±1/2 = u(s (2) i,i±1/2 , ρ n+1 i ).(A.7)

Fig. 2 .
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.white lines represent contours of the z-component of velocity v z = 0.5v jet (0.15c).

Fig. 3 .
Fig. 3. Two-dimensional distribution maps of beam barycenter R G 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 .
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.

Fig. 5 .
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.

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

Fig. 10 .
Fig. 10.Plots related to 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 equation 24.Right: Same as left panel, but displaying U e − U mag histogram.Dashed line plots U e = U mag .

Fig. 11 .
Fig. 11.Projected distance from core to cavity center R vs. Projected semi-minor axis of 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 for models A and C have similar values as model B, we do not show them.

Fig.
Fig.12.Radio synchrotron power P radio versus jet mechanical power estimated from X-ray cavity system P cav = 4pVt −1 c (adopted from Laura Bîrzan's PhD thesisRafferty 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 equation 20.

Table 1 .
Numerical Models Model β gas,jet M A B jet [µG] L x × L y × L z [kpc] N x × N y × N z

Table 2 .
Jets and ICM common setup parameters Meanwhile, because the jet does not suffer the kink mode, the barycenter for model C is within 1 kpc in simulation time, i.e., the jet propagates straight.

Table 3 .
Radio power and the amount of PdV work for the simulation results