AGN jet feedback on a moving mesh: gentle cluster heating by weak shocks and lobe disruption

ABSTRACT

1 INTRODUCTION

Galaxy clusters, which due to the nature of hierarchical structure formation, build up through successive mergers of smaller systems over cosmic time to become the most massive, gravitationally bound systems in the Universe. Comprising a massive dark matter halo, menagerie of galaxies, and the hot intracluster medium (ICM), they represent a rich ecosystem in which to study a range of astrophysical phenomena, as well as being complimentary cosmological probes (see e.g. Allen, Evrard & Mantz 2011; Kravtsov & Borgani 2012; McNamara & Nulsen 2012, for reviews). The ICM, which radiates energy through thermal bremsstrahlung, provides a key channel through which clusters can be observed in the X-ray band, with its exact properties being determined by a combination of astrophysical processes and cosmology. It was originally thought that left unabated, high central radiative losses would result in the formation of a cooling flow (Fabian 1994), leading to large quantities of cold gas dropping out of the ICM and hence massive molecular reservoirs and high star formation rates (SFRs) within cluster centres. However, cluster cores typically show only moderate molecular gas content (McNamara et al. 2014; Russell et al. 2014, 2017, 2019; Fogarty et al. 2019; Castignani et al. 2020) and low SFRs (McNamara & Nulsen 2007; Donahue et al. 2010; Cooke et al. 2016; Fogarty et al. 2017), although there can be exceptions (e.g. Crawford et al. 1999; Egami et al. 2006; Von Der Linden et al. 2007; Mittal, Whelan & Combes 2015; Fogarty et al. 2017). Additionally, direct observational signatures of the cooling, such as X-ray emission lines below 1 keV, are absent from numerous observations (e.g. Ikebe et al. 1997; Makishima et al. 2001; Peterson et al. 2001, 2003; Tamura et al. 2001; Böhringer et al. 2002; Lewis, Stocke & Buote 2002; Matsushita et al. 2002) and ultraviolet (UV) detections of O VI lines suggest lower than expected cooling rates (e.g. Oegerle et al. 2001; Bregman et al. 2006; Donahue et al. 2017).

Instead, heating by active galactic nuclei (AGN) feedback, possibly in combination with other physical processes (for example, thermal conduction, Ruszkowski & Begelman 2002; Zakamska & Narayan 2003; Voigt & Fabian 2004; Conroy & Ostriker 2008; Bogdanović et al. 2009; Ruszkowski & Oh 2010, 2011; Yang & Reynolds 2016a), is very likely to be able to keep the ICM warm and prevent the formation of the aforementioned cooling flows (see e.g. McNamara & Nulsen 2007, 2012; Fabian 2012, for reviews). Specifically, accretion on to the central black hole (BH) can produce jets that inflate lobes of relativistic gas in cluster cores (Binney & Tabor 1995; Omma et al. 2004; McNamara et al. 2005; Fabian et al. 2006; Sijacki & Springel 2006a; Cattaneo & Teyssier 2007; Forman et al. 2007; Sijacki et al. 2007; Dubois et al. 2010, 2012), which exhibit themselves as cavities in X-ray observations (Fabian et al. 2000, 2011; McNamara et al. 2000; Heinz et al. 2002; Forman et al. 2007). Along with radio emission, such cavities (or bubbles) are found to be common in cool core clusters, i.e. those with short central cooling times (Burns 1990; Dunn, Fabian & Taylor 2005; McNamara & Nulsen 2007; Dunn & Fabian 2006, 2008; Sun 2009; Fabian 2012; Hlavacek-Larrondo et al. 2012). Combining this with the observed correlation between the energy content of the lobes (based on PV estimates) and the ICM cooling losses (based on X-ray luminosity, e.g. Rafferty et al. 2006; McNamara & Nulsen 2007; Nulsen et al. 2007; Dunn & Fabian 2008; Fabian 2012) strongly suggests that ICM cooling and heating is coupled to the AGN activity.

While the feedback energetics are sufficient to offset the expected ICM cooling losses, and despite a growing consensus that jet feedback in galaxy clusters should be gentle and continuous (see discussions in McNamara & Nulsen 2007, 2012; Fabian 2012; Werner et al. 2019, for examples), there is still some debate with regard to the mechanism (or combination thereof) through which the mechanical energy of the jets is effectively and (largely) isotropically communicated to the ICM. Perhaps the clearest heating mechanism seen in observations is via the weak shocks driven into the ICM by the lobe inflation (Fabian et al. 2003; Randall et al. 2015) that are also seen in simulations of jet feedback (Yang & Reynolds 2016b; Bourne & Sijacki 2017; Li, Ruszkowski & Bryan 2017; Martizzi et al. 2019). Such shocks should slow as they propagate into the ICM and broaden into sound waves that could transport energy to larger radii, with some works interpreting observed ripples in X-ray brightness maps as sounds waves (Fabian et al. 2003, 2005, 2017). Although simulations also exhibit sound wave production (Ruszkowski, Brüggen & Begelman 2004; Sijacki & Springel 2006a; Bourne & Sijacki 2017; Bambic & Reynolds 2019), the volume over and rate at which they heat the ICM will depend upon the ICM viscosity set on micro-physical scales (Ruszkowski et al. 2004; Sijacki & Springel 2006b; Zweibel et al. 2018). Other mechanisms provide less direct or obvious observational signatures, such as instability-driven mixing (Hillel & Soker 2016, 2017), cavity heating (Churazov et al. 2002; Bîrzan et al. 2004), jet-driven turbulence (Zhuravleva et al. 2014; Zhang, Churazov & Schekochihin 2018), cosmic ray (CR) production (Sijacki et al. 2008; Pfrommer 2013; Ehlert et al. 2018), and gas circulation (Yang & Reynolds 2016b).

Given that a number of these processes do not leave a directly or easily observable signature, the role of simulations is critical to determine how they manifest themselves in cluster observables and thus ultimately resolve what processes dominate the heating budget. Previous simulation works typically sit within one of two categories. On the one hand, cosmological simulations self-consistently capture the buildup and evolution of galaxy clusters over cosmic time and hence by construction exhibit the cluster properties and dynamics expected from large-scale structure formation theory (e.g. Rasia et al. 2015; Bahé et al. 2017; Barnes et al. 2017a, 2017b; Hahn et al. 2017; McCarthy et al. 2017; Barnes et al. 2018; Henden et al. 2018). While such simulations provide realistic cluster environments, their large spatial and temporal dynamic range means that they typically lack sufficient resolution in cluster cores, close to the central BH, to include detailed modelling of the AGN feedback process. On the other hand, previous works including physically more realistic models of jet injection that self-consistently capture the lobe inflation and its interaction with the ICM are restricted in their dynamic range by the requirement of high resolution within the jets and lobes and hence adopt idealized ICM setups (e.g. Reynolds, Heinz & Begelman 2001; Basson & Alexander 2003; Krause 2003; Omma et al. 2004; Vernaleo & Reynolds 2006; Gaibler, Krause & Camenzind 2009; Hardcastle & Krause 2013, 2014; English, Hardcastle & Krause 2016; Bourne & Sijacki 2017; Weinberger et al. 2017). Up until recently, only a handful of examples of an AGN jet feedback model included in a cosmological environment exist in the literature (Heinz et al. 2006; Morsony et al. 2010; Mendygral, Jones & Dolag 2012). Our previous study (Bourne, Sijacki & Puchwein 2019) builds upon these early works by combining very high-resolution jets and uniquely a dedicated refinement scheme to ensure continued high resolution in the jet lobes and at the lobe-ICM interface, within a fully cosmological cluster that was evolved with state-of-the-art physical subgrid models. Here, we further expand on this work by including jets of different power and perform a detailed analysis of the lobe evolution and cluster heating processes with the aim to address the crucial issues of how the jet energy is communicated to the ICM and if it is sufficient to hinder cooling flows.

The paper is structured as follows. In Section 2, we outline the numerical implementation and initial condition (IC) generation and properties. Section 3 presents our main results and is split into sections covering lobe evolution (Section 3.1), how energy is communicated to the ICM (Section 3.2), and the impact on cluster properties (Section 3.3). Finally, in Section 4, we discuss our results and limitations, and in Section 5, we provide a short summary of our key results.

2 NUMERICAL METHOD

2.1 Simulation overview

Simulations are performed using the moving mesh code arepo (Springel 2010), which employs a finite-volume approach to solve for the hydrodynamics on an unstructured Voronoi mesh that moves with the flow. Our ICs are identical to those used in our previous work, Bourne et al. (2019), and comprise a cosmological zoom simulation of a galaxy cluster evolved from redshift z = 127 down to z ≃ 0.1 using subgrid models almost identical to those employed for the original Illustris project (Genel et al. 2014; Vogelsberger et al. 2014; Sijacki et al. 2015), barring a small difference in the radio-mode AGN feedback model, and we point the reader to Section 2.5 for a more detailed account of the generation of the ICs, along with the initial properties of the cluster.

Subsequently, instead of using the radio feedback model that was employed up until z = 0.1 (Sijacki et al. 2007, 2015) where AGN lobes are not inflated self-consistently, we activate the jet feedback model presented in Bourne & Sijacki (2017) and discuss in Section 2.2 to study a range of different jet powers. Building on Bourne et al. (2019), as well as a fiducial jet power of 3.9 × 10⁴⁴ erg s⁻¹, we perform three additional simulations with low (7.8 × 10⁴³ erg s⁻¹) and high (1.94 × 10⁴⁵ erg s⁻¹) power jets and one with no jets for comparison. All jets are active for 20 Myr and their properties are summarized in Table 1. We follow the subsequent jet lobe evolution for a further ∼80 Myr. The jet model relies on the super-Lagrangian refinement (SLR) technique of Curtis & Sijacki (2015) combined with additional modifications, which are discussed in Section 2.3. We note that only the most massive BH in the cluster of interest is active during our simulations, and unlike the original simulations used to produce the ICs, for simplicity we do not include the effects of an AGN radiation field in the heating and cooling prescription. Additionally, while during the generation of the cluster ICs, we followed the supernovae wind model used in the Illustris simulations, where the mass loading, |$\eta _{\rm w}=\dot{M}_{\rm w}/\dot{M}_{\rm *}$|⁠, is a function of the halo 1D velocity dispersion (see Vogelsberger et al. 2013, for model details); for simplicity, we instead use a fixed value of η_w = 1 for the post z = 0.1 simulations presented here. All other subgrid models and parameters remain unchanged.

2.2 Jet feedback model

Two tracers are used to track jet material. The jet tracer, |$f_{\rm J}=m_{\rm J}^{\rm trace}/m_{\rm cell}$|⁠, is defined for each cell as the fraction of mass within that cell that has been within the jet cylinder during an injection episode, i.e. we set f_J = 1 whenever a cell is in the jet cylinder and the jet is active. The tracer mass, |$m_{\rm J}^{\rm trace}$|⁠, is subsequently advected between cells in line with the total gas mass. We define lobe material as any non-star-forming gas (n < 0.26 h² cm⁻³) for which |$f_{\rm J}\gt f_{\rm J}^{\rm thresh}$|⁠. In Bourne et al. (2019), we noted that while there is some freedom in the choice of |$f_{\rm J}^{\rm thresh}$|⁠, a cell’s internal energy should be dominated by jet material when f_J ≥ T_ICM/T_J, where T_ICM and T_J are the typical ICM and jet-rich lobe material temperatures, respectively. For our fiducial jet power, |$f_{\rm J}^{\rm thresh}=10^{-2.5}$| provides a threshold that coincides with the point at which gas density (temperature) starts to decrease (increase), transitioning from ICM to lobe-like values. However, higher (lower) power jets result in clearly higher (lower) temperature gas by roughly a factor of a few. As such, we adjust values of |$f_{\rm J}^{\rm thresh}$| for different jet powers (see Table 1), noting that this is not an exact art and as such the values chosen reflect the approximate transition between jet and non-jet material. Additionally, we inject a second tracer that tracks the mass injected via equation (2), |$m_{\rm J}^{\rm inj}$|⁠, which like our primary tracer is advected between cells.

2.3 Refinement criteria

2.4 Gas draining

We drain gas at a rate of |$\dot{M}_{\rm in}$| from certain cells within h_bh to ensure energy conservation, in a method similar to Bourne & Sijacki (2017), which is adapted from previous methods (Vogelsberger et al. 2013; Curtis & Sijacki 2015). During a BH time-step of size dt, a mass of |$M_{\rm drain}=\dot{M}_{\rm in}\times dt$| is removed from the cells that lie within a torus-shaped region outside of the jet cylinder, aligned perpendicular to the jet direction and with an opening angle of (π − θ_J), in a mass-weighted fashion. Subsequently, the BH mass is increased by |$\Delta M_{\rm bh}=\dot{M}_{\rm bh}\times dt$| [see equation (4)] and a mass of |$\Delta M_{\rm J}=\dot{M}_{\rm J}\times dt$| is added to the cells within the jet cylinder following the same weighting as the energy and momentum injection. Note that this method is designed to conserve energy, not mass, as a fraction ϵ_r/(1 + η_J) of M_drain is assumed to go into powering the jet during each time-step.

2.5 Cluster initial conditions

The ICs for this work are taken from a cosmological zoom simulations of a galaxy cluster evolved to low redshift, assuming the WMAP-9 cosmological parameters (Hinshaw et al. 2013) Ω_m = 0.2726, Ω_Λ = 0.7274, Ω_b = 0.0456, σ₈ = 0.809, n_s = 0.963, and H₀ = 70.4 km s⁻¹ Mpc⁻¹ = h × 100 km s⁻¹ Mpc⁻¹, as in the original Illustris simulations (Genel et al. 2014; Vogelsberger et al. 2014). The zoom-in region encompasses material that resides within a volume spanning ∼30 h⁻¹ cMpc across. High-resolution gas cells have a target mass, |$m_{\rm cell}^{\rm target}=1.37\times ~10^{7}$| h⁻¹ M_⊙ and adaptive softening lengths equal to 2.5 × r_cell, while dark matter particles have a mass of m_dm = 4.83 ×  10⁷ h⁻¹ M_⊙ with physical gravitational softening lengths set to 2.8125 h⁻¹ kpc below z = 5.

Subgrid models for gas radiative processes (Katz, Weinberg & Hernquist 1996; Wiersma et al. 2009b), interstellar medium (ISM) and stellar physics (Springel & Hernquist 2003; Springel 2005; Wiersma, Schaye & Smith 2009a; Vogelsberger et al. 2013), and BH processes (Di Matteo, Springel & Hernquist 2005; Springel 2005; Sijacki et al. 2007) are almost identical to those employed in the original Illustris project. In producing the ICs, the radio-mode model of Sijacki et al. (2007) is employed during low Eddington rate phases of BH growth. The feedback acts by injecting hot bubbles of gas (to mimic radio lobes) around accreting BHs whenever a BH’s mass increases by a fraction δ_BH. For the model parameters chosen in the original Illustris project, namely δ_BH = 0.15, too much gas was ejected from groups and clusters (Genel et al. 2014). Therefore, in this work we apply a gentler, albeit more frequent bubble injection by setting δ_BH = 0.015, resulting in a gas mass fraction of |$f_{\rm gas}^{500}=14{{\ \rm per\ cent}}$| within R_500, c for our ICs. Key physical properties of the cluster calculated by the on-the-fly friends of friends (FoF) and subfind algorithms (Davis et al. 1985; Springel et al. 2001; Dolag et al. 2009) at z ≃ 0.1 are provided in Table 2. Our cluster IC is visualized in the central panel of Fig. 1, which shows the temperature encoded by colour and column density encoded by brightness projected through a 15 h⁻¹ Mpc on a side box centred on the central BH of the cluster of interest. Additionally, the radii r_500, c = 763 h⁻¹ kpc and r_200, c = 1178 h⁻¹ kpc are indicated by the dotted–dashed and dashed lines, each of which contains total masses of 2.82 × 10¹⁴ h⁻¹ M_⊙ and 4.14 × 10¹⁴ h⁻¹ M_⊙, respectively. The figure shows the main cluster as the large region of hot ICM sitting at the nodes of a major cold filament linking to a smaller cluster in the top left of the image as well as many smaller systems along filaments in other directions that contain many galaxies.

Figure 1.

Overview of cluster environment and jet evolution. The large central panel shows a column-density/temperature map of the main cluster, projected through a box of 15 h⁻¹ Mpc on a side, with the colour indicating temperature and brightness indicating the gas density. On the left-hand side are three smaller stacked panels that show a face-on temperature projection of the cold disc structure constituting the BCG at the centre of the cluster (top), side-on column density map of the same structure (middle), and a temperature projection of a cold substructure that falls into the cluster centre and interacts with the jet (bottom). The temperature projection of the BCG is overlaid with velocity streamlines while the density projection shows line-of-sight velocity contours, illustrating the rotation of the cold structure. To the right of the main panel is a grid of six mock observations of different power jets (low to high power, from top to bottom) at ∼20 Myr. The left-hand column shows composite X-ray images, in which the different colours represent different energy bands and the right-hand column shows X-ray/radio composite images in which the blue shows the high-energy X-ray band and the pink shows the synchrotron emissivity as a proxy for the radio band.

Open in new tabDownload slide

The cluster and epoch were chosen as they show no signs of a recent AGN outburst and hence provide a relatively ‘clean’ environment into which we can inject the jets. Multiple gas phases exist in the simulations including ISM, warm, ICM, and jet lobe material. The top- and middle left-hand panels of Fig. 1 show temperature and density projections of a ∼4 × 10¹⁰ M_⊙ cold gas structure, face-on and edge-on, respectively. Its disc-like structure presents distortions that are particularly evident in the density projection and has similar properties to the cold gas observed as optical emission line nebula in the centres of some galaxy clusters (e.g. Hamer et al. 2014, 2016, 2019). The structure rotates about the potential minimum of the cluster, as illustrated by line-of-sight velocity contours shown in red (away) and blue (towards) in the middle left-hand panel, and the streamlines of the in-plane velocity shown in the top left-hand panel.

The chosen jet direction along the z-axis is not perfectly perpendicular to the plane of the cold disc, although this outcome is exacerbated by the distorted disc shape. The jet propagation and lobe properties are subsequently altered by jet–ISM interaction (Bourne et al. 2019). Additional cold material exists beyond the central cold disc, which contributes to the population of structures that pass near the brightest cluster galaxy (BCG). The bottom left-hand panel of Fig. 1 shows a temperature projection of an infalling system comprising multiple cold clumps enshrouded by a bow shock of hot gas. The system initially resides ∼65 h⁻¹ kpc from the cluster centre and interacts strongly with the bottom jet lobe at late times.

3 RESULTS

3.1 Lobe evolution

3.1.1 Overview

High jet velocities (peaking at 0.56 c, see Table 1) result in strong shocks and lobe material with temperatures of up to T ∼ 10¹¹ K and densities that can be ∼1000 times lower than the ICM density. This gives rise to observable signatures similar to real clusters such as X-ray cavities surrounded by X-ray bright rims and radio emission from jet material. Mock X-ray and radio observations for different jet powers at ∼20 Myr are shown in the right-hand panels of Fig. 1. X-ray data are produced using pyXSIM (ZuHone et al. 2013) assuming a fixed metallicity of 0.3 solar and that the cluster is at the same redshift as the Perseus cluster (Fabian et al. 2006). Given that we use an effective equation of state to model the ISM and do not include molecular cooling, we are unable to reliably capture the thermal properties of cold dense gas. Therefore, along with setting an upper temperature limit of 64 keV, we apply the temperature density cut of Rasia et al. (2012) to remove gas satisfying |$T_{\rm keV} \lt 3\times ~10^6\rho ^{0.25}_{\rm cgs}$|⁠. RGB composite images are shown in the second column from the right, with each colour representing a different X-ray energy band (0.5–1.2, 1.2–2, and 2–7 keV, in red, green, and blue, respectively). Cooler ICM gas exists in the vicinity of the central disc structure and the previously mentioned infalling substructure, clearly seen in the low energy band. Cavities surrounded by X-ray bright rims are seen in all three cases but are most obvious in the high power jets, which show large, elongated lobes and clear rims around both the top and bottom lobes. On the other hand, the low power jets produce much smaller and rounder lobes, albeit still with clearly visible rims. Interestingly, while still showing rims around both cavities, for the fiducial and high power jets the top lobe rim is thinner and fainter than for the low power jets.

The right-hand column zooms in to the cocoon region (indicated by the dashed lines) and shows the high-energy X-ray band (blue) overlaid by estimated synchrotron emissivity of jet material (pink). We calculate expected synchrotron emissivity³ following Appendix  B of Bicknell et al. (2018), assuming a random magnetic field and that magnetic and relativistic electron energy densities are fixed fractions of the gas internal energy density. There are clear radio bright jet structures along with more diffuse material filling jet lobes. The low power jets have a more distorted morphology and appear bright throughout. The fiducial power jets appear to contain several hotspots, while the high power jets exhibit clear hotspots at their ends. While the high power jets appear FR-II like (Fanaroff & Riley 1974), the low power jets are similar to typical FR-I sources and the fiducial power jets sit somewhat in between, with its extended top lobe but rounder bottom lobe. Ehlert et al. (2018), who performed jet simulations in an idealized environment, found a transition at a similar power to us, although the high-resolution study of Massaglia et al. (2016) found a transition at a lower power of P_J ∼ 10⁴³ erg s⁻¹. While Ehlert et al. (2018) found that they could achieve a transition at lower power with improved resolution, we note that at least some of the morphological features observed in our simulations are driven by the interaction of the jet with the cluster environment, most clearly illustrated by the fiducial power run in which the interaction of the bottom jet lobe with the cold disc drives its morphological difference to the top lobe. The radio images also show fluctuations along the jet axis (especially clear in the high and fiducial power jets), consistent with the internal shocks that occur along the jet axis (for further details, see Section 3.2.2).

3.1.2 Lobe dynamics

The long-term evolution of the jets is shown for all three jet powers in Fig. 2, which comprises volume-rendered images of the jet lobes (⁠|$f_{\rm J}\gt f_{\rm J}^{\rm thresh}$| and n < 0.26 h² cm⁻³) as illustrated by the jet tracer in green, while red colours indicate warm and cold ICM material. Each horizontal collection represents a different jet power and contains eight small panels that, from left to right, show the evolution of the system, with the top and bottom rows showing views that are rotated by 90° about the z-axis with respect to each other. The larger panel in each collection shows a clearer view of the jet lobes at 49.4 Myr along with the in-plane velocity field shown by the coloured arrows. At early times, we see the classical double-lobe structure emerging above and below the central galaxy, with injected jet material residing largely along the jet axis. By the time the jet switches off, there are clear asymmetries between the bottom and top lobes for all three jet powers, with the top lobes typically being elongated and rather narrow, while the bottom jet lobes are shorter and more round, somewhat similar to Abell 4059 (Heinz et al. 2002) and Abell 2052 (Blanton et al. 2001). These differences are driven by the interaction of the bottom jet lobes with the cold gas (Bourne et al. 2019).

Figure 2.

Volume-rendered images of the lobes as illustrated by the jet tracer in green, as well as ISM and warm phase material shown in red. The images are organized into three main groups for each jet power, going from low-to-high power from top to bottom, respectively. Each group consists of eight small panels, which show four different times (left to right) for two line-of-sight orientations (rotated by 90° about the z-axis with respect to each other), and a large panel to the right of each group, which shows an enlarged version of the panel at 49.4 Myr in the lower row, with velocity field lines overlaid. The figure illustrates how cold material interacts with the jet material, in particular, the southern lobe, which is disrupted by the passage of a cold structure. The figure also shows the relative survivability of the lobes of different powered jets, with higher power lobes being more resilient to the cluster weather and surviving longer.

Open in new tabDownload slide

An alternative comparison of the lobe dynamics for different jet powers is shown in Fig. 3. The plot shows the evolution of the mass-weighted average position of lobe material for the top and bottom lobes in the low (red), fiducial (green), and high (blue) power jet simulations, with averages calculated using the total cell masses over all cells with |$f_{\rm J}\gt f_{\rm J}^{\rm thresh}$| for the relevant jet power. Starting at the origin, the lines get darker with time up to a maximum value of 100 Myr, with each segment representing the value from a single snapshot and the segment width being the mass-weighted standard deviation of the z − y positions of cells within the jet lobes. All three jet powers show somewhat similar behaviour as the lobe material initially travels almost directly vertically along the trajectory of the jets, with the distance travelled increasing with jet power. Once the jet turns off, the top lobe material is displaced to the right by the ram pressure of bulk ICM motions, which is illustrated by the velocity streamlines in the large panels on the right-hand side of Fig. 2. For the fiducial and high power jets, this displacement accompanies the continued buoyant rise of the top lobes. However, for the low power jet, the mean z position of the top lobe begins to decrease. Comparing vertical and horizontal motions, we find that the horizontal motion becomes increasingly more important than the buoyant motion of the top lobe for lower power jets. After ∼40–50 Myr, a shredded cold substructure directly interacts with the bottom lobe first compressing and then disrupting it. We find that the average trajectory of the bottom lobe material in the low and fiducial power jets is completely reversed, while for the high power jet, its progress appears to stall. The jet lobe sizes and lifetimes increase with jet power while the low power jet lobes are completely absent from the 100-Myr panel, and very little of the fiducial power jet lobes remain (mostly the top lobe), having been completely destroyed by the ICM weather and mixed in with the ambient medium, the top lobe of the high power jet run still appears prominently along with the diffuse remains of the bottom lobe. Overall, these two figures illustrate the significant role that the inherent ICM velocity field can have on the trajectories of jet lobes. Such displacement of lobes due to ICM bulk motions has been seen in previous works (Sijacki, Springel & Haehnelt 2009; Bourne & Sijacki 2017) and could provide an explanation as to why relic lobes appear along different trajectories without invoking large-angle jet precession (see e.g. Dunn, Fabian & Sanders 2006a; Babul, Sharma & Reynolds 2013).

Figure 3.

Overview of the dynamical evolution of the jet lobes for different jet powers, showing the mass-weighted mean position of each jet lobe for the low (red), fiducial (green), and high (blue) power jet runs. Points are calculated for every snapshot and become darker with time, which, on average, are at 1.3-Myr intervals and in total span a 100-Myr period. The width of each point shows the dispersion about the mean lobe cell positions, with thicker lines indicating a more extended lobe.

Open in new tabDownload slide

3.1.3 Lobe energetics

Tracking the energy content of lobes and the PdV work done to inflate them can provide insight into ICM heating mechanisms and allow comparison with observational studies that estimate jet energetics from lobe properties (e.g. Rafferty et al. 2006; McNamara & Nulsen 2007; Nulsen et al. 2007; Dunn & Fabian 2008; Fabian 2012). The main panels in Fig. 4 show lobe energy evolution for each jet power, normalized by the total injected energy (⁠|$\int _{0}^{20\:{\rm Myr}}P_{\rm J}dt$|⁠). The combined lobe thermal plus kinetic energy is shown by the solid black line; the individual thermal, total kinetic, and turbulent kinetic energies are shown by the blue, red, and green solid lines, respectively; and the total injected energy is shown by the dotted black line. Estimated lobe losses due to PdV work and mixing are shown by the dashed cyan and magenta lines, respectively. Quantities are calculated following the method in Bourne et al. (2019) assuming jet bulk velocity cutoffs given by |$v_{\rm J}^{\rm mean}$| defined in Table 1. The grey shaded region indicates the period over which the jet is active and the horizontal dotted lines indicate 1, 5, 10, and |$50{{\ \rm per\ cent}}$| of the total injected energy. Additionally, the three lower panels show the evolution of other lobe properties of interest, namely the ratios of the lobe energy to injected energy (E_lobe/E_inj, solid line), tracer mass in the lobe to total tracer mass (⁠|$M_{\rm J}^{\rm lobe}/M_{\rm J}^{\rm tot}$|⁠, dashed line), total lobe mass to maximum lobe mass (⁠|$M_{\rm lobe}/M_{\rm lobe}^{\rm max}$|⁠, dot–dashed line), and lobe volume to maximum lobe volume (⁠|$V_{\rm lobe}/V_{\rm lobe}^{\rm max}$|⁠, dotted line).

Figure 4.

Evolution of jet lobe energy content is shown in the top panel, split into the thermal (blue), total kinetic (red), and turbulent (green) components. The total injected energy is shown by the dotted line, while the total (thermal + kinetic) energy content of jet lobe material (⁠|$f_{\rm J}\gt f_{\rm J}^{\rm thresh}$|⁠) is shown by the solid black line. The dashed cyan and magenta lines show estimated lobe energy loss due to PdV work and mixing, respectively. While the jet is active, approximately half of the total injected energy remains in the lobe and approximately half goes into the ICM. The grey shaded region indicates the period over which the jet is active. The thermal component dominates the lobe energy, while the turbulent component dominates the total kinetic energy. In the bottom panel, we show the evolution of four mixing-related quantities, namely the lobe energy normalized by the total injected energy (solid), the injected jet mass within the lobe normalized by the total injected mass (dashed), the total lobe mass normalized to its maximum value (dot–dashed), and the lobe volume normalized by its maximum value (dotted).

Open in new tabDownload slide

During the inflation phase, we find that the combined thermal and kinetic energy of the lobes generally accounts for approximately half of the injected jet energy (see also e.g. Hardcastle & Krause 2013; Bourne & Sijacki 2017; Weinberger et al. 2017), with values of 51.7, 48.9, and |$46.3{{\ \rm per\ cent}}$| for the low, fiducial, and high power jet lobes, respectively, shortly before the jets switch off. Large departure from this occurs only in the first few Myr of jet injection, where the fraction of injected energy in the lobes peaks at |$70-80{{\ \rm per\ cent}}$| before significant expansion and mixing have occurred. For all three jet powers, the thermal component dominates the energy budget because of the shocks driven by the fast kinetically dominated jets rapidly heating the lobe gas. The residual kinetic energy accounts for only |$\sim 5-6{{\ \rm per\ cent}}$| of the jet energy residing in the lobes by 20 Myr, which itself is dominated by the turbulent component that accounts for |$\sim \! 70-80{{\ \rm per\ cent}}$| of the lobes’ total kinetic energy.

Because lobes retain only half of the jet energy during inflation, the rest must be communicated to the ICM. Based on our simulations, we find that mixing is subdominant during the inflation phase, contributing between |$8{{\ \rm per\ cent}}$| of the injected energy from the low power jet down to |$2.6{{\ \rm per\ cent}}$| for the high power jet. On the other hand, the PdV work done to inflate the lobes accounts for between |$37.4{{\ \rm per\ cent}}$| and |$42.5{{\ \rm per\ cent}}$| of the injected energy from low to high jet powers, respectively. The amount of PdV work done on the ICM and how this energy manifests itself within the ICM depends on how rapid or ‘explosive’ the lobe inflation is (see section 5.4 of McNamara & Nulsen 2007). If inflation is slow (adiabatic), the ICM can be compressed and sound waves can be driven. If inflation is fast, shocks are driven into the ICM, generating entropy as kinetic energy is irreversibly dissipated as heat. Fig. 4 shows that the fraction of injected energy going into PdV work increases with jet power, in line with theoretical expectations.

At the end of inflation, the ratio of the instantaneous PV estimate, i.e. using the instantaneous lobe pressure and volume, and integrated PdV estimate are PV: ∫PdV = 0.82, 0.72, and 0.64 for the low, fiducial, and high jets powers, respectively. Therefore, PV ∼ ∫PdV holds only if the lobe pressure remains roughly in equilibrium with the ICM pressure, i.e. for slow adiabatic inflation of the lobes. This assumption breaks down at higher powers and impacts our ability to estimate jet powers and the energies needed to inflate lobes from observations. The instantaneous lobe enthalpy, |$H=E_{\rm lobe}^{\rm therm}+PV=\gamma PV/(\gamma - 1)$|⁠, is often taken as a measure of the jet energy needed to inflate a lobe. We find that the ratio of |$H:E_{\rm J}^{\rm inj}=0.77$|⁠, 0.73, and 0.68 for the low, fiducial, and high power jets, respectively, i.e. the lobe enthalpy underestimates the energy needed to inflate the lobes by at least |$23{{\ \rm per\ cent}}$| in the low power case and |$32{{\ \rm per\ cent}}$| for the high power jets, suggesting that lobe enthalpy sets a lower limit on the required energy.

While the partition of energy for different jet powers is quite similar during lobe inflation, differences become more apparent after jets switch off. All runs show a sudden drop in the kinetic energy component as the jet action stops, accompanied by a drop in the thermal energy component, which is no longer replenished by shocked jet material. Significant differences in the fraction of energy lost due to mixing begin to appear between the different jet powers, with the lower power lobes being more susceptible to mixing than the higher power lobes. Although it has been suggested that jet-driven instabilities can effectively mix lobe material with the ICM (Hillel & Soker 2016, 2017), we previously showed that by implementing appropriate refinement criteria our method results in limited instability and numerical mixing (Bourne & Sijacki 2017). We therefore stress that the efficient mixing seen at late times is predominantly driven by ICM and substructure motions disrupting the lobes. By 100 Myr, we estimate that while only |$23.4{{\ \rm per\ cent}}$| of the high power jets' energy has been mixed into the ICM, this jumps up to |$71{{\ \rm per\ cent}}$| in the case of the low power jets, suggesting that the lower power jet is more susceptible to the cluster weather. While different values of |$f_{\rm J}^{\rm thresh}$| chosen for each jet power can impact the mixing rate, a more detailed mixing analysis in Section 3.2.3 shows that our conclusion is still valid. The stark differences in the levels of mixing lead to stark differences in the amount of energy retained by the jet lobes, with the low power jet lobes retaining only |$2.4{{\ \rm per\ cent}}$| of the injected energy by 100 Myr, while the high power jet retains |$28.6{{\ \rm per\ cent}}$|⁠.

Finally, while the cumulative PdV remains flat after the high power jets switch off, it decreases for lower power jets, suggesting that work is done on the jet lobes. From the lower panels, peak lobe volumes occur between a few to ∼15 Myr after the jet action halts, while peak lobe masses occur later. These peaks occur earlier for lower power jets and the subsequent declines appear steeper when compared to the higher power jets. Continued entrainment and hence gain in mass even after lobe volumes begin to decrease suggest that the lobes are compressed, leading to the decrease in cumulative PdV losses. The effect is particularly evident around ∼40 Myr, where lobe thermal energy shows a distinct bump and corresponding dip in PdV losses. This is due to the infalling substructure (see Section 3.1.2) impacting the bottom lobe. However, compression is not the only mechanism that can effect our PdV estimate. Once lobes begin to lose significant amounts of mass, the corresponding volume of the lost material will also impact our dV estimates and as such, we suggest that our estimates begin to break down as the system evolves beyond the mass peak.

3.2 Energetics: from jet to environment

3.2.1 Where does the energy go?

Differences in the energy content of any two systems are a combination of energy directly injected by the jet and differences in energy gains/losses that arise from the other physical processes, such as radiative cooling, star formation, and stellar feedback that are themselves impacted by the jet action. Therefore, the energy difference between a jet and non-jet run could exceed the total energy input of the jet, as is the case for the high power jet run (see Table 3). In addition to the non-linear nature of the physical processes and their inter-dependant evolution, the stochastic nature of simulations and the effect that small perturbations (such as launching a small-scale jet) can have on the system⁴ can also result in differences between runs (see e.g. Keller et al. 2019). While such differences appear as only very small fluctuations when compared to the whole system, they can be large compared to the total injected jet energy.

To explore where the jet energy ultimately resides in the system, we compare the evolution of the fractional differences in mass and energy profiles of the clusters with and without jets in Fig. 5. The rows, from top to bottom, show the fractional difference in mass, thermal energy, and kinetic energy as a function of radius, and from left to right for the low, fiducial, and high power jets, respectively. Values calculated for low/high power jets have been doubled/halved to fit all three jet powers on the same y-axis and all curves have been smoothed using a Gaussian filter with a standard deviation of 1.5 h⁻¹ kpc to reduce noise. The different coloured lines refer to different times, while the solid and dotted lines show values calculated with and without jet material included, respectively. All quantities are calculated within radial bins of fixed width dr = 1.5 h⁻¹ kpc and the small vertical lines at the bottom of each panel define an influence radius, r_out, calculated as the outer radius after the main peak in δU/U, i.e. the radius at which it returns to zero.

Figure 5.

Changes in radially averaged quantities between jet and non-jet runs. From top to bottom, the rows show the radial fractional difference in mass, internal energy, and kinetic energy, with each jet power, from low to high, shown in each column from left to right, respectively. The profiles are calculated at different times and shown in different colours, as given by the key, and are presented both with and without jet lobe material included, by the solid and dotted lines, respectively. Additionally, the low/high power run quantities have been doubled/halved to allow all three jet power runs to be plotted on the same y-axis and we have defined a time-dependent influence radius, r_out, indicated by the small ticks at the base of each panel and defined by the outer edge of the main δU/U peak.

Open in new tabDownload slide

The qualitative behaviour is similar for all fractional changes, albeit with higher power jets exhibiting more well-defined, larger (noting the re-normalization), broader, and longer lived changes. The inflation and subsequent rise of the lobes sweep up and displace gas to larger radii, as indicated by the positive peak in δM/M moving outwards and flattening with time and reduced values at small radii. Similarities in the δM/M curves with and without jet material illustrate that the displaced mass is dominated by non-jet material, with higher power jets displacing more mass. At early times, the main peaks exist at |$r\lower.5ex\hbox{$\,\,\buildrel\lt \over \sim \,\,$}r_{\rm lobe}^{\rm max}$|⁠, where |$r_{\rm lobe}^{\rm max}$| is the maximum radius of lobe material; however, at later times (particularly for low and fiducial power jets) the peaks occur beyond |$r_{\rm lobe}^{\rm max}$| and may be related to the weak shock/sound wave generated by the initial lobe inflation. While in some cases, particular for low power jets, small fluctuations are present beyond the main peak, they are very subdominant compared to the apparent influence of the jet.

Mass displacement by jet feedback leads to a net gain in the potential energy of baryonic material within the central region of the cluster. To quantify and minimize the effect that small fluctuations in δM/M seen at large radii may have on the net difference in potential energy of a given jet run compared the non-jet run, ΔE_pot, we perform averaging over spheres. Specifically, we calculate ΔE_pot within multiple spheres that have outer radii ranging from r_out to 1.25 × r_out at 1.5 h⁻¹ kpc intervals (matching the bin widths in Fig. 5) and take the standard deviation in these values as the associated error.

Average values of r_out and |$\Delta E_{\rm pot}/E_{\rm tot}^{\rm inj}$| at 90.2 Myr are shown in the second and final columns of Table 3, respectively. The associated error is large for the low power jets and decreases for higher powers, although this may be expected, given that the values are a percentage of the total injected energy. However, despite the large error for the low power jet, the average value is still consistent with that of the fiducial and high power jets, which implies that systems including jets exhibit potential energies that are larger than those of non-jet systems by |$\lower.5ex\hbox{$\,\, \buildrel\gt \over \sim \,\,$}50{{\ \rm per\ cent}}\, E_{\rm inj}^{\rm tot}$|⁠.

Previous simulations found smaller values, for example in (Bourne & Sijacki 2017), we found that for a non-radiative run |$\sim \! 16{{\ \rm per\ cent}}$| of the inject energy ends up as potential energy (see also Hardcastle & Krause 2013); however, these simulations only considered values while the jet is active and did not capture the buoyant rise and further displacement of mass after the jets switch off. On the other hand, Weinberger et al. (2017) performed simulations that tracked the post-inflation rise of the lobes and found that |$\sim \! 40{{\ \rm per\ cent}}$| of the jet energy goes into the potential energy of the system, much closer to the values found here. While it has been suggested that most of the jet energy ultimately ends up in potential energy of the system (e.g. McNamara & Nulsen 2007; Mathews 2009; Mathews & Guo 2011), it is important to remember that we measure the net difference between two systems (jet versus non-jet) at one instance in time, capturing the combined effect of the lobe inflation displacing mass outwards and the lack of jet heating possibly resulting in mass flowing inwards. Instead, measuring the change in the potential energy within r < 150 h⁻¹ kpc between t = 0 and t = 90.2 Myr separately for each run, we find net losses in potential energy irrespective of the inclusion of jet feedback, because ultimately the jets only act to reduce the net inflow of gas (see also Guo, Duan & Yuan 2018).

Considering all gas (solid lines), the similarities between the troughs in δU/U and δM/M indicate that at small radii the drop in energy is due to the removal of mass. However, larger peaks in δU/U compared to those in δM/M indicate that the gains in thermal energy here are not caused by mass redistribution alone and must also include heating due to the jet action. When jet material is not included, there are differences between the δU/U curves, which are more pronounced for higher power jets at earlier times and at smaller radii, where the energy in the jet material makes a significant contribution to the total energy budget. This fits with the picture drawn in Section 3.1.3, whereby during inflation, the jet lobes retain a significant fraction of energy within the cluster centre. Thereafter, this energy is transferred to the ICM at ever larger radii as the lobes rise and become disrupted by the cluster weather. As for δM/M, the prominent peaks correspond well to the outer regions of the jet feedback influence and any fluctuations beyond this region are small compared to the total thermal energy of the environment.

Differences in thermal energy content between jet and non-jet runs, |$\Delta E_{\rm therm}/E_{\rm inj}^{\rm tot}$|⁠, are estimated in the same manner as the potential energy and are listed in the third column of Table 3. The associated error is largest for the low power jet; however, it is considerably smaller than for the potential energy estimate. At face value, we find that the difference in thermal energy between jet and non-jet runs, as a percentage of jet power, increases with jet power. In other words, more powerful jets are more efficient at heating the ICM. This makes sense physically for a couple of reasons. First, higher power jets drive stronger shocks into the ICM and hence not only convert more of their energy into thermal energy but also heat gas to higher temperatures that has longer cooling times and therefore remains in the system longer. Secondly, short ICM cooling times mean that thermal energy in lobe material mixed into the ICM would radiate more quickly than if it were retained in the lobes. From Fig. 4, we found that higher power jet lobes retain a larger fraction of the injected energy at late times and the longer that energy can be retained in the low-density, high-temperature jet lobes, the more thermal energy the system is likely to retain.

δK/K shows far more fluctuations and although for high jet powers the peaks behave similarly to the fractional δM/M and δU/U, with larger peaks forming at small radii early on before moving to larger radii and flattening with time, this behaviour is less obvious at later times in the fiducial power jet run and is essentially non-existent in the low power jet run. Additionally, fluctuations seen at radii beyond r_out sometimes appear rather large, e.g. the red bump in the low power run just below 50 h⁻¹ kpc. This suggests that differences driven by stochasticity of the simulation (Keller et al. 2019) (for example, the possibility of substructures on slightly different orbits in different runs) potentially have a greater effect on the kinetic energy compared to other quantities.

With this caveat in mind, we present |$\Delta E_{\rm kin}/E_{\rm inj}^{\rm tot}$| at 90.2 Myr in Table 3 but caution that although by this time fluctuations beyond r_out appear small, larger early time fluctuations may now reside within r_out and contribute to our estimates of ΔE_kin. This could explain why the low power jet run seemingly contains a larger kinetic fraction than the fiducial and high power jet runs, which both contain |$\sim \! 8-9{{\ \rm per\ cent}}$| of |$E_{\rm inj}^{\rm tot}$|⁠. On the other hand, if these values are solely driven by jet feedback, a possible physical explanation is that higher power jets result in lower kinetic energy differences (and higher thermal energy differences) because they thermalize a larger fraction of their kinetic energy through shocks and are more effective at inhibiting/thermalizing bulk flows within the ICM (for example, see right-hand panels of Fig. 2, where the ICM velocity fields are more readily disrupted by higher power jet lobes). However, to definitively confirm whether or not this picture is correct, it will be necessary to perform simulations covering a wider range in parameter space in terms of both jet power and cluster environment, which is beyond the scope of present work.

3.2.2 Shocks and sound waves

In this section, we look in detail at the work done by inflating lobes against the ICM, discuss where and when shocks and sound waves are important, and how they compare for differ power jets. We use the shock finder of Schaal & Springel (2015), where cells need to meet three key conditions to constitute the shock zone. Namely, the finder must detect compression (∇ · v < 0), guard against misclassifying weak shocks and contact discontinuities (requiring |$\Delta \log T\geqslant \log [T_2/T_1]|_{\mathcal {M}=\mathcal {M}_{\rm min}}$| and |$\Delta \log P\geqslant \log [ P_2/P_1]|_{\mathcal {M}=\mathcal {M}_{\rm min}}$|⁠), and avoid detecting spurious shocks such as in shear flows (requiring ∇T · ∇ρ > 0).

Fig. 6 shows dissipation-weighted Mach number projections for low-to-high jet powers and at different times. Projections are performed through a distance of one-fifth of the box width. Schaal & Springel (2015) set a minimum Mach number of |$\mathcal {M}_{\rm min}=1.3$| when applying the criteria outlined above, and our results using this setting are shown in the yellow–orange colour maps. Additionally, we attempt to find weaker shocks using |$\mathcal {M}_{\rm min}=1.1$|⁠, as shown by the grey-scale colour maps. The footprint of the jet lobe material is shown by the solid grey contour. At early times, the jets undergo internal shocks and drive bow shocks, inflating lobes of hot plasma that continue to expand into the ICM. The initial rapid expansion of the lobes slows and shocks driven into the ICM weaken and broaden into sounds waves. At late times, strong shocks are again visible at the base of the bottom lobes as the substructure mentioned previously moves through the ICM and interacts with the bottom lobe.

Figure 6.

Dissipation-weighted Mach number projections, showing different jet powers, low to high, in each column from left to right, respectively, and at different times, increasing down each column. The projections are performed through boxes with a depth equal to 1/5th the box width for two different shock finder setups that correspond to the standard finder in yellow/orange and weak shock finder in grey-scale. The figure illustrates that higher power jets typically produce stronger shocks for longer. It is also worth noting that in all runs the weak shock finder detects a feature propagating away from the jet lobes that is consistent with a sound wave driven into the ICM.

Open in new tabDownload slide

As jet power increases, the quantity and strength of shocks increase, with the low power jets showing very few strong shocks at early times only, and high power jets showing strong shocks throughout the jets’ lifetime and even for a short period after the jets switch off. All jet powers show weak shocks, which at later times detach from the lobes and travel at roughly the local sound speed in all directions, suggesting that the shocks have broadened into sound waves. The Mach numbers of shocks driven into the ICM in our simulations are consistent with those observed in real galaxy clusters (e.g. Fabian et al. 2006; Forman et al. 2007; Graham, Fabian & Sanders 2008; Croston et al. 2011; Sanders et al. 2016) and the behaviour found between jet power and shock strength is in agreement with other simulations (e.g. Ehlert et al. 2018).

Shocks can play two distinct roles, namely (i) thermalizing the kinetic jet and inflating the lobes, and (ii) heating the ICM due to shocks driven by the lobe expansion. To understand the first of these processes, we consider shocks that occur within jet material (based on the jet tracer definition), which we call ‘internal’ shocks, versus those that occur in non-jet material. In Fig. 7, we plot dissipation-weighted Mach number projections in grey for the period over which the jet is active. Additionally, we add points to indicate the locations of internal shocks colour coded by their Mach number. This is done for all three jet powers using the results from the |$\mathcal {M}_{\rm min}=1.3$| shock finder and are shown separately in each row as labelled.

Figure 7.

Detection of internal shocks within the jet lobe material. Grey-scale maps show emission-weighted Mach number projections (using |$\mathcal {M}_{\rm min}=1.3$|⁠), while the coloured points show the locations of shocks detected within jet lobe material. Each row shows a different jet power run while each column shows a different time spanning from t = 0 to 19.8 Myr at ∼1.3 Myr intervals. Additionally, the dotted line shows the horizontal intersection of the central black hole. Higher power jets drive more internal shocks that typically have higher Mach numbers.

Open in new tabDownload slide

Strong internal shocks typically occur at earlier times, with only a handful having |$\mathcal {M}\gt 6$| and the majority of shocks are not particularly strong, i.e. many have |$\mathcal {M}\lt 2$| and most have |$\mathcal {M}\lt 3$|⁠. However, the fact that the shocks have relatively low Mach numbers, despite the very high velocities observed in the jets, is because of the high sound speeds within the hot lobes, which in itself is a result of the earlier strong shock heating, meaning that Mach numbers are typically low.

More shocks occur as the jet power is increased. For the low power jets, it seems that internal shocks appear somewhat sporadically and do not extend along the entire jet (in some cases ending well before the bow shock). However, for the fiducial and high power jets the internal shocks are typically seen along the whole jet and particularly in the high power jets reach the region of the bow shocks (although appear confined to the jet axis). Note further that while a large number of shocked cells in the bow shocks are not classified as being within the jet material, they could still be important in providing thermal energy to the jet lobes via mixing (see below).

Fig. 8 shows how the total shock energy produced by the jets evolves. The low, fiducial, and high power jet runs are shown from left to right, respectively, while the top and bottom rows show ‘internal’ shocks and all shocks, respectively. Schaal & Springel (2015) estimate shock energy using Mach numbers defined through the shock temperature jumps. However, due to the possibility of confusing contact discontinuities or weak shocks with strong shocks between the hot lobes and the shocked ambient medium, which exhibit large temperature jumps, we additionally calculate shock energies using pressure jumps. This approach generally results in lower Mach number estimates and hence lower shock energies. Therefore, in Fig. 8 we present shock energies as coloured bands that span between these two sets of estimates. We separate out dissipative heating, which results in the irreversible conversion of kinetic energy to heat and adiabatic compression, as shown by the red and green bands, respectively. The blue band shows the combination of the two, calculated using the minimum Mach value of |$\mathcal {M}_{\rm min}=1.3$|⁠, which due to the reversible nature of adiabatic processes does not reflect the total thermal energy gained by the system but rather the total thermal energy produced via shocks over the lifetime of the simulation. The dark grey band additionally shows the total shock thermal energy budget if we instead assume |$\mathcal {M}_{\rm min}=1.1$|⁠. Finally, we include only cells with f_J > 10⁻¹⁵. Having tested a range of values, we find that this choice suitably limits the detection of shocks driven by processes at large distances and removes the late-time sound wave highlighted in Fig. 6, without compromising the detection of shocks close to the jet lobes.

Figure 8.

Evolution of the shock energy budget compared to the jet energy budget. The top row shows the evolution for internal shocks only, while the bottom row shows the evolution for all shocks in the vicinity of the jet and lobes (i.e. f_J > 10⁻¹⁵). In each case, the energy budget has been split up into the dissipation (red), adiabatic compression (green), and total (blue). An additional band (light grey) shows total energy fractions if we include weak shocks, as calculated by using the shock finder with |$\mathcal {M}_{\rm min} = 1.1$|⁠. The bands are calculated by detecting shocks via pressure (lower bound) or temperature (upper bound) jumps. All panels show the cumulative injected jet energy (black line) and the minimum required energy to inflate the jet lobes (magenta line). Internal shocks make a greater contribution to the total lobe energy as jets become more powerful, but even the high power jet requires gas that has been shocked in the bow shock to contribute to the lobe thermal energy. Additionally, the fraction of shock heating (and the strength of the shocks involved) available to heat the ICM increases with jet power.

Open in new tabDownload slide

Jets initially comprise almost exclusively of kinetic energy; however, from Fig. 4, the vast majority of this energy must rapidly thermalize. The minimum amount of thermalization required is equal to the thermal energy in the lobes minus that of entrained material, plus the energy needed to inflate the lobes. We estimate this as |$E_{\rm lobe}^{\rm req}=\int dH - dE_{\rm ent}^{\mathrm{ therm}}$|⁠, where |$dH=dE_{\rm lobe}^{\rm therm}+PdV$| is the differential lobe enthalpy, and |$dE_{\rm ent}^{\mathrm{ therm}}=\bar{u}_{\rm ICM} dm_{\rm ent}$| is the estimated thermal energy of the entrained ambient material. This estimate is calculated until shortly after 20 Myr, at which point the jets stop and |$E_{\rm lobe}^{\rm req}$| reaches its maximum value and along with the cumulative injected jet energy is shown in Fig. 8 by the magenta and black lines, respectively. As in Fig. 4, all quantities are normalized by the total injected energy for easy comparison between different jet power runs.

The top row shows the total thermal energy injection by internal shocks; we see that while quantitatively they differ, the general behaviour is similar for all jet powers; initially, we measure limited thermal energy injection via shocks, before a rapid increase after a ∼few Myr, followed by a more steady increase in the total energy injection, which eventually plateaus as the jet action halts and no further internal shocks are detected. Additionally, we find that the total shock thermal energy (blue band) is dominated by adiabatic compression as opposed to direct dissipation of kinetic to thermal energy in the shock, which is a result of the majority of the shocks having Mach numbers below 3 (roughly the point at which the contribution due to adiabatic compression is equal to that due to dissipation).

None of the jet powers provides sufficient thermal energy via internal shocks to meet the required conversion of jet kinetic energy into lobe thermal energy and PdV work. If we include weak shocks, the high power jets come closest, with internal shocks accounting for at most |$\sim \! 74{{\ \rm per\ cent}}$| of the required energy, while weak power jets only account for up to |$\sim \! 25{{\ \rm per\ cent}}$| of the required energy through internal shocks. This suggests that shocks detected outside of lobe material must contribute to the thermal content of the lobes via mixing, which becomes increasingly more important for lower power jets. We do not expect turbulent dissipation to play a significant role here as the turbulent energy content of the lobes is only a small fraction of the total injected energy (see Fig. 4).

Therefore, considering the bottom row, which shows the total thermal energy injection for all shocks in the vicinity of the jet lobes, we see that for all three jet powers, there is a significant increase in the amount of energy available. Qualitatively, the energy injection rates show similar behaviours to the internal shocks, although adiabatic compression becomes even more important. The importance of weak shocks (grey band) increases significantly, and thermal energy production can persist beyond 20 Myr as the lobes continue to expand and interact with the ICM. We also note that at ∼20 Myr, shocks driven by the infalling substructure discussed previously intersect the shocks that have travelled furthest from the jet lobes. We attempt to exclude shocks driven by the substructure from our analysis, although this becomes difficult at later times (light grey shaded region in the figure) when the bottom lobe and substructure drive shocks into each other, which, for example likely explains the up-kick in shock energy seen at late times in the fiducial run.

This being said, even without additional shocks due to this substructure, it is clear that the total jet-driven shocks can explain the required conversion of jet kinetic to thermal energy and after the initial ∼5 Myr also provide additional thermal energy to the ICM. Note that the weak shock contribution becomes increasingly more important at low powers and that the total amount of thermal energy generated via shocks increases with jet power. This analysis paints a picture in which a combination of internal shocks and material passing through the jet-driven bow shocks contribute to the thermal energy budget of the jet lobes, with internal shocks becoming increasingly more important as the jet power increases. Additionally, shocks driven into the ICM by the propagation of the jet and inflation of the jet lobes can contribute to thermal energy in the ICM, with weak shocks being more prominent for low jet powers.

The bottom row of Fig. 8 shows that weak shocks are important for injecting thermal energy into the ICM for all jet powers. Without weak shocks, the total shock energy budget at 19.8 Myr is between |$18-27{{\ \rm per\ cent}}$|⁠, |$33-50{{\ \rm per\ cent}}$|⁠, and |$72-102{{\ \rm per\ cent}}$| of |$E_{\rm inj}^{\rm tot}$| for the low, fiducial, and high power jets, respectively. When weak shocks are included, we estimate that the equivalent of between |$63-94{{\ \rm per\ cent}}$|⁠, |$92-131{{\ \rm per\ cent}}$|⁠, and |$122-168{{\ \rm per\ cent}}$| of |$E_{\rm inj}^{\rm tot}$| contributes to shock thermal energy production in total for the low, fiducial, and high power jet runs, respectively. While we could expect these to exceed |$100{{\ \rm per\ cent}}$|⁠, given that the same energy can be recycled via multiple shocks resulting in double counting of the same energy,⁵ the upper values, particularly for the high power jets, seem excessive. Even if we assume an extreme case in which all of the jet kinetic energy is thermalized through shocks and the resulting PdV work done by the inflation of the lobes on the ICM also goes into thermal energy, this sets an upper limit of |$\sim \! 145{{\ \rm per\ cent}}$|⁠. In reality, the ‘detached’ sound waves seen in Fig. 6 would themselves contain a fraction of this energy, and so the true upper limit for pure thermal energy injection via shocks would be somewhat smaller than this. As such, measuring shock thermal energy budgets similar to or exceeding this value would indicate that we are detecting spurious weak shocks/sound waves, capturing shocks from processes other than the jet inflation (e.g. the infalling substructure), and/or preexisting ICM kinetic energy is additionally being thermalized via shocks at the lobe/ICM interface. The truth is likely to be a combination of these, and as such the values we quote here should be taken as strict upper limits on the shock energy budget.

Taking into account the thermal energy required from shocks to inflate and thermalize the lobes, |$E_{\rm lobe}^{\rm req}$|⁠, and subtracting it from the total shock thermal energy budget, we find that when weak shocks are included, by 19.8 Myr, the equivalent of between |$0-21{{\ \rm per\ cent}}$|⁠, |$15-54{{\ \rm per\ cent}}$|⁠, and |$43-90{{\ \rm per\ cent}}$| of |$E_{\rm inj}^{\rm tot}$| for the low, fiducial, and high power jet runs, respectively, is potentially available to be communicated to the ICM.⁶ We stress, however, that for the reasons discussed above, these values are upper limits, which especially for the high power jets are likely to overestimate the amount of thermal energy injected into the ICM via jet-driven shock. However, a clear picture still emerges showing that higher power jets can contribute more thermal energy to the ICM via shocks as a fraction of their total energy, and as already discussed above, this is facilitated via stronger (high |$\mathcal {M}$|⁠) shocks. Combining these results provides validation for the picture presented in Section 3.2 in which systems with more powerful jets can retain more thermal energy (as a fraction of injected energy), and conversely, the energy communicated to the ICM via lobe expansion of weaker jets more likely leads to relatively gentle displacement of the ambient medium rather than violently shock heating it.

3.2.3 Mixing

At late times, mixing has the potential to become a dominant energy transfer mechanism, where it can be comparable to or even exceed our PdV estimates. However, there are some caveats when studying mixing in numerical simulations and trying to untangle physical and numerical effects. While Eulerian codes inherently generate entropy through mixing, they tend to suffer from numerical over-mixing (Mitchell et al. 2009; Springel 2010; Hopkins 2015). On the other hand, an inability to accurately capture fluid instabilities means that mixing is inefficient/suppressed in SPH codes (Dolag et al. 2005; Agertz et al. 2007; Springel 2010; Sijacki et al. 2012) unless additional numerical prescriptions are included (e.g. Price 2008; Wadsley, Veeravalli & Couchman 2008; Read, Hayfield & Agertz 2010; Read & Hayfield 2012). However, the moving mesh nature of arepo means that it can capture mixing without suffering from significant over-mixing (Springel 2010; Sijacki et al. 2012). Additionally, through the use of appropriate refinement criteria over-mixing of jet lobes in idealized cluster environments through shredding (e.g. Reynolds et al. 2005) is also avoided when appropriate resolution is afforded to jet material (Bourne & Sijacki 2017; Weinberger et al. 2017; Ehlert et al. 2018).

This being said, an implicit assumption within our simulations is that the fluid within a resolution element is mixed instantaneously, which although results in entropy generation when fluids of different specific entropies combine, neglects physical processes driving mixing below the resolution scale, and could potentially result in shorter time-scales than expected physically (with obvious extra levels of complexity arising when considering mixing of relativistic and thermal plasmas in the presence of magnetic fields). For example, while turbulence and fluid instabilities (e.g. Kelvin–Helmholtz and Rayleigh–Taylor) are captured within our simulations, the physical scale on which the instabilities and turbulence that drive mixing occur is limited by our resolution. Interestingly, however, our previous work found that although improved resolution resulted in reduced mixing between lobes and the ICM, the trend was converging and overall differences in the lobe masses and energy content were small, thus suggesting that at high enough resolution further improvements lead to not very significant changes in the mixing rates (Bourne & Sijacki 2017).

Bearing this in mind, in an attempt to illustrate the relative importance of mixing for different jet powers, in Fig. 9 we plot the volume probability distribution as a function of jet tracer normalized by the total volume, dV/(Vdlog₁₀f_J), i.e. the fraction of the total volume in the distribution with a jet fraction between log₁₀f_J and log₁₀f_J + dlog₁₀f_J, from low to high jet power, left to right, respectively. The distributions are shown at different times by different coloured curves (see legend) and the regions of jet tracer space occupied by the lobe material are shown by the grey shaded regions bounded by the solid, vertical black line, while vertical dotted lines indicate the tracer threshold values used for the other two jet powers.

Figure 9.

Evolution of the jet tracer volume probability distributions with time as a proxy of mixing efficiency. From left to right, the panels show the volume probability distribution as a function of f_J (and normalized by total volume) for the low, fiducial, and high power jet runs, respectively. Each panel contains distributions at five different times and the values of |$f_{\rm J}^{\rm thresh}$| are indicated by the vertical lines (solid for the corresponding jet power, dotted for other jet powers). The figure shows that mixing is more effective for lower power jets, where the peak in the distribution moves further to the left (lower jet tracer fraction) and sooner, when compared to higher power jets.

Open in new tabDownload slide

As jet material mixes with the ambient medium, it is diluted, reducing its value of f_J. With time, the amount of jet material mixing to lower values of f_J increases, and a larger fraction of the gas volume is occupied by more dilute jet material driving the peak of the distribution to the left. This behaviour is seen for all three jet powers, with the distributions for the three earliest times remaining within the jet lobe boundary, albeit with the peaks still moving to the left. However, after ∼30 Myr, differences begin to appear; the peak has shifted well to the left for the low power jet by 49.4 Myr, while the same can be said of the fiducial jet power only at 90.2 Myr. On the other hand, at 90.2 Myr, the peak in the high jet power distribution is still on the boundary between jet and non-jet material. Putting aside the differences in jet tracer thresholds used for each jet power run, the late-time peaks occur at much smaller values of f_J for lower power jets, indicating that the low power jets are more readily mixed. This is in agreement with the findings in Section 3.1.3. While we expect that the lower power jet lobes are more susceptible to the cluster weather and therefore are more readily mixed, we should also highlight that similar conclusions have been found previously in more idealized environments (e.g. Ehlert et al. 2018).

3.3 Impact on cluster properties

3.3.1 Overview and radial profiles

To give an overview of how the jets impact the cluster gas at early and late times, Fig. 10 shows 2D radial histograms of the gas entropy (grey) at 19.8 and 90.2 Myr in the top and bottom rows, and for the non-jet, low, fiducial, and high power jets from left to right, respectively. The red and blue colours in the jet run plots show positive and negative residuals when the non-jet profile is subtracted from the corresponding jet-run profiles. Jet lobe material is shown by the green contour lines. The high-entropy jet lobes are visible in the top row, with higher entropy gas existing in the lobes of higher power jets, which heat the gas to higher temperatures and inflate larger, lower density lobes. Considering the residuals in the top row, it is also clear that jets remove material irrespective of entropy at smaller radii and also remove low-entropy material at intermediate radii. At later times, high-energy jet lobes reach larger radii but contain typically lower entropy gas compared to younger lobes, while the low power jet lobes appear as only a small bump above the background entropy level. The impact of the jets on the ambient medium can still be seen at late times as coherent regions of red and blue, where the radial entropy distribution of material is different to that of the run with no jet, typically showing that the jets remove mass from low radii and low-to-intermediate entropy regions of phase space. All jet powers cause noticeable increases in the ICM entropy compared to the run with no jet feedback, with the highest power jet resulting in longer lived features that extend to larger radii.

Figure 10.

Two dimensional radial entropy histograms, shown for jets of increasing power from left to right and at t = 19.8 Myr and 90.2 Myr in the top and bottom rows, respectively. The total mass-weighted histograms are shown in the grey-scale, while residuals from subtracting the non-jet profile from each jet profile are shown in red (+ve) and blue (−ve). Additionally, the jet lobe content is shown by contours of jet material mass in green. The figure shows that the jets inflate lobes of high entropy gas, which move out through the cluster atmosphere and modify the entropy of ambient medium by generally increasing it.

Open in new tabDownload slide

In using a self-consistent cosmological simulation including subgrid models for gas radiative cooling and heating, star formation, evolution, enrichment, and stellar feedback, our resultant galaxy cluster contains a plethora of features not found in idealized setups. Our simulations capture the cluster weather driven by processes such as mergers and cosmic accretion. Additionally, the simulations capture a range of gas phases that are relevant to galaxy and cluster evolution. In what follows it is useful to consider the impact that the jet feedback has on each of these different phases and how gas transitions between them. We define four key phases (see a phase diagram in Appendix  A) whereby the gas is split into jet material, which we have already defined as non-ISM gas with |$f_{\rm J}\gt f_{\rm J}^{\rm thresh}$|⁠, ISM, which is any gas with densities above ρ_crit = 8.5 × 10⁵ h² M_⊙ kpc⁻³ and is considered star-forming, the ICM, which is defined using the temperature density cut, |$T_{\rm keV} \gt 3\times ~10^6\rho ^{0.25}_{\rm cgs}$|⁠, of Rasia et al. (2012), and finally the warm phase, which sits between the ICM and ISM phases and is defined as |$T_{\rm keV} \lt 3\times ~10^6\rho ^{0.25}_{\rm cgs}$| and ρ < ρ_crit.

To provide a general overview of the impact of the jet on different gas phases, cluster radial profiles are shown in Fig. 11. The figure is arranged as follows: the top, middle, and bottom rows show the radial density (mass-weighted), temperature (volume-weighted), and entropy profiles, respectively, while from left to right, the columns show profiles calculated at different times during the jet evolution. The different jet powers are shown in different colours (as shown in the legend at the bottom of the figure) while different line styles indicate different gas-phase combinations. Specifically, the solid lines show pure ICM gas (ICM profile), the dashed lines show the combined ICM and warm-phase profiles (Warm–Hot profile, WH), while the dotted lines show the combined profile of the ISM, warm and ICM phase (Cold–Warm–Hot profile, CWH), noting that the jet material is excluded from all profiles. Additionally, in the first column, we show the initial cluster profile (before jet switches on) with the black lines. Before discussing these profiles, we want to stress two key points: first, the initial cluster profiles are the result of subgrid models typical of cosmological simulations and not a reflection of the new jet feedback model we present here, and secondly, the jet feedback we implement assumes a fixed jet power and lifetime and is in no way coupled to the cluster properties, as such these profiles are not presented as strict comparisons to observations but rather as an analysis tool for understanding the cluster heating process.

Figure 11.

Radial profiles for density (top), temperature (middle), and entropy (bottom) at different stages of the jet evolution, moving forward in time from left to right. Different line colours correspond to different jet powers, while the different line styles correspond to different gas phase combinations, namely ICM (solid), ICM + Warm (dashed), and ICM + Warm + ISM (dotted). Note that jet material is neglected in all profiles. Although differences can be seen in the profiles due to the jet action, these are largely confined to earlier times and for higher power jets. The ICM profile varies very little between different jet powers, while all profiles look very similar by 90.2 Myr.

Open in new tabDownload slide

First, a brief look at all of the profiles shows that beyond the central ∼30 h⁻¹ kpc, the differences in profiles become negligible. This is to be expected as the total injected energy is equal to the initial ICM thermal energy content within radii of ∼10 h⁻¹ kpc, 18 h⁻¹ kpc, and 30 h⁻¹ kpc, for the low, fiducial, and high power jet runs, respectively. This means that any energy injected into the ICM at larger distances very quickly becomes small compared to the total ICM energy contained within this region. It is also worth pointing out that at these distances, the profiles are dominated by the ICM. Considering the profiles just before the jets turn off, at 19.8 Myr, jets of all powers have transferred |$\sim \! 50{{\ \rm per\ cent}}$| of the jet energy to the local environment and all present differences in the profiles compared to the run with no jet. Specifically, the jets remove gas, reduce the density, and heat all phases within the central 1–2 kpc region, which becomes dominated by jet material (not shown on this plot).

Beyond this central region, the CWH and WH profiles show large density enhancements as well as decreases in temperature and entropy when compared to the non-jet run, especially in the case of the fiducial and high power jets, while the ICM profiles show much weaker reactions to the jet action with mild density enhancements and temperature/entropy decreases for the fiducial and high power jets only, with the ICM in the low power jet run showing almost no difference when compared to the non-jet run. The differences between the jet/non-jet run CWH and WH profiles can be driven by two processes: (i) the compression of gas as the inflation of the jet lobes pushes it out to larger radii and/or (ii) the removal of hot ICM and/or warm phase gas. The former process by construction increases the density and temperature of the gas (although ultimately the increased density results in shorter radiative cooling times, with gas temperatures subsequently falling), while the latter process means that the cooler phases make a more significant contribution to average profile quantities giving the appearance of gas densities/temperatures increasing/decreasing. As we show below in Section 3.3.2, while the jet does compress all gas phases during the inflation (although this peaks before 20 Myr), a significant amount of hot ICM gas is also cleared out in the fiducial and high jet power runs, and so it is a combination of these two processes that drive the differences in the radial profiles.

Considering only the CWH and WH profiles, at later times, after the jet action has halted, while the high power jet run still exhibits differences, the lower power jet profiles become increasingly more similar to those of the non-jet run. Specifically, at 43.3 Myr, all jet powers still show significant differences within the central ∼2 kpc in the CWH profiles, although beyond this both the low and fiducial jet power runs are very similar to the non-jet run and only the high power jet shows large differences for both the CWH and WH profiles. By 68 Myr, the differences are even smaller, the low and fiducial power jet run profiles show only minor differences to the non-jet run profiles for all radii, while the high power jet profiles show only a significant difference in the CWH profiles. Finally, by 90.2 Myr even the high power jet run shows only very small differences compared to the non-jet run.

Importantly, if we consider the ICM entropy profiles, we see that the low power jet results in a mean entropy excess in the cluster core even at 43.3 Myr, while an excess still exists in the 68-Myr profile for the fiducial jet power and the high power jet run exhibits a not insignificant bump in the entropy profile out to ∼20 kpc at 90.2 Myr when compared to the run without a jet. The jet feedback increases the mean ICM entropy, which in the case of the high power jet persists long after the jet has switched off. Comparing to Fig. 5, we find that the mean entropy excess appears at smaller radii than the peaks in δU/U. However, this can be understood by noting that although on average higher entropy gas is seen at smaller radii in Fig. 11, there is also a mass deficit at small radii (but a mass excess at larger radii) that, at least in part, explains why the δU/U peak does not match the location of the higher mean entropy.

However, if we were to take each radial profile in isolation, the system evolves such that it relaxes back to a state that exhibits very few ‘scars’ of the jet feedback. As we discuss further in Section 3.3.2, in the cluster core, the ISM gas in the jet runs relaxes somewhat such that the mean density and temperature become closer to those measured in the non-jet run, while ICM material that was initially expelled by the jet action begins to repopulate the cluster core, which combined contributes to the radial profiles becoming more similar at later times. Additionally, within the ICM and particularly at larger radii where it dominates the mass and energy budget, it would require a significant energy injection to make a long-lived and noticeable difference to the ambient density and temperature. This implies that in reality, the jet feedback is quite subtle, meaning that even though the cluster retains a significant fraction of the injected jet energy (see Fig. 5), its impact on the ICM radial profiles become less apparent with time as the energy is distributed over a larger volume and converted into potential energy.

3.3.2 How different phases in the cluster core evolve

We explore the impact of the jet feedback in the cluster centre (r ≤ 20 h⁻¹ kpc) in more detail in Fig. 12. On the left-hand side, we plot the evolution of the ISM (left column), warm phase (middle column), ICM (right column), total mass (top row), mean density (middle row), and mean temperature (bottom row) normalized by either the initial value at t = 0 (dashed) or the value in the non-jet run at the same time (solid curves). The non-jet, low, fiducial, and high jet power runs are shown by the black, blue, red, and green lines, respectively. Additionally, on the right-hand side, we show column density projections of the central region for all four runs increasing in power from no jets to high power jets from top to bottom, respectively, and evolving in time from left to right. Due to the nature of our simulations, which includes processes other than those solely driven by AGN feedback, even the non-jet run (black lines) shows an evolution in all of the properties considered; however, we see that the jet feedback introduces differences in how these properties evolve.

Figure 12.

Evolution of different gas phases within the centre of the cluster. On the left-hand side, we show quantitatively how different gas phases: ISM (left), warm (middle), and ICM (right) evolve with time within the central 20 h⁻¹ kpc of the cluster. We consider how the total mass (top), mean density (middle), and mean temperature (bottom) evolve with respect to the IC (dashed) and the non-jet run (solid curves). On the right-hand side, we show column density projections (through boxes of 40 h⁻¹ kpc on a side) for each of the four runs, increasing with power from top to bottom and evolving in time from left to right, respectively. The figure illustrates how different power jets impact the phases early on, but at later times, the system is able to some extent relax back to similar states, most noticeably for the cold phase and least for the ICM.

Open in new tabDownload slide

First considering the ISM, which dominates the central gas mass budget (⁠|$79{{\ \rm per\ cent}}$|⁠), we note that the structure evolves significantly even in the non-jet run, highlighting that processes including star formation and associated feedback as well as dynamical interactions due to passing substructures can have an impact on the cold gas properties. At early times, the central ISM mass increases primarily due to a substructure, which can be seen to the left of the BCG at 19.8 Myr in the column density projections, entering the region. While the ISM mass remains similar for all simulations, it is clear that the distribution of this material is different. In particular, the jet simultaneously excavates a hole in the centre of and expels material from the cold disc,⁷ with the effect being more extreme for more powerful jets. There are also differences in the mean density and temperature of the ISM. Lobe inflation compresses gas in the central disc, increasing its density and hence shifting warm gas into the ISM phase as well as ‘pushing’ ISM gas along the effective EoS. This results in early peaks in the mean ISM temperature, especially for the fiducial and high power jets, which subsequently drop as the gas cools. This means that during inflation, while the jet does not destroy significant amounts of cold gas, it can redistribute it and change its properties. As we discuss later in Section 3.3.3, the increase in density impacts the SFRs found in the BCG.

After ∼14 Myr, the total mass begins to slowly decline in all runs as a substructure begins to leave the region, and the decline steepens after ∼30–35 Myr as the object finally leaves the central region. At around the same time, small differences in the ISM mass between the runs start to appear and while by ∼40 Myr, the central densities in the three jet runs become similar in value and evolution and they appear to be systematically offset to higher values and peak earlier than the non-jet run. By this time, other substructures have left the central region and differences in the total ISM mass for all four runs become more obvious, with cold material expelled by the jet continuing to leave the central region. However, by late times, we see that while there are small differences in the total ISM mass, the mean ICM density and temperatures are similar for all jet and non-jet simulations and looking at the column density projections, we see that despite the perturbations caused by the jet, the cold gas appears to have a very similar structure, which hides any jet-induced ‘scars’.

The warm phase, which contributes the least to the total gas mass (⁠|$\sim \! 5{{\ \rm per\ cent}}$|⁠), exhibits a continuous drop in mass throughout all four simulations, as it is either heated and joins the ICM or cools on to the ISM. While these processes occur in the non-jet run, the action of the jet seems to enhance them, with jet simulations retaining less warm phase material, as it is either heated or compressed by the jet and transitions to one of the other phases. The density and temperature of the warm phase fluctuate over time, but as in the case of the ISM, there is a clear compression during the jet feedback phase seen as a density and temperature increase, with the increase being larger for higher power jets.

Finally, considering the ICM, in general, the non-jet run exhibits an increase in mass in the central region up until ∼70 Myr, at which point it begins to decline. In terms of the total mass content, which comprises |$16{{\ \rm per\ cent}}$| of the central gas mass initially, we can see that different jet powers can have very different effects. While the low power jet results in a small reduction in the mass of ICM compared to the non-jet run, the fiducial jet results in a decrease in the total ICM material for early times, while the high power jet produces a long-term reduction in the amount of ICM in the cluster centre. This is a direct result of the jets and the lobes they inflate, displacing and clearing out the hot ICM material that once resided at the location of the lobes and can be seen in the early column density projections that show a reduction in ICM material in the central region surrounding the disc. The action of ‘sweeping up’ the ICM compresses it, resulting in clearly defined density and temperature peaks seen during the first ∼20 Myr, as well as the production of the X-ray bright rims shown in Fig. 1. We additionally note that this compression results in shorter gas cooling times within the higher density gas.

However, once the compressed material leaves the central region, there is a drop in the mean density and temperature, which continues for a short time after the jet switches off, with differences compared to the non-jet run, especially for the high power jet, continuing to exist. We also note that shortly after the jet action stops, the ICM starts to replenish with the differences in mass content between the runs becoming smaller with time. Additionally, the increased ICM cooling rates observed during the lobe inflation decrease after the jet switches off. We find that the changes to ICM properties induced by the jet action lead to ICM cooling rates becoming lower than those in the non-jet run, even out to ∼100 h⁻¹ kpc in the case of the high power jet, which will hinder future gas cooling inflows. Overall, this figure paints a picture in which the general evolution of the central ISM is largely determined by processes other than the jet feedback, at least on the time-scale of our simulations. However, the jet can moderate – either increase or reduce – variations in the properties of all phases to varying degrees and, particularly in the case of the higher power jets, reduce the ICM and warm phase gas content. This has the potential to affect the reservoir of hot gas that could cool in the future and can also have additional implications that result from the non-linear evolution of the system, for example through the impact of AGN feedback on star formation, as we discuss next.

3.3.3 Star formation

As shown in Section 3.3.2, the jet can compress gas to higher densities. The resulting increased density in the ISM will have an impact on the SFR of the central BCG, which we now explore. Fig. 13 is split into four panels that from top to bottom show the SFR, the fractional difference in SFRs between jet and non-jet runs, the integrated stellar mass, and the fractional difference in integrated stellar mass between jet and non-jet runs, respectively. The different simulation runs are shown by different colours, as indicated by the legend, while the solid and dashed lines show values within 10 h⁻¹ kpc and 20 h⁻¹ kpc, respectively.

Figure 13.

Star formation within the centre of the cluster, shown for r < 10 h⁻¹ kpc (solid) and r < 20 h⁻¹ kpc (dashed). From top to bottom, respectively, the panels show the total SFR within the corresponding regions, the fractional difference in SFR between jet and non-jet runs, and the integrated SFR and the fractional difference in the integrated SFR between jet and non-jet runs. During the inflation phase, jets can increase the central SFR through the compression of gas, but at later times, the SFRs subside and the net difference in the integrated stellar mass is similar between all runs.

Open in new tabDownload slide

In line with the increase in density of the cold ISM seen in all runs during the first ∼10 Myr, the SFR also increases. While the low power jet provides only a small increase with respect to the non-jet run, the difference becomes larger with increasing jet power, with the high power jets resulting in a |$\sim \! 25{{\ \rm per\ cent}}$| increase over the non-jet run. This is due to the higher power jets being able to compress gas to higher densities, which results in higher SFRs, with similar findings found by some previous works (e.g. Zubovas & Bourne 2017; Mukherjee et al. 2018; Valentini et al. 2020). Despite the increased SFRs during the jet outburst, the total integrated stellar masses after 100 Myr are similar between all runs. This is because even though the SFRs remain elevated in the jet runs even after the jets deactivate, eventually they drop below the SFR in the non-jet run, firstly for the high power jets, followed later by the fiducial and then the low power jets. Differences in SFR will depend on both the amount of cold gas available and the properties of that gas. Considering the left-hand column in Fig. 12, we see that while there are corresponding peaks in SFR and mean ISM density, the size of the peaks is not correlated, i.e. the strong early SFR peak does not correspond to the highest density peak. As such, we expect that the late-time differences in SFR between different runs are governed primarily by the amount of cold ISM available and as such that the behaviour, which drives the differences in the mass of ISM in the central region (top left-hand panel, Fig. 12), is also responsible for the decrease in SFR.

4 DISCUSSION

The art of simulating galaxy cluster populations has come into its own in recent years, with a number of groups performing suites of cosmological-zoom simulations that include a range of different cluster masses (Le Brun et al. 2014; Rasia et al. 2015; Bahé et al. 2017; Barnes et al. 2017a, b; Hahn et al. 2017; McCarthy et al. 2017; Barnes et al. 2018; Henden et al. 2018; Henden, Puchwein & Sijacki 2019a, b). These works have made great progress in being able to reproduce many observed galaxy cluster properties; however, matching thermodynamic profiles, particularly in the cluster cores, and explaining the cool-core/non-cool-core (CC/NCC) dichotomy remain a challenge. For example, overly effective AGN feedback in the Illustris simulations resulted in too much gas being ejected from groups and low-mass clusters (Genel et al. 2014). This issue was resolved for Illustris-TNG, although the CC fraction and its evolution do not match observations due to differences between simulated and observed gas fraction profiles (Barnes et al. 2018), suggesting that the AGN feedback model needs modification. Another symptom of overly efficient AGN feedback was evident in both C-Eagle (Barnes et al. 2017b) and FABLE (Henden et al. 2018), which while being able to match multiple observables found overheating in the cluster centres. Finally, while Rhapsody-g simulations could reproduce the CC/NCC dichotomy, it was at the expense of matching observed X-ray luminosities (Hahn et al. 2017).

The work we have presented here attempts to improve our theoretical understanding of AGN jet feedback by bridging the gap between jet injection and galaxy cluster scales. The focus hence shifts to model observational signatures and behaviour of the jet itself as opposed to its impact on the cluster environment, which we hope to arise naturally from a more physically motivated feedback model. It is important to highlight that while several cosmological works invoke jet mode feedback, this typically mimics the effect of the jet feedback, for example by adding hot bubbles by hand (e.g. Sijacki & Springel 2006a; Sijacki et al. 2007, 2015), or by injecting bipolar kinetic outflows (e.g. Dubois et al. 2010, 2012; Davé et al. 2019). The models are typically (and by necessity) simplified and not designed to track the lobe inflation and their subsequent evolution in detail.

There are a limited number of previous studies that similar to this work attempt to capture the lobe inflation and evolution in a cosmological environment (e.g. Heinz et al. 2005; Morsony et al. 2010; Mendygral et al. 2012). Here, we expand upon these studies by including additional physical models (such as radiative cooling, star formation, and associated feedback) by forcing very high resolution in the jet injection region and crucially by adopting a suite of refinement techniques that ensure that the high resolution (∼120 h⁻¹ pc) tracks the jet lobes at all times to capture the details of the lobe inflation process and its subsequent interaction with the ICM. Importantly, we find in common with these previous studies (see also Sijacki et al. 2008; Bourne & Sijacki 2017) that cluster weather, due to the preexisting turbulence, bulk motions, and orbiting substructures, can play a very significant role in redistributing and/or disrupting jet lobes.

In idealized simulations of jet feedback in hydrostatic atmospheres, jet lobes inflate and rise along the jet axis. Such simulations result in low-density channels along one axis through which feedback energy can easily escape and hence couple poorly with the ICM (Vernaleo & Reynolds 2006). Previously, this issue has been overcome in simulations by including precession (Falceta-Gonçalves et al. 2010; Li & Bryan 2014; Yang & Reynolds 2016a, b) or rapid re-orientation of the jet axis by hand (Cielo et al. 2018) to isotropize and increase the coupling between the jet energy and ICM. However, when cluster weather is included, as found here, the jet lobes can be naturally displaced from their original trajectory by bulk motions in the ICM (see also, Sijacki et al. 2008; Morsony et al. 2010; Bourne & Sijacki 2017). Given the expected difficulty in driving precession of the jets from massive BHs by accretion (Nixon & King 2013), cluster weather provides a possible explanation for observations of successive generations of jet lobes being misaligned in some galaxy clusters (Dunn et al. 2006a; Babul et al. 2013).

Jet lobes in systems with cluster weather are effectively mixed with the ICM (see also, Heinz et al. 2005; Morsony et al. 2010; Bourne & Sijacki 2017). We note that some previous works have suggested that mixing of hot lobe material with the ICM through jet-driven instabilities is an effective method of heating the ICM (Hillel & Soker 2016, 2017), even in a hydrostatic environment. However, we stress that here we are emphasizing the importance of cluster weather-driven mixing, noting that at sufficiently high resolution and by employing appropriate refinement criteria, in Bourne & Sijacki (2017), we found instability-driven mixing to be largely unimportant in hydrostatic environments, even in the absence of magnetic fields (see also, Weinberger et al. 2017; Ehlert et al. 2019). We further note that while our simulations are performed without magnetic fields, which have been shown to affect lobe dynamics, and inhibit small-scale instabilities along the lobe-ICM interface that facilitates mixing (Dursi & Pfrommer 2008; English et al. 2016; Weinberger et al. 2017), we expect that the drastic nature of the cluster weather on lobes (particularly the bottom lobes) would still have a strong impact and drive mixing even when magnetic fields are included.

An important (and related) point to consider is the potential relationship between the dynamical state of the ICM and whether or not a cluster is a CC or NCC. Observations, such as the REXCESS sample, show that there is a lack of CCs in dynamically disturbed systems (Pratt et al. 2010), while Mahdavi et al. (2013) found that while NCC clusters exhibit an X-ray mass bias compared to weak lensing mass measurements, this is not the case for CC clusters, suggesting that CCs are dynamically relaxed. If, as we have found here, cluster weather can disrupt jet lobes and aid in the process of heating the ICM, then we could expect X-ray cavities to be more common in CC clusters, as is the case from observations (Burns 1990; Dunn et al. 2005; Dunn & Fabian 2006, 2008; McNamara & Nulsen 2007; Sun 2009; Fabian 2012). The common occurrence of lobes in CC clusters has often been attributed to a link between cooling, BH feeding, and subsequent feedback; however, there may also be an observational bias in that lobes are more likely to survive and hence are preferentially observed in dynamically quiet environments.

Additionally, the potentially lower efficacy of feedback in CCs may limit its role in disrupting and converting them to NCCs. Even in the simulations presented here in which cluster weather helps liberate feedback energy, it does not have a significant impact on the cold gas content of our cluster. However, the feedback may instead be able to aid in preventing CCs forming in the first place. For example, McCarthy et al. (2008) found that early AGN heating can indeed inhibit the formation of CCs. Alternatively, numerous simulations have found that it is mergers, and not AGN feedback, that disrupt CCs (e.g. Burns et al. 2008; Poole et al. 2008; Rasia et al. 2015; Hahn et al. 2017; Chadayammuri et al. 2020), with Hahn et al. (2017) and Chadayammuri et al. (2020) finding that, similar to the results presented here, AGN feedback is ineffective at significantly altering/removing the cold gas content of the cluster core. This being said, it is also worth noting that Barnes et al. (2018) found from the Illustris TNG cosmological simulations that there is no difference in the dynamical state of CCs versus NCCs, which would instead suggest that mergers play a limited role in converting CCs to NCCs. These differing conclusions illustrate the need for future simulations to combine both realistic cluster environments with physically motivated, high-resolution AGN feedback models.

The simulations presented here consider a single cluster at a single epoch. This has allowed us to garner insights into how AGN feedback can heat the ICM and impact different gas phases within objects of this kind. However, for a full understanding of the role of AGN feedback in shaping galaxy clusters, it will be necessary to consider clusters with a wide range of properties over the whole of cosmic history to understand how the ICM forms and evolves, the development of the entropy excess (Voit & Donahue 2005; Sun et al. 2009), turbulence/cluster weather (Pinto et al. 2015; Hitomi Collaboration 2016; Bourne & Sijacki 2017; Vazza et al. 2017), and magnetic fields, and how each of these impacts the AGN feedback process in both CC and NCC systems (see discussion above). This will also require the inclusion of a self-consistent accretion and feedback model that has not so far been included in our simulations to determine both the jet direction and the production efficiency and close the feedback loop between large-scale gas inflows and central AGN heating.

Additionally, as with many large-scale cosmological simulations, it is not possible to include all physical processes due to a combination of the large dynamic range of the systems considered and current computational or physical limitations. Specifically, we have neglected MHD effects, which have been shown to inhibit instability-driven mixing of jet lobes in idealized simulations and impact lobe dynamics (e.g. Dursi & Pfrommer 2008; English et al. 2016; Weinberger et al. 2017; Ehlert et al. 2018), although these effects are seemingly negated in a large part by implementing sufficient resolution and refinement criteria (Bourne & Sijacki 2017; Ehlert et al. 2018). However, given the growing number of cosmological simulations that successfully include the effects of magnetic fields (e.g. Dolag, Bartelmann & Lesch 2002; Marinacci et al. 2015, 2018; Martin-Alvarez et al. 2018; Katz et al. 2019), it is timely to run a comparison suite in future works, in particular, to check if the cluster weather-aided mixing is as prominent when magnetic fields are included. Neglecting magnetic fields additionally means that thermal conduction is not included in these runs, which, while adding additional levels of complexity to simulations, has been shown to play a potentially important role in solving the cooling flow problem (e.g. Voigt & Fabian 2004; Bogdanović et al. 2009; Ruszkowski & Oh 2010, 2011; Yang & Reynolds 2016a; Kannan et al. 2017).

A further simplification we adopted is in modelling the hot lobes as a non-relativistic thermal gas (γ = 5/3) in which the electrons and protons share a single temperature. The exact nature of the lobe plasma is still uncertain (e.g. Mazzotta et al. 2002; Schmidt, Fabian & Sanders 2002; Dunn et al. 2005; Dunn, Fabian & Celotti 2006b; Bîrzan et al. 2008; Croston et al. 2008; Croston & Hardcastle 2014; Kang, Chen & Wu 2014; Kawakatu, Kino & Takahara 2016), and in reality, the two components may not be in thermal equilibrium and have different temperatures and their distributions may be non-thermal. Uncertainties surrounding how CR energy is communicated to the ICM via streaming and diffusion mean that it is unclear what fraction of the CR energy is actually transferred to the ICM and on what time-scales. Given that we do not include any models for CRs with our simulations, the mixing rates we present may not map directly to those expected in real clusters but rather provide a qualitative mechanism by which cluster weather could enhance mixing of lobe material. In fact, some works have studied the effect of different lobe content, for example by including the CR component (e.g. Sijacki et al. 2008; Weinberger et al. 2017; Ehlert et al. 2019; Yang, Gaspari & Marlow 2019), finding that CRs can impact lobe dynamics, with the additional CR pressure resulting in more buoyancy and easier uplift of ICM gas, and produce ‘rounder’ bubbles compared to kinetically injected jets that are more similar to observed X-ray cavities, suggesting that CRs could have a role to play in ICM evolution.

Finally, we do not include a physical viscosity model and hence any viscous dissipation occurring is purely numerical. As discussed in Section 3.2.2, sound waves could contribute to heating the ICM (e.g. Fabian et al. 2005; Sanders & Fabian 2007; Fabian et al. 2017; Bambic & Reynolds 2019) and, as such, a full understanding of where this occurs and at what rate would depend on the physical viscosity of the ICM (Ruszkowski et al. 2004; Sijacki & Springel 2006b; Zweibel et al. 2018), which itself is still rather uncertain. However, given that arepo has very low numerical viscosity, sound waves are expected to be long lived. As such, the fact that the weak shocks/sound waves in the low power jet run so readily disappear suggests that the cluster weather can at least partially disrupt them, as was also found in Bourne & Sijacki (2017). Therefore, transporting energy over a large volume via sound waves may only be viable in relatively quiet cluster atmospheres or with sufficiently high jet powers.

5 SUMMARY

We have presented results from simulations of very high-resolution jets launched into a fully cosmologically evolved galaxy cluster. The work builds on that of Bourne et al. (2019) by extending the range of jet powers studied and performing in-depth analysis to investigate lobe inflation and dynamics, as well as the jet-lobe energetics and how the jet energy is communicated to and impacts the ICM. A simple picture arises in which rather gentle feedback is able to heat the ICM via a combination of weak shocks and sound waves during lobe inflation, while cluster weather subsequently plays a critical role to effectively liberate lobe energy and transfer it to the ICM. The exact details depend on jet power and our main findings are as follows:

• With appropriate refinement criteria that provide sufficient resolution in the injection region as well as crucially in the jet lobes, it is possible to inflate low-density, high-temperature lobes within a full cosmological zoom simulation of a galaxy cluster. The lobe properties result in mock X-ray and radio appearances that compare favourably to observations of jet lobes in galaxy clusters.

• Analysis of lobe energetics shows that irrespective of the jet power, the lobes retain approximately half of the injected energy during the lobe inflation phase, with the majority of the ‘lost’ energy going into PdV work done by the expanding lobes. This is consistent with previous high-resolution simulations of lobe inflation in idealized clusters.

• Shock heating plays a crucial role in converting the jet kinetic energy to thermal energy to inflate jet lobes. Heating by both internal shocks and the bow shock is necessary to explain the lobe energetics, although their relative importance changes with jet power, with internal shocks becoming progressively more important for higher power jets.

• Shock heating also contributes to ICM heating, with higher power jets providing a larger fraction of the injected jet energy to the ICM via shocks. This being said, weak shocks make up a large portion of the heating (especially for lower power jets), and while ICM Mach numbers increase with jet power, even for the highest power jet explored, they rarely exceed |$\mathcal {M}=3$|⁠, consistent with observations of jet-driven shocks in galaxy clusters.

• Cluster weather displaces jet lobes from their original trajectory, aids their disruption, and drives mixing of hot lobe material with the ICM, especially after the jet has switched off. Mixing is more effective for lower power jets, which are more susceptible to the cluster weather.

• While a significant fraction of the jet energy ends up as potential energy of gas in the ICM, we find that systems with higher power jets retain a higher fraction as thermal energy while low power jets retain a larger fraction as kinetic energy. This is consistent with the picture that higher power jet lobes, which are more ‘explosive’, convert a larger fraction of kinetic to thermal energy via shocks during more rapid lobe inflation, while lower power jet lobes are more likely to simply displace material.

• Overall, the impact on cluster properties is subtle. During lobe inflation, the jets can compress cold, dense ISM gas, and enhance instantaneous SFRs, although this results in only small differences for the long-term integrated stellar mass. However, the jet feedback is able to heat the ICM and gently raise the mean core entropy even long after the jet switches off, with energy in excess of that in the non-jet run being comparable to the total injected jet energy. Given that |$\lower.5ex\hbox{$\,\, \buildrel\gt \over \sim \,\,$}50{{\ \rm per\ cent}}$| of the injected energy ultimately ends up as gravitational potential energy of the system, differences in the radial ICM profiles appear small at later times as the energy is redistributed and the cluster core relaxes. This is consistent with observations of cluster radial profiles that do not show sharp variations even during AGN outbursts, pointing to feedback needing to be a gentle process.

ACKNOWLEDGEMENTS

MAB and DS acknowledge support by the European Research Council (ERC) starting grant 638707 ‘BHs and their host galaxies: co-evolution across cosmic time.’ DS further acknowledges support from the Science and Technology Facilities Council (STFC). This research used: The Cambridge Service for Data Driven Discovery (CSD3), part of which is operated by the University of Cambridge Research Computing on behalf of the STFC DiRAC HPC Facility (www.dirac.ac.uk). The DiRAC component of CSD3 was funded by BEIS capital funding via STFC capital grants ST/P002307/1 and ST/R002452/1 and STFC operations grant ST/R00689X/1. The DiRAC@Durham facility is managed by the Institute for Computational Cosmology on behalf of DiRAC. The equipment was funded by BEIS capital funding via STFC capital grants ST/P002293/1 and ST/R002371/1, Durham University, and STFC operations grant ST/R000832/1. DiRAC is part of the National e-Infrastructure.

DATA AVAILABILITY

The data underlying this article will be shared on reasonable request to the corresponding author.

Footnotes

REFERENCES