AGN jet evolution simulation with GADGET4-OSAKA

TL;DR

This study uses GADGET4-OSAKA to simulate AGN jet evolution, analyzing how numerical methods and parameters affect jet dynamics, energetics, and morphology, providing benchmarks for SPH modeling.

Contribution

It systematically compares SPH and grid-based simulations of AGN jets, highlighting sensitivities to artificial viscosity and establishing benchmarks for SPH-based jet modeling.

Findings

Jet lobe growth follows self-similar scaling relations.

Jet size converges with resolution but is sensitive to artificial viscosity.

Energy partitioning shows significant deviations from idealized models.

Abstract

Active galactic nuclei (AGN) jets are powerful drivers of galaxy evolution, depositing energy and momentum into the circumgalactic and intracluster medium (CGM/ICM) and regulating gas cooling and star formation. We investigate the dynamics of jet evolution in the self-similar regime using the smoothed particle hydrodynamics (SPH) code GADGET4-Osaka, systematically vary jet-launching schemes, artificial-viscosity prescriptions, mass resolution, and jet lifetimes and compare the results with grid-based simulation. Our analysis combines quantitative diagnostics of jet size and energetics with detailed morphological and thermodynamic characterizations from slice maps and phase diagrams. We find that jet lobe growth follows analytic self-similar scaling relations and converges with resolution, but is highly sensitive to the choice of artificial viscosity. While the overall jet size tracks self-similar predictions, the partitioning of thermal and kinetic energy departs significantly from the idealized picture, reflecting enhanced dissipation and mixing, which is consistent with the jet propagation in grid-based simulations. These results establish robust benchmarks for SPH-based jet modeling, provide insight into the physical and numerical factors shaping jet--medium interactions, and lay the groundwork for future studies of AGN feedback in realistic galactic and cluster environments.

Tables

Table 1. Table 1: Summary of simulations in this study.

Fixed Parameters for all runs
Jet Setup			Simulation Configuration
$v_{\mathrm{jet}}$	$15000\mathrm{km}{s}^{-1}$		${\rho }_0$	$1.2\times {10}^{-26}g{\mathrm{cm}}^{-3}$
$P_{\mathrm{jet}}$	${10}^{46}\mathrm{erg}/s$		$T_0$	${10}^{7.2}K$
${\theta }_{\mathrm{jet}}$	$10\mathrm{deg}$		$P_0$	$2.66\times {10}^{-11}g/\left(\mathrm{cm}\cdot s^2\right)$
			$K_0$	$4.24\times {10}^{32}{\mathrm{cm}}^4/\left(g^{2/3}\cdot s^2\right)$
			$(L_x,L_y,L_z)$	$(300,300,1200)\mathrm{kpc}$
			$t_{\mathrm{sim}}$	$100\mathrm{Myr}$
Varied Parameters
Run Name	$t_{\mathrm{jet}}\left[\mathrm{Myr}\right]$	$m_{\mathrm{gas}}\left[{10}^5{M}_{\odot }\right]$	$L_{\mathrm{sep}}\left[\mathrm{kpc}\right]$	AV Model	Launching Scheme
FID	100	1.81	1	TDV, ${\alpha }_{\max }=2$	RG
Artificial Viscosity Variants
StV	100	1.81	1	StV, $\alpha =2$	RG
HVisc	100	1.81	1	TDV, ${\alpha }_{\max }=5$	RG
LLVisc	100	1.81	1	TDV, ${\alpha }_{\max }=0.5$	RG
Launching Scheme Variants
RR	100	1.81	1	TDV, ${\alpha }_{\max }=2$	RR
SR	100	1.81	1	TDV, ${\alpha }_{\max }=2$	SR
Mass Resolution Variants
Hi	100	0.93	0.8	TDV, ${\alpha }_{\max }=2$	RG
LoA	100	6.11	1.5	TDV, ${\alpha }_{\max }=2$	RG
LoB	100	14.48	2	TDV, ${\alpha }_{\max }=2$	RG
LoC	100	48.87	3	TDV, ${\alpha }_{\max }=2$	RG
LoD	100	115.84	4	TDV, ${\alpha }_{\max }=2$	RG
Jet Time Variants
T50	50	1.81	1	TDV, ${\alpha }_{\max }=2$	RG
T20	20	1.81	1	TDV, ${\alpha }_{\max }=2$	RG
Grid-Based Variants
MESH	See Section 3.4

Equations

$${\dot{M}}_{\mathrm{j}}=2\frac{P_{\mathrm{j}}}{v_{\mathrm{j}}^2}.$$

$$L_1=\frac{1}{{\theta }_{\mathrm{j}}}\sqrt{\frac{2}{\pi \rho }\sqrt{\frac{{\dot{M}}_j^3}{2P_{\mathrm{j}}}}}=\frac{2}{{\theta }_{\mathrm{j}}}\sqrt{\frac{P_{\mathrm{j}}}{\pi \rho v_{\mathrm{j}}^3},}$$

$$L_2={\left(\frac{P_{\mathrm{j}}^2\rho }{p^3}\right)}^{1/4}.$$

$$L_{\mathrm{j}}(t)={\left(\frac{P_{\mathrm{j}}{t}^3}{\rho }\right)}^{1/5}.$$

$$L_{\mathrm{j}}(t)=c_1{\left(\frac{P_{\mathrm{j}}{t}^3}{{\rho }_0{r}_0^{\beta }}\right)}^{1/(5-\beta )},$$

$$c_1={\left(\frac{A^4}{18\pi }\frac{({\gamma }_A+1)({\gamma }_L-1)(5-\beta {)}^3}{9\left({\gamma }_L+\frac{A^2}{2}({\gamma }_L-1)\right)-4-\beta }\right)}^{1/(5-\beta )}.$$

$$v_{\mathrm{H}}(t)=\frac{3c_1}{5-\beta }{\left(\frac{P_{\mathrm{j}}}{{\rho }_0{r}_0^{\beta }}\right)}^{1/(5-\beta )}{t}^{\frac{\beta -2}{5-\beta }}.$$

$$\mathrm{d}{t}_{\mathrm{j}}=\frac{m_{\mathrm{gas}}{v}_{\mathrm{j}}^2}{P_{\mathrm{j}}}.$$

$${\dot{M}}_{\mathrm{j}}=2\frac{P_{\mathrm{j}}}{v_{\mathrm{j}}^2}.$$

$$K=\frac{P_{\mathrm{th}}}{{\rho }^{\gamma }}=(\gamma -1)\frac{u}{{\rho }^{\gamma -1}}$$

$${\log }_{10}{K}_{\mathrm{t}}={\log }_{10}{K}_0+f_{\mathrm{t}}({\log }_{10}{K}_{\max }-{\log }_{10}{K}_0),$$

$$v_{\mathrm{H}}=\frac{P_{\mathrm{j}}}{2{\rho }_{\mathrm{A}}{v}_{\mathrm{j}}^2{S}_{\mathrm{w}}}\left(\sqrt{1+\frac{4{\rho }_{\mathrm{A}}{S}_{\mathrm{w}}{v}_{\mathrm{j}}^3}{P_{\mathrm{j}}}}-1\right),$$

$$|z-z_{BH}|>0.8L_J$$

$$\sum {\mathbf{F}}_{\mathrm{H}}={\dot{m}}_j({\mathbf{v}}_j-{\mathbf{v}}_{\mathrm{H}})+{\rho }_A{\bar{S}}_w({\mathbf{v}}_{\mathrm{A}}-{\mathbf{v}}_{\mathrm{H}})|{\mathbf{v}}_{\mathrm{A}}-{\mathbf{v}}_{\mathrm{H}}|=M_{\mathrm{H}}{\mathbf{a}}_{\mathrm{H}},$$

$${\bar{S}}_w\equiv \xi S_w,0<\xi \le 1,$$

$${\dot{m}}_j(v_j-v_{\mathrm{H}})-{\rho }_{\mathrm{A}}{\bar{S}}_w{v}_{\mathrm{H}}^2=0.$$

$$v_{\mathrm{H}}\simeq \sqrt{v_0{v}_j},v_0=\frac{{\dot{m}}_j}{{\rho }_{\mathrm{A}}{\bar{S}}_w}.$$

$$v_{\mathrm{H}}=\frac{v_0}{2}\left(\sqrt{1+\frac{4v_j}{v_0}}-1\right),$$

$${\rho }_j(v_j-v_{\mathrm{H}}{)}^2+p_j={\rho }_{\mathrm{A}}{v}_{\mathrm{H}}^2+p_{\mathrm{A}}.$$

$${\rho }_j=\frac{{\dot{m}}_j}{S_w(v_j-v_{\mathrm{H}})},$$

$$\frac{v_{\mathrm{H}}}{v_j}=\frac{\sqrt{{\eta }_{\rho }}}{1+\sqrt{{\eta }_{\rho }}},{\eta }_{\rho }\equiv \frac{{\rho }_j}{{\rho }_{\mathrm{A}}}.$$

$$\beta \equiv \frac{p_j-p_A}{{\rho }_{\mathrm{A}}{v}_j^2},{\eta }_{\rho }\equiv \frac{{\rho }_j}{{\rho }_{\mathrm{A}}}.$$

$$\frac{v_{\mathrm{H}}}{v_j}=\frac{\sqrt{{\eta }_{\rho }+\beta (1-{\eta }_{\rho })}-{\eta }_{\rho }}{1-{\eta }_{\rho }}$$

$$\Pi \equiv \frac{p_j-p_{\mathrm{A}}}{{\rho }_j{v}_j^2}=\frac{\beta }{{\eta }_{\rho }},$$

$$\frac{v_{\mathrm{H}}}{v_j}=\frac{\sqrt{{\eta }_{\rho }\left[1+\Pi (1-{\eta }_{\rho })\right]}-{\eta }_{\rho }}{1-{\eta }_{\rho }}.$$

Full text

††pubyear: 2026††pagerange: AGN jet evolution simulation with GADGET4-OSAKA–A

1 Introduction

Active galactic nuclei (AGN) are widely recognized as key drivers of galaxy evolution. By injecting vast amounts of energy and momentum into their surroundings, AGNs regulate star formation and shape the thermodynamic state of the circumgalactic and intergalactic media (CGM/IGM). AGN feedback is typically categorized into two complementary modes (see, e.g., Merloni and Heinz 2008; Heckman and Best 2014; Dubois et al. 2016). In the radiative (quasar) mode, rapid black-hole accretion powers luminous, quasi-isotropic winds (Feruglio et al., 2015; Tombesi et al., 2015). In the kinetic (maintenance or radio) mode, more modest accretion launches highly collimated jets that prevent hot gas from cooling and accreting onto the galaxy (Tortora et al., 2009). Together, these feedback channels establish the observed black hole–galaxy scaling relations and shape the star-formation histories of galaxies (Bower et al., 2006; Croton et al., 2006).

AGN jets are among the most powerful astrophysical phenomena, capable of impacting not only their host galaxies but also the surrounding CGM, intracluster (ICM), and even the IGM. These jets, composed of highly energetic particles — primarily electrons and protons — have been extensively studied through multi-wavelength observations (Croston et al., 2018; Turner et al., 2018). X-ray imaging has revealed vast cavities in galaxy clusters (see, e.g., Bîrzan et al. 2004; McNamara et al. 2005; Wise et al. 2007; Hlavacek-Larrondo et al. 2012; Shin et al. 2016; Liu et al. 2019), interpreted as relics of jet-driven outflows, while radio observations trace synchrotron emission from relativistic particles accelerated within the jets (Burke-Spolaor et al. 2017; Gendron-Marsolais et al. 2017; Di Gennaro et al. 2018). Additional indirect signatures, such as X-ray-bright rims around radio lobes (Fabian et al. 2000; see also Prunier et al. 2025 for simulation comparison), reinforce the view that AGN jets can dramatically reshape their environment on both galactic and extragalactic scales.

As AGN jets propagate through the ambient medium, they interact dynamically with the ambient gas, altering the physical properties of both the jets and the surrounding medium. These interactions play a key role in the redistribution of energy, leading to the heating of the CGM/IGM (e.g., Yamada and Fujita, 2001; McNamara and Nulsen, 2007; Dong et al., 2023) and driving large-scale outflows (Lister and Homan, 2005; Guillard et al., 2012; Morganti et al., 2015; Venturi et al., 2021). Such processes are central to galaxy evolution, especially in massive elliptical galaxies, where they suppress star formation by inhibiting gas cooling (McNamara et al., 2006; Guillard et al., 2015). Over cosmic times, this feedback mechanism has had a profound impact on the growth of galaxies and the assembly of large-scale structures. In addition to the galaxy cluster study, the AGN jet also provides a plausible mechanism for the formation of the Fermi bubble around the Milky Way (Cheng et al., 2011; Guo and Mathews, 2012; Yang et al., 2012).

Rapid advancements in cosmological simulations have led to increasingly sophisticated methodologies for modeling the large-scale effects of radio-mode AGN feedback. Early sub-grid models typically employed isotropic “thermal bombs” or broad-cone momentum kicks (e.g., Di Matteo et al., 2005; Springel et al., 2005; Dubois et al., 2013). While effective at quenching star formation, these approaches failed to reproduce the elongated X-ray cavities and narrow radio lobes observed. More recent efforts have introduced collimated kinetic jets, with parameters such as opening angle, duty cycle, and mechanical efficiency calibrated either from small-scale magneto-hydrodynamic (MHD) simulations or from empirical scaling relations (e.g., Weinberger et al., 2017; Bourne and Sijacki, 2017; Cielo et al., 2018; Davé et al., 2019; Perucho et al., 2019). These studies demonstrate that the details of the jet-injection scheme — including mass loading, nozzle geometry, artificial viscosity, and numerical resolution — strongly influence lobe length, morphology of bow shock, and level of metal mixing.

A systematic comparison of these numerical ingredients remains lacking. Notablely, Huško and Lacey (2023) utilizing the smoothed-particle hydrodynamics (SPH) code SWIFT has demonstrated that self-similar jet growth can be reproduced when physical parameters align with theoretical predictions. Their work was conducted under the following specifications: (i) the use of a pre-allocated reservoir of jet particles, (ii) a time-dependent artificial viscosity (AV) model, and (iii) ignoring all physics except hydrodynamics. However, it remains unclear whether these findings can be extended to other modern SPH formulations when the treatment of dissipation and kernel bias varies. Furthermore, the sensitivity of jet lobe evolution to factors such as jet lifetime and gas mass resolution remains only partially understood, despite their relevance to interpreting radio and X-ray observations of relic AGN bubbles.

In this paper, we follow the reservoir-based kinetic-jet setup of Huško and Lacey (2023) (hereafter H23) and implement the method in another modern SPH code GADGET4-Osaka. Our objectives are threefold:

1. Demonstrate the numerical reproducibility of self-similar lobe growth under controlled conditions;

2. Quantify the sensitivity of jet evolution to key numerical and physical setups;

3. Provide benchmark simulations for future model validation of jet feedback implementations.

In this work, we adopt the homogeneous circumgalactic configuration in H23 and a similar particle-reservoir setup for jet launch. In addition to variations in mass resolution, we investigate how viscosity models, jet reservoir arrangement, and jet active time affect jet propagation.

The remainder of the paper is structured as follows. Section 2 reviews the self-similar analytic model that underpins our interpretation of lobe growth. Section 3 describes the implementation of GADGET4-Osaka, including the jet injection algorithm and the diagnostics used to identify jet lobes, and bow shocks. Section 4 presents our results, beginning with the fiducial simulation and following comparisons across all parameter variations. In Section 5, we place these findings in the broader context of previous numerical and observational studies. Finally, Section 6 summarizes our main conclusions and discusses future directions, including the incorporation of additional physics.

2 Analytical model

The jet structure derived from our fiducial simulation is illustrated in Figure 1. Following an outside-in order, the jet structure can be broadly decomposed into three primary components: the bow shock, the jet lobe, and the jet material. The bow shock predominantly consists of shocked ambient gas that has been compressed as the jet propagates, resulting in an elevated density and temperature compared to the unperturbed background. The outer boundary of the shocked ambient gas is marked by a bow shock that separates this shocked region from the undisturbed ambient medium. Within the bow shock, the jet lobe manifests as a cavity of low-density, high-temperature gas resulting from the interaction of the jet with the ambient medium. This region consists of a turbulent mixture of shocked jet material and entrained ambient gas. The tip of the jet lobe, which is known as the ‘jet head’, is the site of jet material mixing with the ambient gas. The jet (material) itself remains collimated along the symmetry axis, terminating at a strong shock near the jet head, where much of its kinetic energy is transferred to the lobe and bow shock.

It is essential to compare our simulated jet lobes with analytical models of self-similar jet evolution (e.g., Falle, 1991; Kaiser and Alexander, 1997; Komissarov and Falle, 1998; Begelman and Cioffi, 1989; Bromberg et al., 2011). In this framework, the jet is launched with a specific half-opening angle, $\theta_{\mathrm{j}}$, and is characterized by key physical parameters including the jet power ($P_{\mathrm{j}}$), launch velocity ($v_{\mathrm{j}}$), and the density ($\rho$), pressure ($p$), and temperature ($T$) of the ambient medium. For a purely kinetic (non-relativistic) jet, the mass-loading rate is determined by the balance between kinetic energy flux and power, yielding

[TABLE]

Here, we neglect relativistic corrections given the galactic scale of simulation and assume symmetry between the two bipolar jets, so that the quoted parameters describe the combined effect of both lobes.

The two characteristic length scales governing the jet propagation (Falle, 1991; Kaiser and Alexander, 1997; Komissarov and Falle, 1998) are given by

[TABLE]

[TABLE]

Here, $L_{1}$ marks the distance at which the mass swept up by the jet becomes comparable to the mass injected by the jet itself. In the early stages, when the jet length satisfies $L_{\mathrm{j}}\ll L_{1}$, the jet remains much denser than the ambient medium and thus propagates nearly unimpeded. This phase, referred to as ‘ballistic’ in H23, is characterized by linear growth of the jet length with time, $L_{\mathrm{j}}=v_{\mathrm{j}}t$.

The scale $L_{1}$ sets the transition from this freely expanding regime to one in which the internal jet pressure begins to balance the pressure of the surrounding shocked gas, leading to partial collimation. In contrast, $L_{2}$ marks the point at which the internal pressure of the lobe approaches equilibrium with the ambient pressure. Once the jet extends beyond $L_{2}$, its expansion slows as the ambient pressure becomes dynamically significant and the advance speed of the bow shock approaches the sound speed of the external medium (Komissarov and Falle, 1998; Huško and Lacey, 2023). Together, these two characteristic scales delineate whether the system evolves in a momentum-driven or pressure-confined regime, and provide essential benchmarks for interpreting the behavior of AGN jets in both simulations and observations.

In the intermediate regime defined by $L_{1}\ll L_{\mathrm{j}}\ll L_{2}$, the jet evolution is expected to follow a self-similar solution. In this phase, both the external pressure and the swept-up mass flux are dynamically negligible, leaving jet power and ambient gas density as the only relevant parameters. Since no intrinsic length or time scale can be derived from these two quantities alone, dimensional analysis implies that the system evolves self-similarly (Sedov, 1959).

In a uniform medium, the jet length grows as the density of the jets increases:

[TABLE]

For a more general ambient medium following a power-law density profile, $\rho(r)=\rho_{0}(r/r_{0})^{-\beta}$, the self-similar jet length evolves as Kaiser and Best (2007)

[TABLE]

where $c_{1}$ is a dimensionless constant defined by

[TABLE]

Here, $A=R_{\ell}/L_{\ell}$ is the aspect radio of lobe; $R_{\ell}$ represents the radius of the lobe, while $L_{\ell}\approx L_{J}$ is the length of the lobe. $\gamma_{A}$ and $\gamma_{L}$ are the adiabatic indices of the ambient gas and the lobe, respectively. In this work, we adopt $\gamma_{A}=\gamma_{L}=5/3$, as we neglect relativistic effects.

Taking the time derivative of Equation (5) yields the self-similar advance speed of the lobe head:

[TABLE]

This framework provides a theoretical baseline against which we can compare our simulated jets to assess the degree of self-similarity and the influence of numerical effects.

3 Numerical method

In this study, we employ the GADGET4-Osaka code (Romano et al., 2022a, b; Oku et al., 2022; Oku and Nagamine, 2024), a form of the SPH code GADGET-4 (Springel et al., 2021), to perform our hydrodynamical simulations. Following the methodology of H23, we include only hydrodynamic interactions among SPH particles and deliberately omit additional physics such as gravity, radiative cooling, magnetic fields, and cosmic rays. This stripped-down setup enables a clean examination of jet propagation in an idealized, homogeneous medium, thereby allowing direct comparison with analytic models and isolating numerical effects from physical complexities.

We describe the SPH methodology implemented in GADGET4-Osaka in Section 3.1, followed by an overview of the simulation setup in Section 3.2. Variations in jet launch schemes are discussed in Section 3.3, and the full suite of feedback variants explored in this study is summarized in Section 3.4. We also introduce comparisons with grid-based simulations by Mizuta et al. (2006); Mizuta and Ioka (2026) in Sec 3.4. Finally, the diagnostics used to quantify jet propagation are detailed in Section 3.5.

3.1 Numerical Schemes

3.1.1 SPH kernel formulation

We adopt the Wendland C4 kernel (Dehnen and Aly, 2012) with bias correction for all SPH interactions, setting the number of neighbors per particle to $n_{\mathrm{ngb}}=200$. This choice offers a good compromise between stability and resolution, while suppressing pairing instability.

Our simulations employ the pressure-entropy formulation of SPH introduced by Hopkins (2013), which offers significant advantages in handling contact discontinuities and shock fronts. The pressure-based formulation is known to better capture fluid mixing and to mitigate numerical surface-tension effects that can artificially suppress Kelvin-Helmholtz and Rayleigh-Taylor instabilities. Therefore, such a formulation becomes the default choice in the Osaka feedback models (Oku and Nagamine, 2024). This choice contrasts with H23, which uses an energy–density formulation, whereas GADGET4-Osaka does not.

Although gravitational effects are not incorporated within our simulations, we still adopted the Jeans pressure floor following the prescription of Hopkins et al. (2011) and Kim et al. (2016). Such a pressure floor is useful for maintaining numerical stability in the scenario of jet simulations characterized by highly supersonic flow dynamics.

3.1.2 Artificial Conductivity and Artificial Viscosity

GADGET4-Osaka incorporates both artificial conductivity and artificial viscosity. Artificial conductivity is crucial for mitigating density and temperature disparities that occur at shock interfaces, which is included with the velocity-based conductivity signal speed (Wadsley et al., 2008). To enhance the resolution of shock waves, a time-dependent limiter (Borrow et al., 2022) is applied to the artificial conduction in these regions.

Artificial viscosity (AV) is critical for resolving shocks and handling velocity discontinuities. In this study, we examine two AV formulations: the standard viscosity (StV) and time-dependent viscosity (TDV) schemes. Both are implemented in Gadget-4 (Springel et al., 2021). The standard viscosity has a constant viscosity coefficient for all SPH particles throughout the simulation, which is intended to capture the shocks around the supersonic flows like the AGN jet. The TDV formulation implemented with the Cullen–Dehnen switch (Cullen and Dehnen, 2010; Hu et al., 2014), by contrast, dynamically adjusts the viscosity coefficient on the basis of local flow conditions, reducing the numerical dissipation away from shocks while preserving the accuracy in supersonic regions.

To further improve the performance of AV in shock-capturing, we apply the higher-order velocity gradient estimator proposed by Hu et al. (2014) and get the particle relative velocity from velocities linearly reconstructed at the particle midpoint (Frontiere et al., 2017; Rosswog, 2020) with the Balsara slope limiter (Balsara, 2004) and the minmod limiter, which sharpens discontinuities and minimizes spurious angular momentum transport. Additionally, we adopt the Balsara switch (Balsara, 1995) to suppress AV in shear-dominated flows, preventing excessive damping of turbulence and vorticity.

3.1.3 Time Integration

The hydrodynamical force is integrated using the second-order predictor-corrector method (Springel et al., 2021) with the timestep limiter (Saitoh and Makino, 2009; Durier and Dalla Vecchia, 2012; Oku and Nagamine, 2024). This ensures proper adjustment of the particle timestep during jet propagation. We adopted a Courant number of $0.3$ in all our simulations.

3.2 Simulation Setup

We follow the initial condition setup by H23. The simulation is performed in a cuboid periodic simulation box with a side length of $(L_{x},L_{y},L_{z})=(300,300,1200)\,\mathrm{kpc}$, where the $z$ axis is extended to accommodate the direction of jet propagation. We emphasize to the readers that the periodic configuration is implemented solely for facilitating the accounting of the energy budget; the jet structures, however, do not traverse the periodic boundary upon the end of the simulation. The total simulation time is set to $T_{\mathrm{sim}}=97.8\,\mathrm{Myr}\approx 100\,\mathrm{Myr}$.

To initialize a uniform and static ambient medium, all gas particles are assigned a constant density of $\rho_{0}=1.2\times 10^{-26}\,\mathrm{g/cm^{3}}$ and a temperature $T_{0}=10^{7.2}\,\mathrm{K}$. The particles are placed on a regular grid and given zero initial velocity to ensure a static condition. The interparticle spacing is determined by the mass resolution. In our fiducial run (FID), the spacing is $L_{\mathrm{sep}}=1\,\mathrm{kpc}$, corresponding to a gas particle mass of $m_{\mathrm{gas}}=1.81\times 10^{5}\,\mathrm{M_{\odot}}$. In the resolution-variant runs (Hi, LoA, LoB, LoC, LoD; see 1), this spacing is adjusted according to the respective mass resolution to explore the effects of resolution on jet propagation.

3.3 Jet launching mechanisms

A key numerical challenge in modeling AGN jet feedback using SPH-based codes arises when injecting energy at a constant power. Continuous energy injection eventually depletes the gas supply in the vicinity of the black hole, leading to the formation of an unphysical low-density cavity around the jet origin (Vogelsberger et al., 2013; Barai et al., 2014, 2016; Weinberger et al., 2017). This artifact stems from the repeated consumption of local gas particles, prompting the black hole to increasingly draw from more distant neighbors, thereby disrupting the local gas structure.

This issue becomes particularly pronounced when the jet operates on short time steps, accelerating the depletion process. Additionally, if launched particles are not properly tagged or tracked, there is a risk that the same particle may be “kicked” multiple times, further exacerbating inaccuracies in the simulation of jet propagation and energy deposition.

In this study, we adopt a reservoir-based jet launching scheme in which a pre-allocated set of gas particles is designated at the initial conditions for subsequent injection. The total number of reservoir particles is determined by total energy of the jet $E_{\mathrm{jet}}=P_{\mathrm{jet}}T_{\mathrm{jet}}$ and the energy of one jet launching event $E_{\mathrm{launch}}=m_{\mathrm{gas}}v_{\mathrm{jet}}^{2}/2$. In our fiducial run, we set the jet-launching duration equal to the simulation time, $T_{\mathrm{jet}}=T_{\mathrm{sim}}=97.8\,\mathrm{Myr}\approx 100\,\mathrm{Myr}$. For comparison, we also examine shorter-duration scenarios with $T_{\mathrm{jet}}=0.5T_{\mathrm{sim}}\approx 50\,\,\mathrm{Myr}$ (T50) and $T_{\mathrm{jet}}=0.2T_{\mathrm{sim}}\approx 20\,\,\mathrm{Myr}$ (T20), in which the jet is active only during the early phase of the simulation.

This reservoir method, originally proposed by H23, offers a computationally efficient alternative to on-the-fly particle selection by avoiding repeated neighbor searches. It is particularly well-suited to our goal of studying jet propagation in a controlled environment. However, in more realistic simulations of self-regulated black hole accretion, reservoir particles would need to be spawned dynamically in the vicinity of the black hole, which would require additional implementation beyond our current GADGET4-Osaka setup.

In this idealized simulation, the AGN jet is launched at a fixed power $P_{\mathrm{j}}=10^{46}\,\mathrm{erg/s}$. Given our choice of a constant launching velocity $v_{\mathrm{j}}=15000\,\mathrm{km/s}$, the time interval between successive jet events, i.e., the injection of a pair of gas particles with mass $m_{\mathrm{gas}}$, is given by

[TABLE]

This expression represents the time required for the AGN to accumulate sufficient kinetic energy to launch one jet event. As a result, $\mathrm{d}t_{\mathrm{j}}$ depends directly on the mass resolution of the simulation. The initial temperature and density of the launched gas are set equal to those of the ambient medium, ensuring consistent thermodynamic conditions at injection.

We find that the spatial configuration of the reservoir particles can also influence jet morphology and propagation. To explore this, we implement three different particle placement schemes: reservoir grid (RG), reservoir random (RR) and Spawning Random (SR). We illustrate the launching schemes in Figure  2.

Reservoir Grid (RG) — In this scheme, the reservoir particles are placed on a Cartesian lattice and then trimmed into a biconical volume defined by the jet’s half-opening angle, $\theta_{\mathrm{j}}$, out to a maximum distance along the $z$-axis. Particles are injected in an ‘outside-in’ order, so the first particles launched lie closest to the bow shock and interact immediately with the ambient medium. Within each radial layer, the launch order is randomised to avoid artificial symmetry. We stress that ‘grid’ refers only to the initial spatial arrangement of reservoir particles; the hydrodynamics remain fully Lagrangian, and no fixed mesh is introduced in the simulation.

Reservoir Random (RR) — Like the RG scheme, the RR model relies on a pre-allocated set of particles, but their initial positions are chosen at random inside the biconical launch volume. This stochastic placement yields a more isotropic sampling of directions, thereby avoiding the drop in particle density toward the axis inherent to a trimmed grid.

In contrast, the trimmed-grid approach of RG reduces the number of particles per layer as one moves toward the center, resulting in a less uniform distribution within the specified opening angle. While randomization in the RR model can enhance spatial uniformity, it also introduces the risk of particles being positioned too closely together, which may affect the simulation’s efficiency and duration.

Spawning Random (SR) — In this scheme, jet particles are ‘spawned’ on the fly rather than drawn from a pre-defined reservoir. At each jet event, we inject one particle pair in opposite directions; the launch vector is chosen randomly within the prescribed half-opening angle, and the injection point is fixed at a distance of 5 kpc from the black-hole position.

A spawning approach is particularly attractive for large cosmological runs that track many black holes, whose jet duty cycles and particle budgets are not known a priori. Nevertheless, the spawned SPH particles in this setup are still pre-allocated in memory because of the technical challenges associated with introducing new gas particles in Gadget-4. We leave the integration of our jet-launching scheme into cosmological simulations to future work, which involves collecting and reusing the memory of removed particles.

All three launching schemes introduce additional mass into the simulation 111Note that, for all three schemes, jet particles are assigned to the computer memory ahead of time, but are not considered in the force calculation before they are launched. , set by the jet’s mass-loading rate,

[TABLE]

In the fiducial run, this yields $\dot{M}_{\mathrm{j}}=141M_{\odot}/\mathrm{yr}$, or $1.41\times 10^{10}M_{\odot}$ over the 100 Myr simulation, which is less than $0.1$ per cent of the total ambient gas mass in the initial condition. Because self-gravity is disabled in these idealised runs, such a modest mass surplus has a negligible dynamical impact.

Prior to launch, reservoir particles are ‘frozen’: they neither interact hydrodynamically nor appear in neighbour lists of active gas particles. The total number of frozen particles equals the expected number of jet events during the prescribed jet-active time, $t_{\rm jet}$. Each particle is activated only when its turn to be launched is reached, at which point it acquires the jet velocity and full hydrodynamic coupling.

3.4 Simulation Variants

To test the robustness and flexibility of our model, we conducted 11 simulations in addition to the fiducial run; their settings are summarised in Table 1. Each variant isolates one numerical or physical ingredient:

• •

Artificial-viscosity variants (StV, HVisc, LLVisc) — The StV employs a standard constant-viscosity model to evaluate the effect of viscosity on jet propagation. The HVisc, and LLVisc runs are still based on the time-dependent viscosity model, but we varied the max allowed viscosity coefficient ($\alpha_{\mathrm{max}}=5,0.5$) to investigate the impact of this parameter.

• •

Launching-scheme variants (RR, SR) — alter the way jet particles are introduced as shown in Section 3.3: RR randomizes the order in which reservoir particles are injected, whereas SR spawns new particles on-the-fly, sampling directions stochastically. These runs test the sensitivity of jet morphology to the injection geometry.

• •

Mass-resolution series (Hi, LoA–LoD) — spans nearly two orders of magnitude in particle mass, with each step differing by a factor of $\sim 3$, allowing us to evaluate how resolution affects lobe growth and bow shock structure in cosmological simulation contexts.

• •

Jet-lifetime variants (T50, T20) — truncate the period of active injection to approximately 50 Myr and 20 Myr, respectively, in order to track the evolution of ‘fossil’ lobes once the power source is switched off.

Lastly, we perform a simulation (hereafter the MESH run) using a three-dimensional grid-based hydrodynamic code (Mizuta and Ioka, 2026), which builds upon the methods described in Mizuta et al. (2006) and is configured to replicate the physical setup of our fiducial run as closely as possible. Because such hydrodynamic simulations are computationally expensive, we set up the computational domain in spherical coordinates $(r,\theta,\phi)$ with logarithmic spacing in the radial direction. Exploiting the symmetry of the system, we simulate only half of the sphere by restricting the azimuthal domain to $-\pi/2\leq\phi\leq\pi/2$ while keeping the full polar range $0\leq\theta\leq\pi$. The jet axis is placed at $(\theta,\phi)=(\pi/2,0)$, corresponding to the Cartesian $x$-axis, in order to avoid numerical issues associated with the polar coordinate singularity. However, for visualization purposes, we transform the coordinates so that the jet is aligned with the $z$-axis, ensuring consistency with the coordinate system adopted in our SPH simulations. The radial computational domain spans from $r=1.78\times 10^{22}\,\mathrm{cm}=5.78\,\mathrm{kpc}$ to $r=4.62\times 10^{24}\,\mathrm{cm}=1.50\times 10^{3}\,\mathrm{kpc}$ to align with the SPH setup. Reflective boundary conditions are imposed at the inner radial boundary as well as at the angular boundaries in $\theta$ and $\phi$. An outflow (free) boundary condition is applied at $r=r_{\rm max}$. The mesh resolution is given by $\Delta r=r\Delta\theta$, with $\Delta\theta=\Delta\phi=0.5^{\circ}$, ensuring approximately isotropic resolution near the jet axis. At the jet launching radius ($r=10\,\mathrm{kpc}$), this corresponds to a spatial resolution of 0.087 kpc. The total number of grid cells is $(N_{r},N_{\theta},N_{\phi})=(640,360,360)$. We adopt an ideal-gas equation of state with a constant adiabatic index $\gamma=5/3$. Although the original code was developed for relativistic jet simulations, it is employed here in the Newtonian regime.

3.5 Measurements

When benchmarking our simulations against analytic models, it is essential to adopt clear, internally consistent measurement definitions. Because we have complete knowledge of the three-dimensional gas distribution and its time evolution, we can construct precise, physically motivated diagnostics. Throughout this study, we employ uniform criteria to identify different structures (jet lobes, bow shocks) of the jet. Using these rigorous, simulation-specific definitions ensures that comparisons with theoretical scalings are both meaningful and reproducible.

3.5.1 Jet lobe

In our simulations, we identify jet-lobe as gas that is simultaneously (i) less dense than the initial ambient medium and (ii) hotter than that medium. In practice, the lobe is most cleanly isolated in entropy space, where it appears as a high-entropy tail well separated from the shocked ambient gas. For a given gas particle, we define the entropy ($K$) with internal energy per unit mass ($u$), mass density ($\rho$), and the adiabatic index $\gamma=5/3$ (Springel, 2005)

[TABLE]

where $P_{\mathrm{th}}=(\gamma-1)\rho u$ is the thermal pressure. We then adopt an entropy threshold

[TABLE]

where $K_{0}=4.24\times 10^{32}\,\mathrm{cm^{4}\cdot g^{-2/3}\cdot s^{-2}}$ is the initial ambient entropy, $K_{\mathrm{max}}$ is the maximum entropy in the domain at the time of measurement, and $f_{\mathrm{t}}=0.5$ sets the threshold halfway (in log-space) between the ambient value and the absolute maximum.

3.5.2 Bow shock

As the jet propagates into the surrounding medium, it drives a strong shock into the initially static ambient gas. Gas swept up by the shock is first accelerated to supersonic speeds, then decelerates to transonic and finally subsonic flow as it expands away from the jet head. The swept-up layer forms a bow shock that envelops the jet lobe and its surrounding gas.

In our diagnostics, this bow shock is most clearly visible in Mach-number maps, where it appears as a sharp jump in $\mathcal{M}$ from the essentially stationary ambient gas to the transonic medium. We delineate the shock front by selecting all gas with $\mathcal{M}>0.0316\quad(|v|>10^{-1.5}\,c_{\mathrm{s},0})$.

For those components, their lengths are defined as half the distance between the two most distant gas particles in opposite directions along the jet axis that satisfy the entropy criterion. The radius measurement is more sensitive to turbulence along the volume boundary. To obtain a stable estimate, we compute the mean cylindrical radius of the outermost 10 per cent of particles, measured relative to the jet axis. This procedure filters out small-scale corrugations and yields a robust measure of the radial extent. For our MESH run, we used a similar criteria; however, the size measurement is derived from pixel masks taken from a slice in the $x$-$z$ plane.

4 Results

We begin by presenting the results of our fiducial simulation run (FID) in Section 4.1, followed by a systematic comparison with the suite of simulation variants in Section 4.2.

4.1 Fiducial jets

This section analyzes the fiducial jet simulation, denoted as FID in Table 1. The set of adopted parameters yields characteristic length scales of $L_{1}=32.9\,\mathrm{kpc}$ and $L_{2}=915\,\mathrm{kpc}$, which define the initial and asymptotic regimes of the jet evolution, respectively. These serve as key reference points for interpreting the dynamics, energetics, and morphology of the jet and its bow shock throughout the simulation.

4.1.1 Time evolution of the jet morphology

Figure 3 presents the time evolution of the gas temperature, density, and entropy slices for the fiducial (FID) run. At $t=20\,\mathrm{Myr}$, the jet reaches a height of approximately 100 kpc, and the morphological features outlined in Figure 1 become evident. The central jet lobe exhibits a low-density, low-temperature, and low-entropy zone near the launch site, surrounded by a higher-temperature envelope. This central region contains newly launched particles with high velocity that have not yet significantly interacted with the ambient medium.

At larger vertical distances, the ejected particles begin to interact with the surrounding gas, converting kinetic energy into thermal energy. The resulting jet head, enclosed by the bow shock, shows elevated temperature and entropy relative to the unshocked ambient medium.

By $t=40\,\mathrm{Myr}$, the overall morphology remains broadly similar to that at $t=20\,\mathrm{Myr}$. However, by $t=60\,\mathrm{Myr}$, the tip of the jet lobe appears sharper on both the density and temperature maps. At $t=80\,\mathrm{Myr}$, this sharpening extends to the bow shock front, leading to deviations from the expected self-similar evolution. At the end of the simulation at $t\sim 100\,\mathrm{Myr}$, the jet extends to a length of 300 kpc, and the shock morphology of the lobe and bow clearly departs from the self-similar analytic profile.

A notable feature appears near the jet origin: while the central region initially retains a thermal state nearly identical to that of the undisturbed ambient gas, it gradually exhibits signs of jet-ambient interaction. By $t=100\,\mathrm{Myr}$, two distinct high-entropy streams emerge near the center. This feature reflects the reservoir setup in the RG launching scheme: particles are arranged in layered grids and launched in randomized order per layer. Inner-layer particles do not start to launch until the outer layers are depleted. Since the reservoir extends to 10 kpc in height, the ambient particles within this region initially remain unperturbed (given a typical smoothing length $h=3.6\,\mathrm{kpc}$ of the ambient gas particles), but are gradually influenced as the inner particles begin to launch.

In terms of quantitative morphology, the jet lobe length and bow shock shape in the FID run generally follow the analytic scaling laws proposed by Kaiser and Best (2007), particularly during the intermediate stages ($t=20$-$60\,\mathrm{Myr}$), where the jet propagation approximately follows the $L_{J}\propto t^{3/5}$ relation. However, at later times ($t\gtrsim 80\,\mathrm{Myr}$), deviations become apparent as the bow shock narrows and the bow shock front steepens beyond the expected aspect ratio. These changes likely reflect non-linear effects such as the build-up of pressure gradients near the jet head, interaction-induced asymmetries, and resolution-dependent turbulence dissipation. The progressive deviation from the self-similar morphology suggests that jet propagation in a stratified medium can enter a kinetically dominated regime once mixing and backflows are suppressed.

The jet exhibits overall symmetry along both the $x$- and $z$-axes, indicating that the RG scheme provides a physically reasonable and stable jet launching strategy.

4.1.2 Gas characteristic in the simulated jets

Figure 4 presents a series of slice plots showcasing key gas properties at the end of the fiducial run ($t=100\,\mathrm{Myr}$). These maps allow us to identify the distinct physical components of the jet structure; namely, the jet material, the lobe, and the bow shock, through their thermodynamic and kinematic signatures.

In the left panel, we show the temperature, density, and entropy distributions. The jet lobe appears as a hot, low-density cavity with markedly elevated entropy compared to that of the ambient medium. This region is the clearest thermodynamic indicator of the jet influence. Notably, there is no sharp boundary between the lobe and the jet material, which together form a continuous structure of entrained and injected gas. Outside the lobe, a broader region of shocked ambient gas is barely visible, enveloped by the bow shock front. The surrounding region exhibits an enhanced density and temperature as a result of the shock gas being compressed and heated by the expanding jet.

The middle panel of Figure 4 displays thermal pressure, ram pressure, and the contrast between the two. These quantities are effective in distinguishing between different jet-driven regimes. The jet material near the axis is characterized by a high ram pressure and moderate thermal pressure, producing a pronounced pressure contrast. This region also exhibits a recollimation pattern, the so-called jet spine, where the ram pressure peaks as the flow tightens along the axis. The jet maintains a nearly uniform, high-velocity stream along the axis until it decelerates abruptly at the jet head. In contrast, the shocked ambient gas surrounding the spine shows ram and thermal pressures of comparable magnitude. In this region, the gas exhibits a bulk velocity that approximates the local sound speed as it moves away from the jet axis, creating a distinct boundary that is differentiated from the surrounding, non-turbulent ambient gas. Although the jet lobe can be traced in this panel, its boundaries are less clearly defined than those in the thermodynamic maps. It appears as an area within the shock front where the thermal pressure and ram pressure are nearly equivalent.

The right panel of Figure 4 presents the vertical velocity ($v_{z}$), the radial velocity ($v_{r}$) and the Mach number. These kinematic diagnostics effectively separate the different dynamical components of the jet system. The ambient medium exhibits negligible velocity, and thus a near-zero Mach number. In contrast, the gas within the shocked boundary moves at approximately sonic speeds. The jet material itself is clearly identifiable as a supersonic stream along the axis, registering the highest values of $v_{z}$ and $v_{r}$. Meanwhile, the jet lobe is particularly well distinguished in the $v_{z}$ and $v_{r}$ maps by its characteristic backflow — regions of negative velocity — indicating material that is being recycled inward from the edges of the lobe. These backflows contribute significantly to the lobe’s pressure support and the redistribution of energy and metals. The Mach number map confirms that the lobe material is near-sonic speed, distinguishing it from the supersonic spine and the sub-sonic ambient gas. The jet material that follows the jet axis is characterized by the highest values of $v_{z}$ and $v_{r}$ in the velocity maps, making it the only supersonic component observable in the Mach number map.

4.1.3 Phase diagrams

Figure 5 presents a set of phase diagrams at $t\sim 100\,\,\mathrm{Myr}$ for the fiducial run, providing thermodynamic and kinematic diagnostics of the gas impacted by the jet. We show the gas distribution in the density-temperature ($\rho-T$) plane (left), thermal pressure–entropy ($P_{\mathrm{th}}-K$) plane (middle), and thermal–ram pressure ($P_{\mathrm{th}}-P_{\mathrm{ram}}$) plane (right). The color scale indicates the mass density of the gas in each phase space bin.

The $\rho-T$ diagram and $P_{\mathrm{th}}-K$ diagram are mathematically related via affine transformations, but we display them separately to highlight different aspects of gas thermodynamics. In both panels, one can identify a strong population aligned with the ambient pressure level ($P\approx P_{0}$), extending over a wide entropy range ($K_{0}<K<10^{3}K_{0}$). This distribution reflects the thermodynamic evolution of the shocked ambient gas and the lobe material. A secondary feature appears along the adiabatic (isentropic) relation, with modest enhancement in entropy beyond $K_{0}$, corresponding to gas undergoing expansion with limited shock heating.

A distinct subpopulation at constant density and elevated temperatures ($T_{0}<T<10^{2}\,T_{0}$) is also visible, associated with jet head where the shocked gas exchanges internal energy with the ambient gas. In contrast, a diffuse ‘delta-shaped’ region in the lower pressure, high entropy quadrant ($P<P_{0}$, $K>K_{0}$) represents the newly launched jet gas that endures a nearly free expansion in the jet spine. Another diffuse component spans the regime $P_{0}<P<10^{2}\,P_{0}$, indicating the presence of partially mixed material at intermediate pressures.

The right panel shows the phase diagram in the thermal-ram pressure plane, which offers a diagnostic of jet-ambient interactions. In the phase diagram, two branches emerge. One passes through the intersection of $P_{\mathrm{th}}=P_{0}$ and $P_{\mathrm{ram}}=P_{\mathrm{launch}}$, whose slope matches the adiabatic process assuming a constant speed. This branch corresponds to recently launched, supersonic jet material, which moves without much deacceleration and follows adiabatic expansion. The second branch follows $P_{\mathrm{th}}\approx P_{0}$ at small $P_{\mathrm{ram}}$ values, but turns to align with $P_{\mathrm{th}}\sim P_{\mathrm{ram}}$ when $P_{\mathrm{ram}}$ goes higher. This population represents the ambient gas: the low $P_{\mathrm{ram}}$ part following $P_{\mathrm{th}}\approx P_{0}$ is the shock heated ambient gas, and the $P_{\mathrm{th}}\sim P_{\mathrm{ram}}$ trend marks the ambient gas in the jet lobe that interacts with the jet particles. The region where two branches meets near $(P_{\mathrm{th}},P_{\mathrm{ram}})\sim(P_{0},10^{-1}P_{\mathrm{launch}})$ marks the transition regime where the shocked jet material mixes with the ambient medium. This structure provides insight into the evolution of pressure support and turbulence inside the bow shock.

4.1.4 Jet lobe decomposition

The findings in sections 4.1.2 and 4.1.3 motivate us to decompose the jet into the approach detailed in Section 3.5. We apply a filter to the snapshot using the criteria for the shocked gas ($\mathcal{M}>0.0316$) to remove ambient unshocked gas, then further classify the gas based on two criteria: gas origin (from jet ejection or ambient) and jet speed (subsonic or supersonic, determined by Mach number). Figures 6 and 7 are the slice plots and the phase diagrams for jet decomposition, respectively.

The jet particles with $\mathcal{M}>1$ (second column in Figure 6; red points in Figure 7) correspond primarily to the particles ejected in the jet spine. These particles are dispersed across the $\rho$-$T$ plane; near the jet axis, a substructure with low entropy, low temperature, high ram pressure and high Mach number is observed. These particles have minimal interaction with the ambient gas, resulting in negligible kinetic energy loss. At larger radii, jet particles exhibit higher entropy and lower ram pressure, indicating interaction with the surroundings and a reduction in speed. A gradient in the Mach number is observed along the temperature axis on the phase diagram, as anticipated: initially, jet particles possess high velocity (high Mach number) and low temperature, then decelerate and heat up upon interaction. From the $P_{\mathrm{th}}$-$K$ phase diagram, it is clear that this population is situated at the isobaric component and the ‘delta region.’ The $P_{\mathrm{th}}$-$P_{\mathrm{ram}}$ plot shows the inclusion of the adiabatic component, as depicted in Figure 5.

Jet particles with $\mathcal{M}<1$ (third column in Figure 6; yellow points in Figure 7) predominantly coincide with the structure of the jet lobe. This group is characterized by high temperature, low density, and high entropy, suggesting the low-density volume evacuated by the jet particles. Notably, the vertical velocity of these particles presents a backflow with speeds of several thousand km/s. The particles also have a notable negative radial velocity. These negative velocities suggest that low-Mach-number jet particles have reached the tip of the jet and are being pushed back because of the gas pressure. In the phase diagram, these particles sit along the same equal-pressure line as the $\mathcal{M}>1$ component. However, there is a marked difference in the Mach number, particularly at the low-entropy end: a location on the $\rho$-$T$ diagram can simultaneously exhibit $\mathcal{M}>10$ and $\mathcal{M}\sim 0.1$.

The component of ambient gas with $\mathcal{M}>1$ (fourth column of Figure 6; green points in Figure 7) interacts with high-speed jet particles or participates in jet propagation. These particles exhibit similar thermal conditions as the jet particles but have lower Mach numbers. Compared with their spatial distribution to the jet core, these particles form the shock ‘head,’ which is characterized by low entropy, high temperature, and high density, in addition to the jet lobe. In the phase diagram, the ‘jet lobe’ aligns with $P=P_{0}$, similar to the jet gas with $\mathrm{M}<1$, while the ‘jet head’ is closer to the original ambient gas $(\rho_{0},T_{0})$ and is located in the low-entropy region of the $P_{\mathrm{th}}$-$K$ diagram.

The most substantial content in the bow shock is the ambient gas with $\mathcal{M}<1$. In Section 4.1.3, we identified three primary components, with the most prominent having an entropy just above $K_{0}$.

Combining the aforementioned analysis, we can trace the evolution of jet particles in phase space. Initiated at $(\rho_{0},T_{0})$, the jet particle has a high velocity, as evidenced by an elevated Mach number. Initially, dependent on the previous cavity scale, the particle experiences near-free expansion until reaching an appropriate smoothing length, putting it in the ‘delta region’. During this process, heat convection with surrounding particles and interaction-driven deceleration occur, resulting in a Mach number gradient (see the top panels of Figure 7). The particle then gradually approaches the isopressure line $P=P_{0}$ due to the interaction with jet lobe material. As the high-speed particle collides with the hot jet head, its speed drops. Meanwhile, density and temperature drastically increase as a result of shock compression, converting kinetic energy into thermal energy. This process shifts the particle toward the upper right phase diagram into another diffuse region. In real space, particles rebound from the tip, residing within the jet lobe. Here, the particle decelerates to subsonic speeds, losing internal energy to the ambient gas. The final state is achieved when the local pressure equilibrium is met, aligned with the $P=P_{0}$ branch and $\mathcal{M}<1$.

In the phase diagram of ambient gas, the primary shaping factors are shock compression and jet particle interactions. Unlike the jet particle, the shock front moves through the ambient gas, causing an increase in density, temperature, and entropy, resulting in an iso-entropy component at $K$ slightly above $K_{0}$. The ambient gas’s iso-density component arises from heat exchange at the jet head, where super-heated jet particles cool after shock compression. Ambient gas particles enclosed by the shock front reach pressure equilibrium with their surroundings, forming a component along the iso-pressure line.

4.1.5 Jet size evolution

Figure 8 presents the time evolution of the characteristic spatial scales in our fiducial jet simulation, including the lengths and radii of the lobe and bow shocks. These are compared against the analytic self-similar model of Kaiser and Best (2007) and the simulation results of H23. We also included the measurement of MESH run as green lines in the figure, which is discussed later in Section 4.3.

We define the extent of the lobe using an entropy-based criterion (§3.5.1), which reliably tracks the thermalized, low-density cavity produced by jet activity. Overall, our measured lobe length (solid red line) follows the theoretical scaling $L_{J}\propto t^{3/5}$ reasonably well across most of the simulation, confirming that the the expansion of bow shock broadly adhere to self-similar expansion. However, at late times ($t\gtrsim 70$ Myr), a mild deviation appears: the lobe length curve flattens slightly, diverging by $\sim$10% from the analytic curve. This may reflect residual turbulence, finite-resolution effects, or asymmetry in the lobe structure near the termination phase.

In contrast, the simulation by H23 (blue solid line) displays slightly longer jet lobes at late times, exceeding both our results and the theoretical expectation. This may be due to differing numerical viscosity, launching schemes, or entropy thresholds used to define the jet-ambient interface (cf. §4.2). Notably, our lobe structure remains more compact, suggesting enhanced mixing or greater dissipation of kinetic energy into heat (see Fig. 9).

The evolution of the lobe radius ($R_{\ell}$, dot-dashed red line) initially lags behind the analytic prediction and the H23 result for the first $\sim$20 Myr. During this early phase, the radial width closely tracks the bow shock radius ($R_{s}$), hinting at resolution limitations that suppress early lateral expansion. After $t\sim 20$ Myr, however, $R_{\ell}$ begins to converge toward the expected scaling, albeit consistently remaining slightly narrower than the theoretical or H23 profiles.

As with lobe length, delayed growth in lobe width also suggests a departure from self-similar evolution. An increase in the lobe aspect ratio (length-to-radius) was predicted analytically by Alexander (2002) to arise from gravity-driven perturbations. In our view, however, this mechanism cannot explain the analogous behavior observed in our simulations, since gravity is not included. The grid-based simulation of Gaibler et al. (2009) found that the aspect ratio increases for a light jet propagating in a uniform ambient medium. A direct comparison between our results and those of Gaibler et al. (2009) is not straightforward, since our setup involves a heavy jet. However, we observe that the light jet exhibits a steadily increasing aspect ratio, whereas the heavier jets show a rapid initial rise followed by much slower growth, consistent with the evolution trend seen in our simulation around $t\simeq 20\,\mathrm{Myr}$.

The bow shock structure, traced via entropy and pressure discontinuities (§3.5.2), reveals a coherent and egg-shaped expansion of the shock front. The bow shock length ($L_{c}$, dashed red line) remains modestly ahead of the jet lobe throughout, indicative of a persistent contact discontinuity between the forward shock and the jet head. The radial expansion of the shock front ($R_{s}$, dotted red line) maintains a nearly constant aspect ratio relative to $L_{s}$, suggesting that the shock front as a whole evolves self-similarly, even when the internal lobe structure deviates due to mixing, turbulence, or dissipation.

Overall, these results validate the applicability of analytic jet models for describing large-scale bubble growth but also highlight the need for high-resolution simulations and entropy-based diagnostics to capture deviations from self-similarity at late times.

4.1.6 Energy composition evolution

Figure 9 shows the time evolution of the energy partition in the fiducial run, decomposed into kinetic (dot-dashed line), thermal (dashed line), and their measured sum (thick dotted line). For reference, the total injected energy from the jet source is shown as a solid line. The red color is used to represent our results, while the blue color indicates the rescaled 222The total energy budget in H23 is $\sim 0.05$ dex higher than the prediction of $P_{\mathrm{jet}}=10^{46}\,\mathrm{erg/s}$, so we rescaled the energy budget measured in H23 to eliminate the difference. Rescaling has a minor effect on our conclusions. results from H23. Similar to Figure 8, we also indicate the energy budget in MESH run, which will be discussed in Section 4.3.

We find that total energy is conserved to within 1% throughout the simulation, as confirmed by the near-perfect agreement between the measured (thick dotted) and injected (solid) curves. The bottom panel quantifies this small discrepancy, with $\Delta E_{\rm gas}/E_{\rm jet}$ fluctuating around unity by less than 1%. This level of consistency confirms that our numerical setup captures shock heating and momentum deposition with high fidelity, with negligible artificial dissipation or numerical loss.

A key distinction from the findings of H23 emerges in the time-dependent breakdown between thermal and kinetic energy. In their simulation (blue lines), the thermal-to-kinetic ratio remains nearly constant throughout the full runtime, indicating a steady energy-conversion efficiency during the jet’s evolution. Notably, the kinetic energy measured by H23 remains about two times higher than the thermal energy in the whole simulation. In contrast, our run shows a clear temporal trend: kinetic energy dominates the budget during the initial $\sim$20 Myr, but the thermal component overtakes thereafter and continues to grow more rapidly.

The energy budget was also discussed in previous grid-based jet simulations. In Gaibler et al. (2009), the fraction of thermal energy was found nearly unchanged throughout the simulation, hinting the self-similarity of jet evolution; however, even for the heaviest jet in their simulation, more than half of the energy ejected is in the form of thermal energy, which is in agreement with our kinetic-thermal energy partition at late stage of the jet evolution. Hardcastle and Krause (2013) discussed the energy budget in their Figure 9 and 10. They reported that the contribution of thermal energy to the overall energy budget exceeds that of kinetic energy. The gravitational potential energy remains around $10\%$ of the total energy, implying the role of gravity stays secondary in this stage of jet evolution.

4.2 Jet morphology across parameter variations

We compare the fiducial simulation (FID) with a series of controlled variants to assess how various numerical and physical parameters influence the morphology, energetics, and evolution of AGN jets. The tested variations include: (i) the jet particle launching mechanism, (ii) artificial viscosity scheme, (iii) jet active duration (duty cycle), and (iv) mass resolution. Simulation setups are summarized in Table 1, and the resulting differences are illustrated through slice plots (Fig. 10, 11, 12, 13) and the time evolution of lobe lengths (Fig.14).

Launching mechanism

Figure 10 compares three particle launching schemes: the grid-based reservoir model (FID), randomized launch from the same reservoir (RR), and stochastic spawning at runtime from a fixed radius (SR). While all produce broadly consistent macroscopic morphologies with a central jet spine, expanding lobes, and surrounding gas, subtle yet significant differences arise in their forward propagation.

The FID and RR runs yield similar lobe lengths ($\sim$300 kpc), indicating robustness to random particle ordering as long as the launching geometry is preserved. The SR run, however, lags behind at $\sim$250 kpc. This 15-20% deficit arises from its injection method: by placing new particles at $r=5$ kpc rather than 10 kpc (as in FID and RR), the SR jet delays the delivery of momentum to the bow shock. Using $v_{j}\sim 10^{4}$ km/s, the time delay is roughly $0.3$ Myr, which translates to a shorter jet during the linear growth phase that persists into the self-similar regime.

Entropy maps further reveal that the SR lobes are slightly more isotropic and rounded at the tip, consistent with reduced axial thrust and broader angular energy dispersion. These differences highlight that even when the total energy and momentum is conserved, the specific details of how the particles are initialized can affect the collimation, directionality, and eventual observables, particularly in subgrid models of feedback.

Artificial viscosity and jet coherence

The contrast among FID, StV, HVisc, and LLVisc, as shown in Fig. 11, demonstrates the critical influence of numerical viscosity. FID uses the time-dependent Cullen–Dehnen switch to suppress viscosity in shear flows and only activates it near shocks ($\alpha_{\rm AV}$ ranging from 0.01-2). In contrast, StV adopts a constant, maximally dissipative $\alpha_{\rm AV}=2$, producing smoother but overly laminar flows. As a result, the length of the StV jet lobe reaches $\sim$40-50% greater than the FID case by $t=100$ Myr. We interpret this finding to suggest a high thermal pressure behind the jet head in the StV setup. Excess thermal pressure arises from extra thermal energy production due to constant viscosity in solenoidal flow, unlike TDV, where artificial viscosity is applied solely during fluid compression. At the same time, the high viscosity may effectively clear the ambient gas and carve the low-density region along the jet spine, thereby reducing kinetic energy loss along the jet spine. Eventually, the higher viscosity allows more energy to be directly channeled into headward propagation, resulting in narrower and more collimated lobes.

A similar rationale can be applied to the HVisc and LLVisc runs. In the HVisc run, the increased artificial viscosity behaves similarly to the StV run on emptying the ambient gas along the jet spine, leading to a deviation from self-similar growth around $t\sim 20\,\mathrm{Myr}$, similar to the behavior seen in the StV run. Throughout the simulation, the viscosity in the TDV model might decrease when no shock is detected, whereas it is constant in the standard viscosity model. This difference results in the sublinear jet lobe expansion in HVisc, unlike the steady linear growth of the jet lobe in the StV run. In contrast, in the LLVisc run, the initial lower resistance in the ambient medium allows the jet particle to achieve the longest lobe among the viscosity models at the early stage of the jet evolution. However, by $t\sim 20\,\mathrm{Myr}$, the low viscosity and overmixing between the jet and the ambient hinder the jet head’s ability to sweep material in the desired direction, ultimately resulting in the shortest jet lobe among the AV variants at the end of the simulation.

Among the grid-based jet simulations, the mere example that contains artificial viscosity is by Krause and Camenzind (2001). Interestingly, although the artificial viscosity suppressed the KH instability in their simulations, it also widened the jet, seemingly the opposite of our results. Despite differences in our setups and theirs, artificial viscosity drives bulk motion across the velocity discontinuity in both situations.

Jet duty cycle

Figure 12 presents the time evolution of jets with different injection durations: full-duration FID ($t_{\rm jet}=100$ Myr) and truncated cases T50 ($t_{\rm jet}=50$ Myr) and T20 ($t_{\rm jet}=20$ Myr). Despite their early termination, key features, such as jet spines and lobes, persist for several tens of Myr. At $t=21$ Myr, the T20 run still shows a collimated spine; likewise, T50 remains coherent at $t=51$ Myr. However, these features dissipate about 30 Myr after shut-off.

In T20 and T50, lobe growth stalls once momentum injection ceases, but shock front expansion continues. Vertical-velocity and Mach-number slices confirm that the forward flow decays, while backflows and side expansion dominate. Entropy and pressure maps show that cavities remain morphologically distinct even when the axial thrust is gone, consistent with ’ghost bubble’ interpretations in cluster observations (e.g., Churazov et al., 2001; Reynolds et al., 2005).

Sudden decreases in lobe length seen in Figure 14 (at $t\sim 65$ Myr for T20 and $\sim 80$ Myr for T50) result from our entropy-based lobe identification. When particles near the jet front cool or mix into the ambient medium, their entropy drops below threshold, causing them to fall out of the lobe definition and trigger apparent contraction. This illustrates how tracer definitions can significantly impact quantitative morphology metrics during decaying phases.

Resolution dependence.

Figures 13 and 14 compare simulations spanning two orders of magnitude in mass resolution ($1.8\times 10^{5}$ to $1.15\times 10^{7}\,M_{\odot}$). Despite this large range, the global lobe morphology and final extent are relatively robust, with all runs reaching $\sim$300 kpc by $t=98$ Myr. The lobe length tracks the analytic scaling $L_{J}\propto t^{3/5}$ within $\sim$10%, except for early-time deviations in LoD and late-time acceleration in Hi.

Internal structures, however, degrade significantly at lower resolution. The coherent spine and sharp bow shock seen in Hi and FID become poorly defined in LoC and LoD, where numerical viscosity dominates and mixing increases. The vertical velocity and Mach number drop substantially, reflecting reduced thrust efficiency. Entropy maps show broader, rounder lobes, with less directional elongation and greater transverse dispersion.

In grid-based jet simulations, Rosen et al. (1999) reported an increase in mass entrainment as the spatial resolution was raised. Similarly, Krause and Camenzind (2001) showed that using a finer grid permits more complex instabilities to form at the jet head, leading to a broader jet and a reduction in jet speed, which contrasts with the behavior seen in our simulation. In our case, increasing the resolution in SPH simulations enhances the velocity contrasts, which in turn strengthens the effect of artificial viscosity in these high-resolution SPH runs, contrary to those grid-based simulations.

Importantly, FID shows weaker sensitivity to resolution than StV due to its AV switch: since AV activates only in converging flows, excess diffusion is avoided in marginally resolved regions. In contrast, the StV scheme applies uniform damping, leading to underresolved feedback and exaggerated lobe variations at low resolution.

These findings suggest that while feedback effects can be studied at typical cosmological resolutions, finer details such as metal entrainment, turbulence spectra, or jet stratification require high-resolution runs and controlled AV schemes.

4.3 Grid-based Counterpart

In this section, we present the results from MESH run and compare them with our SPH-based simulations.

The jet morphology is presented in Figure 15. Because of the higher spatial resolution of MESH run, we can resolve the structures in much finer detail. The jet material, defined as the region with the Mach number $\mathcal{M}>2$, extends to $z\simeq\pm 500\,\mathrm{kpc}$. Its overall transverse size is smaller than in the SPH run, but the sharpened jet head follows the same general behavior as in the SPH jet. Within the jet material, there is a strongly supersonic segment between the launch point and $z\simeq\pm 200\,\mathrm{kpc}$; beyond this distance, the flow remains supersonic albeit much milder. A similar trend is observed in the SPH simulation, but it becomes evident only at higher mass resolution. A lobe is still present outside the jet material, but its morphology and dynamics differ from those in the SPH case. The lobe is generally narrower, and vortex-like features develop in the middle of the lobe, in contrast to the more inflated and puffy lobe seen in the SPH runs. The bow-shock morphology broadly agrees with that of the SPH simulations; the bow-shock front is sharply pointed, indicating a late-time acceleration of the jet head. As a result, the overall morphology more closely resembles that of the StV or HVisc runs with enhanced artificial viscosity than that of the fiducial SPH setup.

The evolution of the jet size in MESH run is shown in Figure 8. For both the jet lobe and the bow shock lengths, MESH run follows the self-similar model and matches the fiducial SPH run to $t=50\,\,\mathrm{Myr}$; beyond this time, the jet advances more rapidly than predicted by the self-similar solution, consistent with the sharp bow shock visible in Figure 15. This late-time acceleration also appears in the fiducial run, but is substantially weaker; only the StV or HVisc runs exhibit a comparable level of acceleration (Figure 14). The measured lobe radius exhibits larger fluctuations relative to the self-similar expectation, which can be attributed to the higher spatial resolution of the mesh grid and the limitations of using 2D slice data for measurements.

In Figure 9, we also present the energy budget of MESH run. We find that the evolution of the thermal-kinetic energy partition in MESH run is generally consistent with the fiducial SPH run: the kinetic energy initially dominates but is eventually overtaken by the thermal component. However, in MESH run this transition occurs later at about $t=50\,\,\mathrm{Myr}$, compared to $t=30\,\,\mathrm{Myr}$ in the FID run. In addition, the thermal-to-kinetic energy ratio in MESH run remains lower than in the FID run at the end of the simulation.

5 Discussion

In this study, we conducted simulations of AGN jets with an idealized ICM setup. The primary focus of our simulations is the convergence of hydrodynamics when varying conditions of jet launching implementations, artificial viscosity configurations, mass resolutions, and jet launching durations. To accomplish this, we performed 13 simulation variants in addition to the fiducial one (Table 1). In Section 5.1, we compare our simulation results with self-similar predictions. Section 5.2 is a comparative analysis of our simulation with H23, wherein we sought to closely replicate their simulation setup. In Section 5.3, we compared our SPH simulation with the grid-based counterpart. Finally, Section 5.4 explores the connection between our jet-time simulation variants and the Milky Way Fermi bubble.

5.1 Comparison with analytical solution

In this section, we compare our simulation results with analytical predictions. Figure 8 shows the evolution of the jet lobe in our fiducial simulation alongside the self-similar model of Kaiser and Best (2007). Overall, the simulation exhibits good agreement with the self-similar theory, although some deviations are present in the detailed growth trend of the lobe.

As shown in Figure 14, the jet lobe initially grows in a sub-self-similar manner, following $L\propto t^{\alpha}$ with $\alpha<3/5$ for $0<t<40\,\,\mathrm{Myr}$. Beyond this point, the growth rate increases, and the lobe length eventually catches up with the self-similar prediction by the end of the simulation.

When examining the jet lobe evolution across different artificial viscosity models, we find that all simulations—except for the LLVisc run, which maintains low viscosity throughout—deviate from self-similar growth around $t\sim 50\,\,\mathrm{Myr}$.

Combining this fact with the evolution of the jet lobe of other variants of artificial viscosity, we find all these simulations, except LLVisc, which has a low viscosity throughout the simulation, deviate from self-similar jet lobe growth around $t\sim 50\,\mathrm{Myr}$.

An alternative approach to modeling jet head motion, proposed by Begelman and Cioffi (1989), balances the jet’s thrust against the ram pressure of the ambient medium. In contrast to the self-similar framework, this model explicitly incorporates the effect of ambient pressure, making it particularly applicable to the later stages of jet evolution. Specifically, it becomes relevant once the jet length surpasses the characteristic scale $L_{2}$ (Equation 3), beyond which the ambient pressure is no longer negligible. This transition marks the onset of a pressure-confined regime, where the expansion decelerates and the lobe morphology becomes increasingly governed by pressure equilibrium with the surrounding environment.

Assuming that the jet eventually attains a terminal velocity while propagating through a homogeneous ambient medium, one can estimate the lobe’s advance speed by equating the jet thrust with the ram pressure of the surrounding gas (see Appendix A for details). This balance yields the following expression for the jet head velocity:

[TABLE]

where $\rho_{\mathrm{A}}$ is the ambient gas density, $S_{\mathrm{w}}$ is the working surface area at the jet head, $v_{\mathrm{j}}$ is the velocity of the jet material, and $P_{\mathrm{j}}$ is the jet thrust (momentum flux). This expression captures the transition between the low-thrust and high-thrust regimes and explicitly incorporates the geometrical effect of the jet-ambient interface through $S_{\mathrm{w}}$.

Unlike the self-similar expression in Equation 7, the lobe propagation speed derived in Equation 12 is not explicitly time-dependent. Nevertheless, it may still exhibit an implicit time dependence in practice, arising from the gradual evolution of the working surface area, spatial variations in the ambient density, or other environmental inhomogeneities that affect the jet-ambient interaction over time.

As the working surface area $S_{\mathrm{w}}$ is difficult to measure in our simulation, we try to estimate it from approximate relations with other well-defined quantities. One may approximate $S_{\mathrm{w}}$ with the cross-sectional area of jet head $\pi R_{J,\mathrm{head}}^{2}$. To measure the radius, we first identify the jet material as jet-launched particles with $\mathcal{M}>2$ and measure its length $L_{J}$. We then define the jet head particles as the jet material located within the outermost $20\%$ of this length.

[TABLE]

Finally, the jet head radius $R_{J,\mathrm{head}}^{2}$ is measured with the method in Section 3.5. This allows us to evaluate the jet head velocity using Equation 12 with a physically motivated estimation.

Figure 16 compares the jet-head velocity measured in the FID simulation with the self-similar (Equation 5) and pressure-equilibrium analytic (Equation 12) predictions. During the first $\sim 20\,\,\mathrm{Myr}$ the simulation tracks the self-similar solution closely.

The jet head velocity drops below the self-similar value at about $t\simeq 20\,\,\mathrm{Myr}$, but it speeds up again around $t\simeq 50\,\,\mathrm{Myr}$. This behavior suggests that the working surface is becoming more focused, consistent with the pressure equilibrium expected from the jet head radius.

The sharpening of jet head is presumably due to the reconfinement geometry of the jet material (e.g. Komissarov and Falle 1998), which is governed by the properties of jet and the external pressure: the jet flow first recollimates and then reconfines, producing a thin, spindle-like morphology along the jet axis. In Figures 3 and in the animated version of Figure 4, one can clearly see that the jet head exhibits a flat cross section in the middle of the simulation, but becomes reconfined by the end. This interpretation is also consistent with the variation in jet lobe lengths discussed in Section 4.2, where the choice of viscosity configuration modifies the degree of jet refinement and, in turn, influences the advance speed of the jet head.

5.2 Comparison with SPH-based Simulations

In this study, we retained most choices for jet simulation parameters in our fiducial run to remain consistent with the fiducial setup of H23. Our fiducial simulation reproduces many of the features reported by H23, although several notable differences emerge.

Most prominently, in terms of jet morphology, we do not observe the pronounced oscillatory re-collimation pattern characteristic of their fiducial self-similar run. Instead, re-collimation in our simulation occurs only once, producing a pencil-shaped jet spine. Meanwhile, we notice that the shape of jet spine in our SPH simulation is in good agreement with other well-established grid-based simulation. We argue that the distinct re-collimation patterns in fact arise from differing approaches to jet launching: H23 appears to inject a swarm of particles simultaneously, producing a recurring sequence of bubbles along the jet axis. Such a launching scheme may break the steady-state condition and axisymmetry assumptions commonly adopted in the literature (e.g., Bodo and Tavecchio 2018; Gourgouliatos and Komissarov 2018). In our simulation, we adopted a more randomised and smooth launching scheme, which helps the jet remains steady and reaches a good agreement with the grid-based simulation.

Despite this difference in the jet spine, the morphology of the jet lobe and bow shock is broadly consistent with H23. In both cases, the lobe begins as a quasi-spherical, overpressured bubble that gradually stretches into a slender, cigar-shaped cavity with a sharp head, while its radial growth remains relatively modest (Fig. 8). This evolution reflects a distinct morphological transition, as discussed in Section 5.1. The bow shock outside, characterized by only mild enhancements in temperature and density, preserves an approximately self-similar aspect ratio throughout the evolution in both our simulations and those of H23.

Quantitatively, we do find a number of differences in the jet evolution measurements (Figures 8 and 9) as follows:

Jet structure sizes

The evolution of jet structure sizes of our fiducial simulation is compared to H23 in Figure 8. The difference is mainly due to the measurement at an early time of the simulation. The length of the jet lobe measured in H23 is initially shorter, gradually catching up with the prediction of theory and our measurement around $t<20\,\mathrm{Myr}$. The radius of the jet lobe of our measurement deviates from the theoretical prediction and H23 during $t<20\,\mathrm{Myr}$, which is likely a resolution-related issue, given $\sim 10\,\mathrm{kpc}$ the scale of the jet launching reservoir.

At the same time, we also recall the difference in our lobe identification criterion (Section 3.5.1) from H23. H23 classified lobe material solely by temperature, labeling any gas particle that ever exceeded $5\times 10^{8}\,\mathrm{K}$ (their fiducial threshold) as part of the lobe. A temperature-only criterion blurs the distinction between genuine lobe gas and shocked ambient gas because both phases can reach temperatures much higher than the ambient value. In contrast, our entropy-based selection effectively separates the two components by incorporating density information. Jet lobe gas is simultaneously hot and underdense, whereas shocked ambient gas is hot but overdense. We argue that this combined density–temperature signature provides a cleaner and more physically motivated identification of the lobe, though its performance is limited when the resolution effects play an important role.

Energy partition

Because our idealised setup disables gravity, cooling, magnetic fields, and cosmic rays, the only energy source is the kinetic power injected by the jet, which is subsequently converted to thermal energy via shocks and shear. In a perfectly self-similar system, the time derivatives of kinetic and thermal energy should scale in fixed proportion, so that their individual changes ($\Delta E_{\rm kin}$ and $\Delta E_{\rm th}$) remain in a constant ratio, as reported in H23. Figure 9 shows that this is not the case in our run: kinetic energy dominates up to $t\simeq 20\;\mathrm{Myr}$, after which thermal energy grows more rapidly and overtakes the kinetic component. The two curves cross and the thermal energy dominates at the late stage of the simulation, whereas the ratio of thermal and kinetic energy remains roughly 1:2 in the SWIFT simulations of H23. The difference may arise from the choice of SPH formulation, where we adopted the pressure-entropy one while for H23 it was the energy-density scheme.

At late times, the kinetic energy curve in our simulation gradually flattens toward late times, suggesting a saturation in momentum injection as the jet transitions from a momentum-driven to pressure-driven phase. The bow shock expands more isotropically, with buoyant uplift dominating over bulk axial thrust. This late-stage evolution again underscores the importance of feedback modeling that allows for mode transitions, e.g., from ballistic jets to thermally buoyant bubbles, which can influence the long-term energy coupling to the CGM.

Overall, although the large-scale morphology agrees with the analytic self-similar theory, the internal energy partition does not follow the constant-ratio behavior found by H23. This highlights the sensitivity of kinetic-to-thermal conversion – and thus observable X-ray cavity energetics – to seemingly minor implementation details.

5.3 Comparison with Grid-based Simulation

In this study, we investigate the evolution of jets in an idealized environment using the SPH simulation code GADGET4-Osaka, originally developed as a general-purpose tool for cosmological simulations. This approach contrasts with many previous works that ran jet simulations using a dedicated mesh-based hydrodynamic simulation code, motivating us to perform a mesh simulation with identical setup. In this section, we discuss the comparison between SPH and grid-based simulations.

Following the results presented in Section 4.1 and 4.3, we conclude that our SPH code agrees well with the mesh-based code on jet propagation and jet energy budget, albeit the lobe structures are not well-resolved in SPH runs. We find that the jet lobe in our SPH run does not exhibit resolved gas motions, such as vortices, that are observed in grid-based simulations. This is likely a consequence of the low-density nature of the jet lobe: because SPH is fundamentally mass-weighted, there are too few particles in this region to accurately capture its internal kinematics. Despite the limited resolution of the lobe’s internal structure in our SPH calculations, the method nonetheless reproduces the jet material properties and length evolution, as well as the morphology of the bow shock, provided that the artificial viscosity parameters are tuned appropriately. In addition, the SPH framework can readily incorporate a wide range of additional physics—including gravity, radiative cooling, and sub-grid physics such as star formation and cosmic rays, making it particularly well suited for jet simulations in more realistic setups and better aligned with the requirements of galactic astrophysics studies.

A complementary perspective is offered by the clumpy-ISM study of Dutta et al. (2024). Using analytic arguments calibrated with 3-D grid-based simulations, they show that once the effective density ratio between clouds and diffuse gas exceeds $\lambda\approx 20$, the jet head stalls, the bow shock inflates quasi-isotropically, and most of the mechanical power dissipates locally (their Eq. 27 and Figure 3). Our idealised runs occupy the opposite extreme, $\lambda\simeq 1$, where the jet remains momentum-driven and follows the self-similar $L_{J}\propto t^{3/5}$ growth until late-time turbulence steepens the head advance. The two studies therefore bracket the range of ambient inhomogeneity likely encountered by real AGN. Incorporating a clumpy CGM into GADGET4-Osaka — while retaining the low-viscosity, reservoir-based launch scheme validated here — will allow a direct SPH test of Dutta et al. (2024)’s dissipation criterion and of how numerical viscosity modulates the anisotropic vs. isotropic transition.

5.4 Jet remnant evolution

Our truncated-jet runs (T20, T50) demonstrate that an over-pressurised cavity remains visible for tens of Myr after energy injection ceases, with the bow shock continuing to expand buoyantly. Rescaled to Sgr A*, a $P_{\rm j}\sim 10^{41}$ erg s*-1* jet active for $20\!-\!30$ Myr would leave a $10^{56}$ erg thermal bubble comparable to the eROSITA shell, while a younger, more energetic burst could account for the sharper, higher-luminosity Fermi edges. We conclude the aspect ratio of near-spherical eROSITA bubble morphology can be achieved with a shorter but more energetic burst. At the same time, the variation on the jet properties (e.g., lighter collimated jet in Zanni et al. 2003; higher external Mach number and larger opening angle in Krause et al. 2012) or the environment setup (e.g., a steeper density profile or clumpy halo in Dutta et al. 2024) may also contribute to the morphology of the bubble. Incorporating complexities—together with cosmic-ray transport—will be the next step in linking idealised jet simulations to multi-wavelength observations of Galactic outflows.

6 Conclusions

We have performed a suite of hydrodynamic simulations of kinetic AGN jets using GADGET4-Osaka SPH code over a duration of $\sim$100 Myrs, following the work by Huško and Lacey (2023). Our fiducial simulation not only successfully reproduces the principal characteristics predicted by classical analytical models, but also exhibits good agreement with grid-based simulations at limited mass resolution, while retaining sufficient flexibility to incorporate additional physical processes. In a uniform ambient medium, the jet lobe length follows the expected self-similar scaling $L_{\rm j}\propto t^{3/5}$ to within 10% for most of the evolution. Deviations appear at late times when ambient pressure becomes significant, at which point the jet-head velocity transitions smoothly from the self-similar regime to a pressure-equilibrium regime. This transition is well captured by the momentum-balance model derived in Appendix A. To further characterize the jet structure, we decomposed it into distinct components and examined their thermodynamic evolution in the $\rho$–$T$ phase diagram.

We have also explored the sensitivity of jet morphology to various numerical schemes. Switching from a pre-allocated particle reservoir to on-the-fly spawning at a finite distance shortens the final lobe length by 15–20%, due to delayed and more dispersed momentum injection. Nonetheless, the large-scale morphology remains broadly consistent. Artificial viscosity prescriptions play a more significant role: adopting a constant, maximally dissipative viscosity or a high-ceiling time-dependent model enhances jet collimation and extends the lobe by 40–50% relative to the fiducial model. In contrast, a low-viscosity ceiling leads to early mixing and deviation from self-similar growth as early as $t\sim 50\,\,\mathrm{Myr}$.

When jet power is switched off at $t=20$–$50\,\,\mathrm{Myr}$, the resulting fossil lobes decelerate and contract, as high-entropy material gradually mixes into the ambient gas. These results highlight how tracer definitions and feedback history affect late-time lobe morphology. We also tested convergence with mass resolution: varying the particle mass by two orders of magnitude yields final lobe lengths within 10% of the analytic expectation. However, fine structural features—such as bow-shock sharpness, backflows, and metal entrainment—are significantly degraded at lower resolution or with constant-$\alpha$ viscosity, demonstrating the importance of high resolution and adaptive AV schemes for capturing CGM/ICM-relevant physics.

Looking forward, our study lays the foundation for several promising directions. Embedding jets in stratified cluster atmospheres with gravity, cooling, magnetic fields, and cosmic rays will be essential for capturing realistic buoyancy, mixing, and non-thermal pressure support. Coupling this feedback model to self-regulated black hole growth in cosmological simulations will enable studies of jet duty cycles and their broader impact on galaxy populations. Synthetic X-ray and radio maps informed by our detailed lobe structure could provide valuable comparisons to upcoming observations from Athena, SKA, and other facilities. Finally, the implementation of Lagrangian tracers, subgrid turbulence models, or differentiable GPU-accelerated hydrodynamics will help reduce numerical uncertainties and enable field-level inference of jet parameters.

In summary, our results demonstrate that GADGET4-Osaka is capable of reproducing analytic jet scalings, identifying key numerical dependencies, and serving as a robust platform for next-generation studies of AGN feedback in galaxy formation and evolution.

Acknowledgments

We thank the anonymous reviewer for the precious suggestions on the improvement of this manuscript. We are grateful to Filip Huško, Cedric Lacey for the discussion about their work on jet simulation. C.D. thanks Enrico Garaldi for a useful discussion on the artificial viscosity schemes. KN is grateful to Volker Springel for providing the original version of GADGET-4, on which the GADGET4-Osaka code is based. Some of the numerical computations for this study were carried out on the computational cluster idark at Kavli IPMU, the Cray XC50 at the Center for Computational Astrophysics, National Astronomical Observatory of Japan, the SQUID at the Cybermedia Center, Osaka University as part of the HPCI system Research Project (hp240141, hp250119), and the HOKUSAI BigWaterfall2(HBW2) at RIKEN. This research was supported by FoPM, WINGS Program, the University of Tokyo. This work is supported in part by the MEXT/JSPS KAKENHI grant number 20H00180, 22K21349, 24H00002, and 24H00241 (K.N.) and 25H00675 (A.M.). K.N. acknowledges the support from the Kavli IPMU, World Premier Research Center Initiative (WPI), UTIAS, the University of Tokyo. R.C. acknowledges in part financial support from the start-up funding of Zhejiang University and Zhejiang provincial top level research support program.

Data Availability

The snapshots of GADGET4-Osaka jet simulation, as well as the analysed data products, will be shared upon request to the authors. The animated visualisation of the simulation is available on YouTube: https://www.youtube.com/@KNastro-cosmos.

Appendix A Lobe head motion

We model the advance of the jet (lobe) head into the ambient medium by balancing the jet thrust against the ambient ram (and, when needed, thermal) pressure at the working surface (e.g., Bourne and Sijacki, 2017; Begelman and Cioffi, 1989). Considering one side of the jet, Newton’s law across the head gives

[TABLE]

where $\dot{m}_{j}$ is the jet mass outflow rate, and $\bar{S}_{w}$ is the effective working surface area through which the head couples to the ambient gas. We parametrize

[TABLE]

with $S_{w}$ the geometric area and $\xi$ a dimensionless efficiency factor that accounts for obliquity/porosity of the working surface.

Note that in this analysis, only one side of the jet (or lobe) is considered, i.e. $\dot{m}_{\mathrm{j}}=\dot{M}_{\mathrm{j}}/2$. $\vec{v}_{\mathrm{H}}$ is the velocity of the lobe head, $\vec{v}_{\mathrm{A}}$ and $\rho_{\mathrm{A}}$ are the velocity and the density of the ambient gas respectively. Upon being ejected from the black hole, jets undergo minimal deceleration while travelling through the low-density environment within the jet lobe cavity. However, once they encounter the ambient medium, the jet rapidly decelerates, reaching a terminal velocity.

Assuming a static ambient medium ($\boldsymbol{v}_{\mathrm{A}}=\boldsymbol{0}$), a fixed jet direction, and a quasi-steady head ($\boldsymbol{a}_{\mathrm{H}}=\boldsymbol{0}$), the axial momentum balance reduces to

[TABLE]

Define $v_{0}\equiv\dot{m}_{j}/(\rho_{\mathrm{A}}\bar{S}_{w})$. When $v_{j}\gg v_{\mathrm{H}}$ (the classical limit; Bourne and Sijacki 2017; Begelman and Cioffi 1989), Eq. (A3) gives

[TABLE]

Without taking $v_{j}\!\gg\!v_{\mathrm{H}}$, the exact solution of Eq. (A3) is the positive root of the quadratic,

[TABLE]

which satisfies $0<v_{\mathrm{H}}<v_{j}$.

To include thermal pressure, work in the head frame and balance the normal stresses (e.g., Matzner 2003),

[TABLE]

Relating the jet density at the working surface to the mass flux,

[TABLE]

and neglecting thermal pressures ($p_{j},p_{A}\to 0$) recovers the standard thrust-ram result (Marti et al., 1994),

[TABLE]

When thermal pressure differences matter, define

[TABLE]

Combining Eq. (A6) with these definitions yields a quadratic for $x\equiv v_{H}/v_{j}$ whose physical root ($0<x<1$ for light jets) is

[TABLE]

which reduces exactly to Eq. (A8) for $\beta\to 0$.

If one prefers to normalize the pressure jump by the jet ram pressure, let

[TABLE]

and write Eq. (A9) equivalently as

[TABLE]

Useful limits.

(i) $\beta\to 0$: $v_{\mathrm{H}}/v_{j}\to\sqrt{\eta_{\rho}}/(1+\sqrt{\eta_{\rho}})$ (Eq. A8). (ii) $\eta_{\rho}\to 0$ at fixed $\beta$: $v_{\mathrm{H}}/v_{j}\to\sqrt{\beta}$, i.e. the head speed is set by jet overpressure against ambient inertia. (iii) Large overpressure $\Pi\gg 1$: $v_{\mathrm{H}}/v_{j}\simeq\sqrt{\eta_{\rho}\,\Pi}$ until $v_{\mathrm{H}}$ approaches $v_{j}$.