© R. S. Beckmann et al. 2019

1. Introduction

One of the most striking features of the nearby Perseus cluster, NGC 1275, is the extended filamentary Hα emission nebula in its centre (Lynds 1970; Heckman et al. 1989; Crawford & Fabian 1992; Conselice et al. 2001; Hatch et al. 2007; Fabian et al. 2008). Harbouring up to 5 × 10¹⁰ M_⊙ of cold gas (Salomé et al. 2006), this emission nebula has a filamentary morphology, with individual filaments up to 40 kpc long and only 70 pc wide (Conselice et al. 2001; Fabian et al. 2016). Within the extended, filamentary Hα emission, dense clumps of molecular gas have been observed (Salomé et al. 2006; Lim et al. 2012), and some filaments show signs of star formation (Fabian et al. 2008; Canning et al. 2010, 2014). Larger observational samples show that the Perseus cluster is not the only object to house such Hα emission nebulae, with many massive galaxy clusters showing similar features (Crawford & Fabian 1992; Heckman et al. 1989; McDonald et al. 2010, 2012; Olivares et al. 2019) in their centre. Where does this gas come from, and what causes its characteristic filamentary morphology?

Finding cold gas in cluster centres is not unexpected. As cooling times in the intra-cluster medium (ICM) of massive galaxy clusters are short, a massive cooling flow of the order of 100 − 1000 M_⊙ yr⁻¹ is expected to develop in the cluster centre (Fabian 1994). However, observed star formation rates (SFRs) in clusters are of the order of only 1−10% of the naive cooling rate inferred from X-ray observations (McDonald et al. 2018). Clusters must therefore contain a heating source which prevents overcooling and slows down star formation. Many clusters show evidence for extended jets powered by active galactic nuclei (AGN), which are inflating large cavities in the ICM whose power is sufficient to offset cooling (McNamara & Nulsen 2007; Rafferty et al. 2006; Fabian 2012). Via the self-regulation cycle, which consists of cold gas feeding the AGN, which in turn powers a jet, which then inflates cavities that heat the ICM, AGN are expected to play a decisive role in determining the cooling and star formation properties of the cluster (see McNamara & Nulsen 2007; Fabian 2012, for a review). This picture of self-regulation cycles from AGN jets is getting increasing support from hydrodynamical simulations both in an idealised (Cattaneo & Teyssier 2007; Gaspari et al. 2011; Li & Bryan 2014a) and in a cosmological context (Dubois et al. 2010).

The cospatiality of the Hα emission nebula with the AGN jets and bubbles suggests that the AGN might not only control global cooling properties of the cluster but also be more directly responsible for the morphology of the existing dense gas (Salomé et al. 2006; Russell et al. 2017; Vantyghem et al. 2017, 2018; McKinley et al. 2018; Tremblay et al. 2018). The often complex line-of-sight velocity field of the nebula in Perseus also suggests that this gas is not merely free-falling, or rotationally supported (McDonald et al. 2012; Gendron-Marsolais et al. 2018), but most likely interacts with the turbulence injected by the AGN jets and buoyantly rising bubbles (Fabian et al. 2003; Hatch et al. 2006; Revaz et al. 2008). However, with only line-of-sight velocity information, the three-dimensional velocity pattern of gas is difficult to ascertain.

Simulations by McCourt et al. (2012) and Sharma et al. (2012) show that even for a globally thermally stable ICM (required to avoid overly strong cooling flows) dense gas can condense out of the hot ICM via local thermal instabilities for sufficiently low values of t_cool/t_ff. Here, t_cool is the local cooling time

(1)

where n_i, n_e and n are the ion, electron and total number density respectively, T is the temperature and Λ the cooling rate. The free fall time is

(2)

where g is the local gravitational acceleration and r is the radius from the cluster centre. Condensation into multi-phase can take place when locally t_cool/t_ff <  1, but it is also observed for larger values of the radial t_cool/t_ff profile due to the turbulence and inhomogeneities injected by uplifting hot gas from the cluster centre via AGN driven feedback processes (Voit et al. 2017; Voit 2018). It has been confirmed observationally that molecular gas is observed at the minima of t_cool/t_ff profiles (Hogan et al. 2017; Pulido et al. 2018; Olivares et al. 2019), with some of these authors stressing that only t_cool determines condensation rates as the growth of linear perturbations is largely independent of the geometry of the gravitational potential (Choudhury & Sharma 2016).

Simulations show that the turbulence injected by AGN feedback can cause the local thermal instabilities predicted by McCourt et al. (2012), but struggle to reproduce the observed morphologies, with dense gas having either a very clumpy morphology (Li & Bryan 2014b; Yang & Reynolds 2016a) or settling into a massive central disk (Gaspari et al. 2012; Li & Bryan 2014a,b; Prasad et al. 2015). While the latter has is observed in some clusters, such as in Hydra-A (Hamer et al. 2014; Olivares et al. 2019), only a small central disk is observed in Perseus (Nagai et al. 2019). The dense gas morphology therefore seems to sensitively trace the energy balanced in the ICM.

One feature of these clusters is that the cold gas is expected to rain down on the AGN in a cold and chaotic fashion (Gaspari et al. 2013; Voit & Donahue 2015; Voit et al. 2017), so the cold gas accreted by the supermassive black hole (SMBH) is expected to lack coherent angular momentum, which in turn could lead to a reorientation of the SMBH spin axis over time. In this paper, we investigate the impact of this chaotic dense accretion on the formation of further gas by explicitly tracing the spin of the SMBH, and using this SMBH spin axis as the axis for the AGN driven jet (Dubois et al. 2014). In contrast to existing simulations, which rely on a fixed jet axis with pre-defined precession within a narrow jet cone (Li & Bryan 2014a; Yang & Reynolds 2016a; Ruszkowski et al. 2017; Li et al. 2017; Prasad et al. 2018; Martizzi et al. 2019; Wang et al. 2019), the spin driven approach used here is able to inject turbulence over a larger volume of the cluster centre, and to respond dynamically to the evolving dense gas morphology throughout the simulation.

In this paper, we will investigate the formation and time evolution of dense gas structures in a Perseus-like cluster under the influence of a spin-driven jet, with a particular focus on clump dynamics. The simulations are introduced in Sect. 2. A general overview of results is given in Sect. 3.1, the jet evolution is studied in Sect. 3.2 and the clump properties are investigated in Sect. 3.3. A detailed look at the role of uplifting in clump properties and dynamics is given in Sect. 3.4, and the impact of condensation is studied in Sect. 3.5. A discussion of results can be found in Sect. 4, and conclusions are summarised in Sect. 5.

2. Simulation setup

This paper presents a set of hydrodynamical simulations of isolated galaxy clusters, produced with the adaptive mesh refinement code RAMSES (Teyssier 2002).

2.1. Technical details and refinement

For the simulations presented here, the Euler equations were solved with the second order MUSCL-Hancock scheme, which computes Godunov fluxes using an approximate HLLC Riemann solver and a MinMod total variation diminishing scheme to reconstruct the interpolated variables. The Courant factor was set to a value of 0.8.

The simulation was performed in box of size 8 Mpc with a root grid of 64³, and then adaptively refined to a maximum resolution of 120 pc. Refinement proceeded according to several criteria. We used a quasi-Lagrangian criterion: when a cell contains a mass greater than 3.5 × 10⁹ M_⊙, it was refined (and it was derefined if it contains less than 0.125 this). We also used a Jeans-based criterion: a cell was refined until the local Jeans length was > 4 times the cell size. To refine regions of interest to this work, we also employed two additional refinement criteria. First, the cell containing the SMBH was forced to be refined at the maximum resolution. Second, a passive scalar variable was injected by the SMBH jet with a mass density ρ_scalar equals to that of the gas ρ_gas, which was advected with the gas and marked regions affected by SMBH feedback. The scalar decayed exponentially, with a decay time of t_jet = 10 Myr to ensure that the scalar traces only recent AGN feedback events. After testing different decay times, we confirm that the results do not sensitively depend on this value. To resolve the AGN bubbles, cells were allowed to be further refined when the scalar fraction exceeded ρ_scalar/ρ_gas >  10⁻⁴, equivalent to 92 Myr since the last feedback event, and its relative variation from one cell to another exceeded 10⁻². The latter two refinement criteria ensure that the regions affected by AGN feedback, including the hot, low density bubbles which would de-refine under a purely Lagrangian refinement scheme, remained maximally refined over a reasonable duration of the jet propagation and mixing with hot ICM.

2.2. Initial conditions

The initial conditions for dark matter (DM) and gas consisted of a cored Navarro-Frenk-White profile:

(3)

where r_s = r₂₀₀/c is the scale radius, r_core = 20 kpc is the core radius, ρ_s = ρ_cδ₂₀₀ is the density scaling of the profile, with ρ_c the critical density of the Universe. The rescaling factor , where f(c) = ln(1 + c)−c/(1 + c), rescaled the profile to the radius at which the average density of the profile is 200 times the mean density of the Universe. f_gas is the gas fraction of the halo, here taken to be 15%. The halo had a concentration parameter c = 6.8, and a virial velocity v₂₀₀ = 1250 km s⁻¹.

DM particles had a mass resolution of 3.7 × 10⁸ M_⊙ and were distributed using the DICE code (Perret et al. 2014). The profile was truncated at a radius 2.2 Mpc, for a total DM halo mass of 3.4 × 10¹⁴ M_⊙. DM particle were live and able to move under gravity, allowing the DM potential to respond to the evolution of the cluster core throughout the simulation.

Gas was initiated in hydrostatic equilibrium assuming a gas fraction of 15%, distributed according to the profile of the DM (see Eq. (3)), and then allowed to cool. As part of the initial conditions, turbulence was injected into the gas with a velocity dispersion of 15 km s⁻¹, but no rotation was added to the halo. This small initial velocity dispersion in the hot gas serves to break the symmetry of the initial conditions. Metallicity was initially set to 0.3 Z_⊙ throughout, and the SMBH sinkparticle was placed in the centre of the halo. No stars were added as part of the initial conditions. In order to avoid edge effects, the halo was placed in a sufficiently large box (8 Mpc on a side), and initiated with a gas density of 9.8 × 10⁻⁸ cm⁻³ outside of the truncation radius of the halo.

2.3. Cooling

The metal-dependent cooling of the gas was followed using the tabulated values of Sutherland & Dopita (1993) down to 10⁴ K. The cooling function was extended below 10⁴ K with the fitting functions from Rosen & Bregman (1995). Solar abundance ratios of the elements were assumed throughout, independent of the overall metallicity.

2.4. Star formation and stellar feedback

Star formation proceeded according to a combined density and temperature criterion, with star formation permitted in cells with hydrogen number density of n_H >  1 H cm⁻³ and temperature T <  10⁴ K. The mass resolution of stars was n_Hm_pΔx³/X_H = 5.6 × 10⁴ M_⊙, where X_H = 0.74 is the fractional abundance of hydrogen. The SFR density proceeded according to a Schmidt law , where ρ is the gas density, t_ff is the gas free-fall time, and ϵ_* = 0.1 is the constant efficiency of star formation.

Stellar feedback was included in the form of type II supernovae only. We used the energy-momentum model of Kimm et al. (2015) with each stellar particle releasing an energy of e_*, SN = m_*η_SN10⁵⁰ erg at once after 10 Myr, where η_SN = 0.2 corresponds to the mass fraction of the initial mass function for stars ending up their life as type II supernovae, and m_* is the stellar particle mass. These explosions also enriched the gas with metals with a constant yield of 0.1. Metals were treated as a single species and were advected as a passive scalar.

2.5. SMBH accretion and feedback

AGN feedback from the central SMBH was followed using the model from Dubois et al. (2010) with several modifications that include the self-consistent evolution of the SMBH spin (Dubois et al. 2014) and the spin-dependent feedback efficiency (Dubois et al., in prep.).

A SMBH “sink” particle was placed at the centre of the halo as part of the initial conditions, with a mass of 3.4 × 10⁸ M_⊙. The SMBH was then free to move across the grid throughout the simulation. To compensate for unresolved dynamical friction from the stars within the host galaxy, an analytic drag force was applied to the sink particle according to Pfister et al. (2019). We did not model the equally unresolved gas drag explicitly as the difficulty in measuring the relative velocity between the sink and the turbulent ISM introduces too many numerical artifacts in the SMBH trajectory (see Beckmann et al. 2019). A particular worry was the SMBH getting attached to its own feedback and being ejected from the central galaxy, which we avoided here by not using a sub-grid prescription for the gas drag.

The SMBH accreted according to the Bondi-Hoyle-Lyttleton accretion rate

(4)

where , and are the mass weighted local average density, sound speed and relative velocity between the gas and the SMBH. All quantities were measured within a sphere with radius 4Δx_min centred on the instantaneous position of the SMBH, with the SMBH free to move across the grid. Δx_min is the size of the smallest resolution element of the simulation. Accretion was not limited to the Eddington accretion rate.

The AGN feedback was modelled with jets following the injection method from Dubois et al. (2010). At each feedback event, feedback energy

(5)

(where c is the speed of light) was injected as kinetic energy within all cells contained in a cylinder of radius 0.4 kpc and height 0.8 kpc. The cylinders was aligned with the SMBH spin axis. The efficiency η_MAD is a spin-dependent efficiency obtained from magnetically arrested disc (MAD) simulations from McKinney et al. (2012), which has a minimum at a spin of 0, and a maximum at a spin of 1. The SMBH spin-up rate is taken from the same simulations. The AGN jet was always taken to be aligned with the SMBH spin axis, and the conditions for SMBH-disc alignment in misaligned grid-scale gas angular momentum (with that of the SMBH spin) is obtained by Lense-Thirring considerations (see Dubois et al. 2014, for details). As the spin-axis changes self-consistently throughout the simulation, we did not need to add any explicit precession to the jet, as it naturally arises from the chaotic nature of the cold gas accretion onto the SMBH (Gaspari et al. 2013), which regularly changes the SMBH spin direction over time (see Sect. 3.2).

As mentioned in Sect. 2.1, a passive scalar of density ρ_scalar = ρ_gas was injected within the feedback cylinder at each feedback event, where ρ_gas is the gas density. This scalar then decayed exponentially with a decay time of 10 Myr, allowing cells recently affected by the AGN jet to be identified. Therefore, with the AGN passive scalar quantity, one can define an age for the gas that has been impacted by the AGN, with t_AGN = −10ln(Y_scalar) Myr, where Y_scalar = ρ_scalar/ρ_gas.

2.6. Tracer particles

To follow the dynamical history of gas in the simulation we employ Monte-Carlo based tracer particles from Cadiou et al. (2019). These tracer particles are a significant improvement over classical “velocity”-based tracer particles, in particular in regions with strongly converging flows such as cold gas condensation and gas collapsing under self-gravity. We set up 2 × 10⁸ tracer particles, with each particle tracing a gas mass of 4 × 10⁵ M_⊙. They were initially distributed according to the gas density profile in the initial conditions, out to a radius of 200 kpc.

3. Results

3.1. Cluster evolution

As can be seen in Fig. 1, which shows the gas density, temperature, AGN age (see Sect. 2.5) and the gas radial velocity at various times, the gas in the cluster develops a multi-phase structure with a complex morphology that evolves significantly over the course of the simulation. The hot gas in the intra-cluster medium, which has temperatures in the range 0.09 − 1131 keV (10⁶ − 1.3 × 10¹⁰ K), cools down and condensates into dense clumps and filaments within the central 50 kpc of the cluster, with an average temperature of the dense gas of 4.0 × 10⁻⁴ keV (4.6 × 10³ K). This dense gas falls towards the centre where it feeds the central SMBH and thereby triggers the AGN jet, which, in return, interacts with existing dense gas and stirs turbulence into the hot gas, generating hot outflows with outflow velocities up to 3.5 × 10⁴ km s⁻¹. As the radio jet is oriented along the SMBH spin axis, which in turn is updated according to the chaotic cold accretion onto the central SMBH (Gaspari et al. 2013; Voit et al. 2017), the jet continuously re-orients throughout the simulation (see Sect. 3.2 for details). As a result, the shapes of the jet relics indicated by the “AGN age” also change significantly over time.

Fig. 1. Projections of (from left to right) density, temperature, radial velocity and the time since a cell has been affected by AGN feedback, tAGN, at five different points in time. Radial velocity and temperature are weighted by tAGN. Radial velocity is measured in 3D space with the SMBH at the origin, and negative velocities are inflowing. Contours are based on the plot of tAGN, and are drawn at 10 (solid), 50 (dashed) and 200 (dotted) Myr. The location of the SMBH is marked by a white cross, and black contours in the right hand column denote the outline of dense gas structures. The white dotted line lies along the instantaneous jet axis, which is plotted to be exactly 10 kpc long in 3D space. The shorter it appears, the more it is aligned with the line of sight of the image (here taken to be the z-axis of the box at all times).

More quantitatively, Fig. 2 shows that gas begins to cool after approximately 100 Myr, equivalent to the initial cooling time of gas in the cluster centre as set by the initial conditions. Dense gas, for the remainder of the analysis, is defined to be gas with a maximum temperature of 10⁶ K. By 139 Myr, the dense gas mass first exceeds 10⁹ M_⊙, and the cluster enters a cyclic behaviour where dense gas repeatedly builds up to a total mass in excess of 2 × 10¹⁰ M_⊙ before being reduced to closer to 2 × 10⁹ M_⊙.

Fig. 2. Top panel: time evolution of cluster properties including stellar mass Mstar, SMBH mass MSMBH and gas mass Mgas. Middle panel: AGN luminosity, X-ray luminosity of hot gas within 50 kpc of the cluster centre, and the dense mass again for comparison. Bottom panel: SFR, as well as the dense gas mass again for comparison, for both the fiducial simulation and for a companion simulation without metal cooling for gas with T >  104 K (see text for details). Dense gas is defined to be gas with a temperature at or below Tdense = 106 K, hot, diffuse gas with a temperature above that. White and grey background colours show the heating and cooling dominated regimes of the fiducial simulation.

We have split the evolution of the cluster into two regimes using the total dense gas mass. A cooling dominated regime, when the total dense gas mass of the cluster increases (marked with a grey background in Fig. 2), and a heating dominated regime, when the total dense gas mass of the cluster decreases. The regime of the cluster is evaluated using the smoothed derivative of the mass of dense gas M_gas, dense. The total dense gas mass in the cluster can be reduced in a number of different ways: dense gas can be consumed in star formation, accreted onto the SMBH or destroyed via hot winds or shocks driven by AGN feedback.

AGN activity (see second panel of Fig. 2) is highest during the heating-dominated phase, with maxima in dense gas followed by maxima in AGN activity within 50 Myr or less. These peaks in AGN activity destroy dense gas in the cluster, causing the AGN to enter a low feedback state until the dense gas mass has had time to build up again. Only a small fraction of the gas is directly accreted by the SMBH, as can be seen by the fact that the mass increase of the SMBH mass in the top panel of Fig. 2 is much smaller than the decrease in dense gas mass over the equivalent period of time.

As can be seen in the bottom panel of Fig. 2, the SFR varies strongly over time, following the general trends set by the total dense gas mass in the cluster. There are clear bursts of star formation in the cooling dominated interval. This suggests that a significant amount of the dense gas is directly consumed by star formation. At peaks of up to 1000 M_⊙ yr⁻¹, the SFR of our simulated cluster is extremely high in comparison to observations, which for equivalent mass clusters report SFRs in the range 1 − 100 M_⊙ yr⁻¹ (O’Dea et al. 2008). The dense gas mass, by contrast, falls close to the 10¹⁰ − 10¹¹ M_⊙ observed in Perseus (Bridges & Irwin 1998; Salomé et al. 2006; Mittal et al. 2015). The SFR might be so elevated in comparison to observations because gas is cooling too efficiently to start with, or because gas is being converted too efficiently into stars once cooled. The latter is discussed further in Sect. 4.1.

If gas in the cluster is cooling too efficiently, too much gas is transitioning from the hot, diffuse phase to the dense phase. The X-ray luminosity of the central 50 kpc of the simulated cluster are in the range of 1.2 − 5.3 × 10⁴⁵ erg s⁻¹, with the observed values for Perseus of 1.26 × 10⁴⁵ s⁻¹ (Ebeling et al. 1996) at the lower end of that range. While the initial conditions were chosen to reproduce observed profiles, the emitted X-ray luminosity increases due to the gas cooling in the cluster centre.

One limitation of the simulation presented here is the lack of cosmological context, which means that heating due to turbulence injected by large-scale phenomena, such as galaxy mergers or anisotropic accretion, is absent. In addition, being purely hydrodynamical, the simulation disregards effects such as magnetic fields and other non-thermal energy sources such as cosmic rays, which could heat the gas and provide an extra pressure support against collapse on small scales.

Currently, we rely on equilibrium cooling with an initially uniform metallicity of 0.3 Z_⊙ everywhere, based on observations of the metallicity in the outskirts of Perseus by Werner et al. (2013). By 1 Gyr, the volume weighted metallicity in the central 50 kpc of the hot ICM has risen to 0.36 Z_⊙ due to stellar feedback. While this is higher than the initial value, it still falls below the value of 0.6 observed by Schmidt et al. (2002). One possibility is that equilibrium cooling assumed here over-estimates the contribution of metal cooling at high temperatures. X-rays emitted by the AGN could dissociate metals in high temperature gas, reducing their contribution to cooling (Dubois et al. 2011; Agertz et al. 2013).

If radiative transfer and non-equilibrium processes were included, the hard X-rays emitted by the AGN would be able to photo-ionize some important metal coolants further so that their contribution to cooling is reduced (e.g. Gnedin & Hollon 2012; Segers et al. 2017). As metal line cooling is the dominant cooling channel for gas between 10⁵ and 10⁷ K, shutting down metal cooling would hamper the transition of gas from the hot, diffuse phase to the dense phase. To test this hypothesis, we ran a simulation using a cooling function in which the metal cooling function is modified by a kernel

(6)

which effectively shuts off metal cooling for gas with temperatures above T >  10⁴ K. As can be seen in the bottom panel of Fig. 2, while the initial cooling is delayed in comparison to the fiducial simulation, SFRs remain high even with truncated metal line cooling and the evolution of dense gas is qualitative indistinguishable between the two simulations. We therefore conclude that metal line cooling is not the root cause of the over-cooling reported here. It is more likely that the over-cooling occurs due the absence of non-thermal energies from cosmic rays, which are expected to be able to offset as much as 60% of the thermal cooling in a cluster environment (Pfrommer 2013; Jacob & Pfrommer 2017a,b; Ruszkowski et al. 2017), while only contributing on the percent level to the overall pressure (Reimer et al. 2004; Brown et al. 2011). Due to the large reservoir of heat in cluster outskirts, thermal conductivity in massive clusters can also be an efficient process to bring balance back to the hot cooling gas in the centre of clusters (Narayan & Medvedev 2001; Ruszkowski & Oh 2010; Yang & Reynolds 2016b; Kannan et al. 2017). These avenue of investigation will be explored in future work.

3.2. Jet evolution and turbulence in the cluster

One important difference between the work presented here, and previous works on the subject (Li & Bryan 2014a,b; Yang & Reynolds 2016a; Ruszkowski et al. 2017; Li et al. 2017; Cielo et al. 2018; Martizzi et al. 2019; Wang et al. 2019) is that our jet axis is not fixed throughout the simulation, nor do we add explicit precession. Instead, the spin evolution of the SMBH not only determines the feedback energy but also, crucially, the direction of the jet, as the jet axis is taken to be aligned with the SMBH spin axis, and the SMBH spin axis is continuously updated according to the angular momentum of accreted gas.

Figure 3 shows that the direction of the jet explores the full parameter space of the simulation, repeatedly traversing the full range of both polar and azimuthal angles (0 ≤ θ_jet <  360° and 0 ≤ ϕ_jet <  180°). This is a consequence of the chaotic angular momentum accreted by the SMBH. As can be seen in the bottom two panels of Fig. 3, the angular momentum of the accreted gas varies extremely rapidly, both in θ and in ϕ, as clumps rain down on the SMBH from all directions. As the SMBH spin evolution is a continuous measure, it varies more slowly than the angular momentum of the accreted gas. The only time both the gas angular momentum and the SMBH spin direction settle occurs in the period t = 820 − 950 Myr, when both the SMBH spin and the angular momentum have θ ≈ 90° and ϕ close to zero (the apparent large gap in ϕ at this time is a feature of the coordinate system chosen. ϕ = 2° to ϕ = 178° only represents a rotation of 4°, as both 0° and 180° are aligned with the z-axis of the box). At this time a rotating central gas disc forms around the SMBH, as can be seen in Fig. 4, which drives jet bubbles out to more than 70 kpc from the cluster centre. Dense clumps continue to exist at larger radii, but are preferentially found outside the region recently affected by AGN feedback (see solid grey contours on the image).

Fig. 3. Spin evolution of the SMBH, showing the spin magnitude (top panel), and the two angles defining the SMBH axis (second and third panel) and the angular momentum of the accreted gas at that particular timestep (bottom two panels). The angles are measured in the box frame, and are defined to be the same as in polar coordinates, where θjet is measured in the x − y plane of the box (shown in Fig. 1) and ϕjet is the angle with the z-axis (the line of sight in Fig. 1). Angles are measured in the range 0 ≤ θ <  360° and 0 ≤ ϕ <  180°. Discontinuous jumps from just below the upper end of the range, to just above the lower end of the range, or vice versa, are due to the cyclic nature of the coordinate system. The top three panels show both the fiducial simulation, and a second, identical simulation initiated with a higher spin value. White and grey background colours show the heating and cooling dominated regimes of the fiducial simulation.

Fig. 4. Projection plots at t = 874.5 Myr, showing the central gas disc in the cluster: Top row: density projections of the cluster centre along two different lines of sight, bottom left: composite X-ray image, using the same X-ray bins as in Fig. 5, bottom right: jetscalar weighted temperature projection. Contours mark tAGN = 10 and 50 Myr (solid, dashed). The SMBH location is marked by a black cross, and the jet direction is shown by a white dotted line in the right hand panels only.

Our jets self-consistently produce a three-dimensional distribution of fat feedback bubbles seen in Fig. 5, without the need for adding an ad hoc precession or reorientation of the jet (as done in e.g. Li & Bryan 2014a,b; Yang & Reynolds 2016a; Ruszkowski et al. 2017; Li et al. 2017; Cielo et al. 2018; Martizzi et al. 2019; Wang et al. 2019). Firstly, the jet reorientation due to spin helps to self-regulate the cooling flow in clusters (Cielo et al. 2018) by more uniformly redistributing the energy in the hot gas as long as the reorientation is moderate (i.e. not too close to mimicking isotropic energy input, see Gaspari et al. 2012). Secondly, the reorienting jet has important consequences for the distribution of turbulence in the cluster centre, as over time a much larger volume is directly affected by the AGN jet. However, the bubbles shown here are less round and more broken up than observed X-ray cavities in clusters. This is due to the fact that in the absence of viscosity and magnetic fields, strong Rayleigh-Taylor and Kelvin-Helmholtz instabilities at the bubble surface break up bubbles prematurely and shorten their overall lifetime (Ogiya et al. 2018).

Fig. 5. Synthetic composite X-ray images of the cluster centre, with 0.3−1.2 keV in red, 1.2−2 keV in green and 2−7 keV in blue, to match the image of Perseus in Fabian et al. (2006), towards the beginning, middle and end of the simulation. Each channel is scaled to highlight fainter features. Each image is 100 kpc across.

While we explicitly track the spin evolution of the SMBH, as described in Sect. 2.5, the magnitude of the SMBH spin remains small throughout, as can be seen in Fig. 3, with a maximum spin parameter of 0.08. This is partially a consequence of the model chosen, as the MAD jet model always preferentially reduces the spin of the SMBH. This low spin value in turn has consequences for the jet direction, as the jet axis is aligned with the SMBH spin axis. Due to the low spin value of the SMBH, the chaotic angular momentum of accreted gas (see bottom two panels of Fig. 3), driven by the chaotic infall of the clumps, is able to significantly realign the spin axis throughout the simulation.

As can be seen in Eq. (5), the feedback energy of the SMBH is determined by the feedback efficiency η_MAD, which in turn is determined by the SMBH spin. Due to the consistently low spin-values, the simulation presented here has an average luminosity-weighted feedback efficiency of only 0.046.

To test the consequences of a higher initial spin value of the SMBH, we ran a companion simulation to our fiducial simulation. The only difference between the two was that the companion simulation had an initial SMBH spin value of 0.8. As can be seen in the top panel of Fig. 3, the SMBH spin persistently decreases over the course of the simulation, until it converges with the fiducial simulation after ∼500 Myr. While the spin is high, the jet changes direction very slowly in comparison to the fiducial simulation, as the high angular momentum of the rapidly spinning SMBH makes reorientation more difficult. Once the spin has dropped below 0.4, the jet direction changes more rapidly and the two simulations become statistically indistinguishable. The bubbles remain comparatively fat even in the absence of precession. This is due to the fact that our jets are very light and hot, and therefore over-pressurized in comparison to the background medium. While injected bimodally, the bubbles quickly expand outwards into the surrounding medium. We note that the absence of magnetic fields, whose wound-up helical structure along the jet is expected to keep it confined over kpc scales (see Pudritz et al. 2012, for a review), will have contributed to the fatness of the bubbles. We therefore postpone a comparison between bubble structures in a high spin and a low spin case to future, magnetised simulations.

3.3. Dense gas structures

3.3.1. Quantifying clump morphology

As can be seen visually in Fig. 1, the dense gas in the cluster centre can be found in clumps of a wide range of sizes and shapes. A clump is defined here to be a connected volume of space, for which all cells have a minimum density of 1 H cm⁻³ and a maximum temperature of 10⁶ K. All properties are measured by summing over all cells contained within a given clump. Tracer particles are associated with a particular clump if they are contained within the clump volume at the point of measurement.

To quantify this parameter space, we measured the physical extent of individual clumps using the following methodology:

1. Find the centre of mass for each clump by summing over all cells contained within the clump, treating each cell as a point mass located at the cell centre.

2. Calculate the clump’s mass-weighted reduced inertia tensor using

   (7)

   by summing over all cells n contained within a clump, where x_n, i is the ith coordinate of the nth cell within the clump, measured in the centre of mass frame of the clump. R_n is the nths cells distance from said centre of mass, and m_n is its gas mass.

3. Calculate the physical extent of the major axis r_maj by finding the largest distance between any two cell centres contained within the clump. To this value, Δx_min is added to extrapolate from the cell centres to the cell edges contained within the clumps.

4. Find the axis vectors and axis length ratios using the eigenvalues and eigenvectors of the inertia tensor from Eq. (7).

5. Calculate the median and minor axis length, r_med and r_min respectively, using the axis length ratios from the previous step, and the length of the major axis, r_maj.

6. Calculate the volume of the ellipse defined by the three axes:

   (8)

7. Calculate the volume filling fraction f_V, which is defined to be the ratio of the volume defined by the axis vectors, V_ellipse in Eq. (8), and the sum of the cell volumes contained within the clump:

   (9)

   where V_n is the volume of the nth cell contained in the clump. For solid, round clumps well described by an ellipse, f_V will have a value close to unity. For clumps with a complex morphology, such as bent filaments and three-dimensional networks of filaments and clumps, the volume fraction will be low as the axis vectors used to describe the ellipse mark the total physical extent of the clump along a given axis vector in 3D space, and said ellipse will therefore contain many cells outside the clump.

For further analysis, we split the population of clumps into three categories depending on the length of their major axis relative to the mean major axis of the whole sample, kpc, and the samples standard deviation σ_maj = 1.42 kpc:

1. small clumps have a major axis kpc.

2. big clumps have a major axis length in the range kpc.

3. filaments have kpc.

Some example decompositions according to these criteria can be seen in Fig. 6.

Fig. 6. Example projections of the decomposition of structures into small clumps, big clumps and filamentary structures at three different points in time. Structures are considered to be distinct when not connected in 3D space. From left to right, the plots contain 5, 5 and 9 distinct filaments respectively, partially superimposed due to projection effects.

3.3.2. Clump properties

A variety of bulk clump properties versus axis length are shown in Fig. 7, for the stacked sample of clumps of the whole simulation. As can be seen in column (a), the distribution of major axis lengths ranges from the resolution limit of the simulation to very large, extended objects that have major axes of the order 10 kpc or more. The stacked sample shown here, which contains all objects from all snapshots at all points in time of the simulation, contains 37897 small clumps (87.4%), 4283 big clumps (9.9%) and 1153 filaments (2.7%).

Fig. 7. Clump properties for the whole sample (bottom row) and split into the three structure categories (top two rows). From left to right: clump gas mass Mgas, volume ratio fV, distance between the clump centre of mass and the cluster centre rcentre, bulk velocity vr and gas velocity dispersion within the clump σgas, radial. The probability distributions ϕ in the top row is mass weighted, while the one in the row below is unweighted.

As expected, smaller clumps contain less gas mass (Fig. 7, column (b)), with a minimum gas mass for the current resolution of 5 × 10⁵ M_⊙, and an average value of 1.8 × 10⁷ M_⊙ for small clumps and 1.1 × 10⁸ M_⊙ for large clumps. The population of filaments is much more massive, with an average gas mass of 2.5 × 10⁹ M_⊙. Structures with a mass above 10⁹ M_⊙ are all classified as filamentary. This lower mass limit for gas clumps is determined by the resolution. As we tested with a companion simulation, in which we reduced Δx_min to 30 pc, i.e. a factor 4 smaller than in the fiducial simulation. With this improved resolution, the gas structures fragment further into even smaller clumps, with a new minimum mass of 2.2 × 10³ M_⊙, and a new minimum axis length still approaching the resolution limit. This suggests that the shattering into smaller structures is by no means complete, and with even more resolution, the clumps would continue to break apart, as in the “cloudlet” model by McCourt et al. (2018). However, larger, filamentary structures continued to exist even in the higher resolution simulation.

In terms of shape, smaller clumps have higher values of f_V, so they are indeed much more compact (column (c), Fig. 7). Values of f_V >  1 can occur for compact objects when the axis length for the median and minor axis are under-estimated in comparison to the true extent of the clump which happens mainly for clumps with less than 20 cells. However, the volume of the ellipse used to fit the clump never exceeds that of the sum of the cells contained in the clump by more than 40%. More extended objects have f_V far below unity, which is an indicator of complex morphology. The most clumpy filament produced in this simulation still has f_V <  0.7 so large clumpy structures do not form at any point of the simulation.

Small and big clumps have a similar radial distribution (column (d), Fig. 7) and are preferentially found between 3 and 10 kpc from the cluster centre. Filaments, on the other hand, include both a subsample found at large radii, and a sample of particularly extended structures in the cluster centre, an example of which can be seen in the right hand panel of Fig. 6. This suggests that gas structures merge into larger objects as they reach the cluster centre, consistent with a model in which small structures rain down onto a central massive gas structure. This structure can take the form of a massive gas disk, as for example seen in Li & Bryan (2014a) and briefly also in the simulation presented here (see Fig. 4), or in the form of an extended but not rotationally-supported object such as the one in the right hand panel of Fig. 6, or the gas structures seen in the first, third and fourth snapshot of Fig. 1.

In velocity space, all three populations are similarly distributed (column (e), Fig. 7), with no discernible difference in the unweighted probability distribution of small and medium clumps, as well as filaments. The mass-weighted distribution in the top row shows that all three categories of structures are preferentially infalling (i.e. have v_r <  0). The time-stacked sample of the simulation has an unweighted mean radial velocity of 75 km s⁻¹, with a full width half max of 198 km s⁻¹, where radial velocity is measured in 3D space with the SMBH at the origin. Negative values denote gas falling towards the SMBH. These values are comparable to observed bulk velocities of 100 km s⁻¹ but are at the upper end of observed velocity widths of 100−218 km s⁻¹ for molecular gas in Perseus (Salomé et al. 2008; Hitomi Collaboration 2016; Gendron-Marsolais et al. 2018). By comparison, they fall easily within the range of observed velocity widths for warm ionised gas in massive clusters (Hamer et al. 2016). We note that, in contrast to the observational values, the full width-half max calculated here is calculated across the entire time-stacked sample, not just along the line of sight. While the mean and dispersion values show good agreement with observations, the sample of clumps presented here has an overall larger velocity range than found in cold-gas maps of nearby clusters, which report velocity values across the map in the range of 350 km s⁻¹ at most (Olivares et al. 2019; Gendron-Marsolais et al. 2018).

The velocity dispersion σ_{gas, radial} is defined to be the velocity dispersion of the radial velocities of all resolution elements within an individual clump. It therefore quantifies the range of velocities found within an individual object. Clumpy structures, both small and big, have a low velocity dispersion (column (f), Fig. 7), i.e. a small range of radial velocities, with an average value of just 90 km s⁻¹. The bulk of the filaments, despite major axis lengths of 10 kpc or more, have radial velocity dispersion of less than 200 km s⁻¹ but there is a small population of high-velocity dispersion objects with σ_gas >  200 km s⁻¹, which is preferentially populated by filaments: They make up 28% of the high dispersion objects versus only 2.7% of the total sample.

Dynamically, the clumps are therefore a surprisingly uniform population, despite more than 2 orders of magnitude in size difference, and more than 4 orders of magnitude in mass range. Gas properties across all three populations are also similar, with a temperature range of 10 − 10⁶ K (the latter being the cut-off temperature for the definition of a dense gas structure in this paper), and densities in the range of 1 − 10⁵ H cm⁻³. The bulk of the gas has a temperature around 10⁴ K and a density of 10 − 10³ H cm⁻³. This is not to say that all objects have the same properties at a given point in time, but that all types of objects can be found at all points in phase space at some point throughout the simulation.

The morphology and distribution of objects can vary strongly on a 5 Myr timescale, as can be seen in Fig. 8. Overall, the number of structures at all points in time is dominated by small clumps, which are always the most abundant and make up 87.4% of the time-integrated sample. During some parts of the cooling-dominated phase, they also contain the bulk of the dense gas mass, such as around 400 Myr and at 500−550 Myr. The rest of the time, the bulk of dense gas mass can be found in filaments, despite the fact that they only make up 2.7% of the overall sample by number. Big clumps contain dense gas mass on the order of that contained in the small clumps, but represent 9.9% of the total number of objects.

Fig. 8. Time evolution of (from top to bottom) the number, total dense gas mass, mean gas mass and mean distance to the cluster centre for the three structure categories. Bottom two panels: solid lines show the mean and shaded regions the range from the 10th to the 90th percentile of the distribution. The dashed grey line in all plots shows the AGN luminosity for comparison.

From Fig. 8, strong bursts of AGN feedback are followed by a strong increase in the number of small clumps, as well as an equally strong drop in both the total mass of gas contained in filaments (second panel) and the average mass of gas per filament (third panel). At the same time, the average radial distance between the cluster centre and a clumps centre of mass increases (bottom panel). While the bulk of clumps can usually be found within the central 20 kpc of the cluster, strong AGN outbursts produce clumps at much larger radii, up to 50 kpc from the location of the cluster centre. This suggests a scenario where large objects are being shattered into smaller clumps during their interaction with strong AGN jets, and highlights the importance of the AGN jet not just for slowing down cooling onto the cluster centre but also for the morphology and kinematics of the existing dense gas structures. The details of this interaction will be explored further in the next section.

3.4. Uplifting

Uplifting has been used to explain the unstructured velocity profiles observed in nearby clusters (Pulido et al. 2018; Gendron-Marsolais et al. 2018). When talking about uplifting dense gas in clusters, two different mechanisms need to be distinguished. On the one hand, there is the entrainment of existing dense gas by the AGN driven outflows, which turns previously infalling dense gas into outflowing dense gas, which will be discussed in this section. Alternatively, outflowing dense gas could form via condensation at large radii, when gas is uplifted from the cluster centre by AGN jets, before being deposited at larger radii, where local entropy conditions then allow gas to condense Voit et al. (2017), Voit (2018).

The impact of one interaction between the AGN jet and the dense gas in the cluster centre, namely the outburst at 320−400 Myr, is shown visually in Fig. 9: at t = 323 Myr (top left), the dense gas is predominantly infalling and contained in radially oriented filaments. At this point in time, the filaments contain M_{gas, filaments} = 1.7 × 10¹⁰ M_⊙, i.e. 76% of the total dense gas mass, with an average gas mass per filament of . As the AGN outburst commences, fed by this infalling dense gas (t = 338 − 356 Myr, middle and top right panel), the filaments are broken into small and medium size clumps, and their velocity turns from infalling to outflowing. By 371 Myr (bottom left), gas is predominatly outflowing, and the total mass budget of 8.8 × 10¹⁰ M_⊙ is evenly split between small clumps, medium clumps and filaments. The filaments that continue to exist are much less massive, with an average mass of just .

Fig. 9. a: visual time evolution of one episode of AGN feedback that starts around t = 350 Myr. Only the dense gas is plotted. The colourmap shows the radial velocity of the gas, with negative values denoting infall, with the background colour set to grey for clarity. The location of the SMBH is marked by a cross, and the contours show the extent of the AGN feedback bubbles produced by the feedback event that starts at t = 323 Myr. b: time evolution of dense gas mass contained in the three categories over the same period of time. Vertical grey lines mark the outputs shown in the top panel.

By t = 388 Myr, the gas has reached its largest radial extent for this episode and is beginning to fall back onto the cluster centre in the form of a shower of small, distinct clumps. From 371.9 Myr to 388.5 Myr, the total gas only increases by 5%, from 8.8 × 10¹⁰ M_⊙ to 9.3 × 10¹⁰ M_⊙, but the total number of objects triples as objects continue to break apart, from 244 at 371.9 Myr to 651 individual objects by 388.5 Myr. By this point, small clumps dominate the population, as they represent 94% of objects and contain 64% of the total gas mass, with a further 27% contained in big clumps.

The timeseries of the number of different objects in the top panel of Fig. 8 shows that this behaviour is generic for the cluster presented here. Following a strong feedback outburst, the number of small objects spikes, while the total gas mass and the average mass per filament decrease strongly. At the same time, the average distance for objects of all categories increases as they are ejected from the cluster, with the outermost small clumps being found as far as 40 kpc or more from the cluster centre.

Looking directly at the number of inflowing and outflowing objects, as shown in Fig. 10, strong AGN feedback bursts are followed by a spike in the number of outflowing objects, as larger, filamentary structures are entrained and broken up by the hot winds of AGN feedback and lifted to larger radii. As gas is evacuated from the cluster centre the AGN turns off. The entrained clumps then decelerate under gravity and fall back onto the cluster centre. During this process, they shatter into even smaller components so the number of individual objects continues to increase even after the AGN has become quiescent again. As the small clumps fall back onto the cluster centre, they coalesce and trigger another strong outburst of AGN feedback, which repeats the cycle. The results presented in this paper are similar to work by Yang & Reynolds (2016a), who presented evidence for existing dense gas to be redistributed by the AGN jet. Contrary to their work, the dense gas in the simulations presented here is not indestructible. In our simulations, only 25% of the dense gas survives its interaction with the hot jet. It gets entrained by the AGN driven outflows and lifted to large radii. We note that, with a temperature cut of 10⁶ K, the gas discussed here is equivalent to the ionised dense gas seen in observation, not to the molecular gas. We expect that if we were able to adequately distinguish between ionised warm gas and molecular cold gas, the molecular gas would be much more difficult to uplift by the AGN jet.

Fig. 10. Time evolution of the total number of inflowing and outflowing clumps. The AGN luminosity is shown as a dotted line for comparison. The solid background highlights the event shown in Fig. 9.

This is surprising in the context of work by Klein et al. (1994), who showed that for adiabatic cold structures in hot winds, the drag timescale t_drag ≈ χr_clump/v_wind is always longer than the clump crushing timescale t_cc ≈ χ^1/2r_clump/v_wind, where χ is the density contrast between wind and cold clump, r_clump is the clump radius and v_wind is the relative velocity. It should therefore be impossible to accelerate cold clumps with a hot wind. However, recent work by Gronke & Oh (2018) shows that radiative cooling can replenish the cold clump mass from the hot gas during uplifting and thereby dramatically increase the clump lifetime. Under these assumptions, clumps with radii larger than r_clump >  v_windt_{cool, mixing}/χ, where t_{cool, mixing} is the cooling time in the mixing layer surrounding the cold clumps, should survive the uplifting process, as cooling from the hot to the cold phase replenishes gas faster than cold gas from the clumps is being evaporated. For the simulation presented here, the maximum outflow velocities in the vicinity of clumps is of the order 10⁴ km s⁻¹, the cooling time in the mixing layer around clumps is of the order 0.1 Myr and the density contrast χ ≈ 10⁴. Therefore, clumps with a minimum value of r_clump ≈ 1 pc should survive their interaction with the hot wind, which is much smaller than the smallest cell size of 120 pc. While poorly resolved clumps most likely lack this mixing layer, and are therefore destroyed during the jet interaction, well-resolved cold clouds would be expected to survive their interaction with the hot outflows and become entrained without being destroyed, as shown in Fig. 9. These results are also in agreement with work by Armillotta et al. (2017), who show that the bulk of cold gas in clouds with radii above 250 pc survives being accelerated by a hot wind for 200 Myr. It is however likely that the 25% of dense gas that survives the interaction in our simulations is an overestimate, as work by Sparre et al. (2019) showed that more highly resolved clouds shatter more efficiently during their interaction with hot winds and therefore have shorter overall lifetimes than less resolved clouds.

In comparison to the observed velocity maps for Hα emitting gas in Perseus by Gendron-Marsolais et al. (2018), the velocity maps from our simulation (as shown in Fig. 9) are much more coherent, with clumps either predominantly infalling or outflowing in a given map. In this context we note that the maps in Fig. 9 show an unusual period for our cluster, i.e. the only AGN outburst during which the number of infalling clumps fall almost to zero (see Fig. 10). This episode was chosen for analysis as it illustrates uplifting by AGN feedback particularly cleanly. At other points in time, dense gas can be observed to be inflowing and outflowing at the same time in our simulation, due to the directionality of the jet and the limited width of the jet cone.

It is also important to remember that the observed velocities are line-of-sight velocities, while Fig. 9 shows radial velocities. As can be seen visually in Fig. 11, which shows both radial velocities (top row) and line-of-sight velocities (bottom row) for an inflowing dominated (left column), an outflow dominated (middle panel) and a mixed (right column) point in time, the line-of-sight velocities appear less ordered than the radial velocities. The outflow or inflow dominated nature of the flow (left or middle panels respectively) cannot easily be recovered from line-of-sight velocity maps. This difficulty in distinguishing between flow patterns in the frame of the cluster, and line-of-sight flow patterns, is even more obvious in Fig. 12, which shows the radial velocity probability distributions for the three snapshots in Fig. 11, as well as that for the three line-of-sight velocities (here aligned with the x-axis, y-axis and z-axis of the box respectively). In all three cases, the line-of-sight velocities fail to recover the radial velocity pattern and predict a more gaussian-like pattern with a mean velocity close to zero. The Gaussian distribution of line-of-sight velocities is expected for infalling or outflowing gas distributed roughly spherically around the cluster centre. The chaotic velocity patterns observed in nearby clusters are therefore not necessarily evidence for the absence of coherent radial flows of the gas.

Fig. 11. Density weighted velocity projections of the dense gas at three different points in time. Top row: radial velocity for each snapshot, bottom row: corresponding line of sight velocity (here chosen to be the z-axis of the simulation box).

Fig. 12. Distribution of resolution elements in radial velocity, and line of sight velocity along the x-axis, y-axis and z-axis of the simulation box respectively, for the three snapshots in time shown in Fig. 11.

3.5. Condensation

As first proposed in McCourt et al. (2012), and then shown in idealised simulations by Sharma et al. (2012), dense gas can form out of the hot ICM via local thermal instability, even if the cluster is globally thermally stable. Condensation can happen when locally, t_cool/t_ff falls below 1, and is suppressed for higher values. With sufficient uplifting of gas from the cluster centre, condensation can occur for larger values of the radial t_cool/t_ff profile, up to the range of 10 − 30 (Voit et al. 2017; Voit 2018), as also seen in observations (Hogan et al. 2017; Pulido et al. 2018; Olivares et al. 2019).

In the simulation presented here, we used the tracer particles to estimate the condensation rate of dense gas. As each tracer particle has a unique identification number and traces 2 × 10⁴ M_⊙ of gas mass, the trajectories of tracer particles can be used to track gas flows throughout the simulation. The total mass of gas transferred from the hot, diffuse to dense phase between two simulation outputs can be estimated by counting the number of tracer particles that pass from the diffuse phase to the dense phase between two simulation outputs. The condensation rate Ṁ_condensed is then found by dividing the newly condensed gas mass M_condensed by the time it took to assemble it.

As can be seen in the left hand panels of Fig. 13, our simulation confirms that condensation primarily occurs when t_cool/t_ff <  20. This is somewhat higher than prediction from idealised cooling simulations (Sharma et al. 2012; McCourt et al. 2012), most likely because the hot gas along the jet drives up the spherically averaged cooling time, but in line with observed values (Hogan et al. 2017; Olivares et al. 2019). Profiles of t_cool/t_ff during the cooling dominated phases, which produce the bulk of the condensation, are generally ordered, with a clear minimum around 10 kpc. During heating dominated phases, by contrast, profiles show a much wider range of shapes as gas heated by the AGN rises to large radii in the form of hot bubbles, which significantly increase the cooling time both in the centre and at larger radii. Some condensation continues during the heating dominated phases, and while the condensation remains confined to < 20 kpc from the cluster centres, the values of t_cool/t_ff can be as high as 50 even for actively cooling clusters. We postulate that this continued condensation is due to the multiphase structure of the ICM and the directionality of AGN feedback. Both t_cool and t_ff are calculated for the hot ICM only, and it takes even strong AGN feedback bursts some time to reach large volume filling factors and shut off condensation completely.

Fig. 13. Cluster profiles of the cooling time (tcool) to free fall time (tff) ratio at different snapshots of the simulation. The cluster profiles are sampled each 25 Myr across the full time evolution of the simulation. tcool is calculated for each cell in the simulation, using its instantaneous density, temperature and cooling function as computed by RAMSES. tff is calculated using all mass (DM, gas, stars and the SMBH) contained within a given radius. Profiles are colour-coded by condensation rate (left) or dense gas mass (right), based on the condensation rate and dense gas mass onto clumps at that radius. Snapshots during cooling dominated (top panel) and heating dominated intervals (bottom panel) are plotted separately for clarity.

This hypothesis is confirmed by the condensation time-series in Fig. 14, which shows that condensation is highest towards the minimum of heating-dominated phases and falls to zero as the AGN feedback continues of impact the ICM. Figure 14 also shows that at the end of cooling-dominated phases, condensation occurs preferentially onto filamentary structures, but by the end of heating-dominated phases and the beginning of the next cooling-dominated phases, condensation occurs preferentially onto small and big clumps, in line with the uplifting – shattering – recondensation picture presented in Sect. 3.4.

Fig. 14. Time evolution for the gas condensation rate onto the dense structures in the simulation. See text for how the condensation rate is calculated.

As can be seen in Fig. 14, the total condensation rate of the cluster varies with time, ranging from a minimum of 3 M_⊙ yr⁻¹ at the beginning of cooling dominated intervals to a maximum of up to 1.8 × 10³ M_⊙ yr⁻¹ towards the end of cooling dominated phases. While the bulk of condensation takes place onto filaments, smaller and big clumps dominate when condensation rates are low. As discussed in the context of the clusters SFR in Sect. 3.1, this condensation rate is high in comparison to the observed condensation rate for Perseus, which is in the range of 50 − 100 M_⊙ (Fabian 2012). In future work, we will explore if this over-cooling occurs because of the omission of non-thermal energies from cosmic rays in the work presented here, which are expected to be able to offset as much as 60% of the thermal cooling in a cluster environment (Pfrommer 2013; Jacob & Pfrommer 2017a,b; Ruszkowski et al. 2017).

While the areas of high condensation rate are confined to the minima of the t_cool/t_ff profiles, dense gas can be found over a much wider range of radii (see righthand panels of Fig. 13), and significant amounts of dense gas can also be observed during heating-dominated times. This is due to the fact that existing dense gas free-falls onto the cluster centre from its formation location around 10 kpc, and is uplifted to larger radii due to its interactions with AGN feedback. The location at which dense gas is observed is therefore not a perfect proxy for where it is formed, as the kinematics in active clusters are complex and subject to hysteresis.

This can be seen in more detail when comparing the radial and velocity distributions for stacked samples of newly condensed gas (left panel) and dense gas (right panel) in Fig. 15. While some amount of condensation occurs over the full parameter space of radii and velocities occupied by dense clumps, the distribution in both radius and velocity is different for newly condensed gas and dense gas in general. As shown in both the mass distribution in Fig. 15, and in the probability distributions in Fig. 16, dense gas is preferentially found at the cluster centre, whereas condensation preferentially occurs at larger radii, with a peak of the distribution at 10 kpc. In velocity space, both existing dense gas and new condensation are preferentially infalling, but condensation has a broader distribution towards negative values, with a mean velocity at −155.6 km s⁻¹ for condensation compared to −104.3 km s⁻¹ for dense gas. Overall, only 75.9% of gas is infalling, while 82.1% of condensation occurs onto infalling clumps. This means that while the bulk of newly condensed gas is infalling, with an average condensation rate onto inflowing gas of 2.56 M_⊙ yr⁻¹, there is also evidence for gas condensation onto outflowing clumps, which have an average condensation rate of 1.05 M_⊙ yr⁻¹.

Fig. 15. Phase plot of the total condensation rate (left) and total dense gas mass (right) over a range of radial positions and radial velocities of the clumps. Data shown here is stacked over all clumps at all snapshots of the simulation.

Fig. 16. Probability distribution of clump radius (left) and clump radial velocity (right) weighted by condensation rate and total dense mass respectively.

We therefore conclude that condensation occurs preferentially onto infalling clumps within the radial range of 5 − 15 kpc, but approximately a fifth of all condensation occurs onto outflowing gas.

4. Discussion

In this paper, we have studied the formation, evolution and destruction of dense gas in the centre of a Perseus-like cluster, under the influence of a spin-driven AGN jet. We have particularly focused on the role played by uplifting and condensation in the kinematics and morphology of the dense gas.

4.1. Cooling and star formation

As reported in Sect. 3.1 and shown again in Fig. 17, the SFR of the cluster (solid line) is much higher than observed values for Perseus, such as for example the 71 M_⊙ yr⁻¹ measured by Mittal et al. (2015; dotted line). As discussed in Sect. 3.1, this could be due to an overly high cooling rate of the gas, or because too much of the resulting dense gas is turned into stars. To compare the star formation efficiency of the cluster with observation, Fig. 17 shows both the total dense gas condensation rate Ṁ_cond from Fig. 14 and a naive cooling rate, defined as

(10)

Fig. 17. Cooling flow rates and SFRs (top panel), as well as the resulting star formation efficiencies (bottom panel) using two different measures for the cooling rate, Ṁcond and Ṁcool,naive. Solid lines show time evolution, while dashed lines show time averages. Dotted and dashed-dotted lines show observational values for Perseus from Mittal et al. (2015).

following McDonald et al. (2018), where M_gas(r <  r_cool) is the total gas mass contained within the cooling radius r_cool, which in turn is defined to be the radius at which the cooling time profile, t_cool(r) = 3 Gyr.

As can be seen in Fig. 17, both the time series and the average value for Ṁ_cool,naive are a factor 2−4 higher than the observed value, except during the disc-dominated stage between 850 and 950 Myr. By contrast, the dense gas condensation rate, Ṁ_cond shows significant variablity but has a time-averaged value that is close to the naive observed cooling rate. It is also noticeably lower than the naive cooling rate, suggesting that reheating by the AGN keeps the majority of cooling gas from cooling efficiently and prevents it from condensating into dense gas.

Looking at the resulting star formation efficiencies (ϵ_cool, bottom panel), the average value of ϵ_cool = SFR / Ṁ_cool,naive = 0.19±0.27 is only slightly higher than the ϵ_cool = 0.16 reported by McDonald et al. (2018), but the scatter on this value is large. The error on ϵ_cool given here is equal to one standard deviation of the distribution. Looking at the efficiency of converting dense gas into stars, ϵ_cond = SFR / Ṁ_cond = 0.76 ± 1.37 means that the majority of dense gas is turned into stars. This shows that despite individual cold clumps loosing as much as 75% of their mass during interactions with strong feedback episodes such as the one shown in Fig. 9, only about a quarter of the total dense gas is returned to the hot phase in this manner. Destroying dense gas once it has condensed is therefore not an efficient mechanism to regulate star formation in the cluster. Given the large variation in ϵ_cond, the instantaneous SFR is not a reliable tracer of Ṁ_cond, the cold gas formation rate of the cluster. It is possible that we overestimate the SFR in dense gas, as we use a comparatively simple density-based star formation recipe of the form , which does not take the effects of small-scale turbulence into account, and could therefore be too efficient for the context shown here (Shi et al. 2011; Salomé et al. 2016).

Figure 17 also shows that in general Ṁ_cool,naive > > Ṁ_cond, so the vast majority of gas that cools out of the hot phase does not reach the dense phase, and is instead reheated by the AGN before condensing fully. As already known from the classic cooling flow problem and confirmed here more quantitatively, Ṁ_cool,naive is therefore not a good tracer of the overall cooling budget of the cluster as Ṁ_cond/Ṁ_cool,naive = 0.38 ± 0.27 on average, with the error again denoting a standard deviation.

4.2. Filament lifetimes

One notable result of our simulations is that extended gas structures form preferentially during comparatively AGN quiet times, and are readily destroyed in the interaction with AGN feedback. While this interaction between dense gas and AGN feedback is one of the requirements for effective self-regulation of cooling in the cluster, it also means that the lifetime of filaments is limited by the length of AGN duty cycles. While as much as 25% of the dense gas mass survives the interaction with the hot, AGN driven outflows, larger structures are broken into smaller structures in the process. The result is a volume-filling distribution of small clumps, which are at first outflowing and then fall back onto the cluster centre. Such a clumpy morphology of the dense gas is not supported by observations, which show more extended, filamentary structures (Conselice et al. 2001; Fabian et al. 2006). Two possible explanations come to mind.

One possibility is that the dense filaments are too readily destroyed in our simulations. If physical processes not modelled here, such as notably magnetic fields, could support the filaments against fragmentation, they might survive their interaction with the AGN jet and retain their extended morphology for longer. This hypothesis is supported by work on the survival rate of isolated clumps accelerated by hot, magnetised winds (Shin et al. 2008; McCourt et al. 2015; Xu & Lazarian 2018), which show that magnetised winds draw spherical clouds out into extended, filamentary structures instead of evaporating them or breaking them into smaller clumps. From this point of view, we over-estimate the fragmentation rate of dense filaments into the hot ICM.

The other possibility is that we underestimate the ability of AGN feedback to destroy dense clumps, for example by under-resolving the mixing layers at the outer clump surface (Gronke & Oh 2018), or simply due to lack of resolution to follow the fragmentation process to smaller scales. This hypothesis is supported by our high-resolution companion simulation, which showed that the fraction of dense gas that survives this particular uplifting event falls from 25% at a resolution of Δx_min = 120 pc to 19% at a resolution of Δx_min = 30 pc. The fact that the minimum clump size remains at the resolution limit shows that this process is by no means converged, and higher resolution would likely lead to even smaller clumps and even lower dense gas survival rates. This question has been investigated further by McCourt et al. (2018), who report that for individual clouds accelerated by a hot wind, even a sub-pc scale resolution is insufficient for fragmentation to converge. Based on work by Armillotta et al. (2017), the survival rates for small gas clumps in hot winds is very low, which suggests that we would expect the gas currently contained in our small, compact gas clumps to break into an even large number of even smaller clumps until it evaporates entirely and mixes back into the ICM. From this point of view, we are under-estimating the fragmentation rate of small clumps, as well as under-estimating the ability of the AGN to evaporate dense gas.

4.3. The width of filaments

Throughout this paper, we have shown that extended dense gas structures readily form in the cluster centre. While our filamentary dense gas structures show maximal extents of 1−10 kpc, in agreement with observations (Conselice et al. 2001), many of our structures appear much wider than the observed 70 pc. Resolution will play a role in determining the width of the filaments, particularly for very thin filaments which currently have a width close to the resolution limit, such as the long, thing structures seen in the left two panels of Fig. 6. A comparison simulation with higher resolution of Δx_min = 30 pc, run for only a span of 15 Myr, produced thinner filaments than the fiducial simulation at 120 pc. However, many filaments seen in the fiducial simulation, such as for example the extended structures in the right hand panel of Fig. 6, are well resolved at the current resolution and therefore not influenced by improvements in resolution.

One process not modelled here, which is thought to play an important role in the morphology of filaments, is anisotropic thermal conduction along magnetic field lines. In the presence of anisotropic thermal conduction, in combination with magnetic fields, the characteristic thermal collapse length scale (the field length) becomes much larger along field lines than perpendicular to it (Field 1965), as thermal energy is preferentially redistributed along field lines. Collapse therefore preferentially occurs perpendicular to magnetic field lines, smearing spherical collapse out along magnetic field lines. Isolated simulations have shown that in the presence of magnetic fields, local thermal instabilities do indeed produce more extended filamentary gas structures (McCourt et al. 2012; Ji et al. 2018; Xu & Lazarian 2018) compared to more clumpy dense gas for the purely hydrodynamical runs. While this process could help smear dense, round clumps into long, extended filaments, it is unlikely to make the existing filaments thinner. Understanding why the filaments reported here take their particular shapes, and how their morphology might change in the presence of magnetic fields and cosmic rays, will be the subject of future work.

Another limitation of our work is that with many structures shown here at the resolution limit of the simulation, it will be impossible to resolve the detailed internal structure observed for filaments, which consist of dense molecular clumps surrounded by an Hα envelope (Salomé et al. 2006, 2011). With a more complex internal structure and gas dynamics, we would expect the energy balances of filaments to change, with as of yet poorly understood consequences for their morphology.

5. Conclusions

In this paper, we have investigated the formation and evolution of dense gas in the centre of a Perseus-like cluster under the influence of a spin driven AGN jet, using hydrodynamical simulations.

We showed that:

1. Under the influence of the AGN jet, the cluster undergoes repeated cycles of cooling dominated phases, when dense gas builds up in the cluster centre, and heating dominated phases, when the total amount of dense gas decreases. Cycle lengths are on the order of 100 Myr, but show significant variation (Sect. 3.1).

2. For low SMBH spin values, the chaotic cold accretion onto the cluster centre is able to continuously reorient the spin axis, with characteristic reorientation timescales of the order of 10 Myr, allowing the jet to sweep out the full parameter space in both polar and azimuthal angle (Sect. 3.2).

3. The morphology of dense gas is highly variable throughout the simulation, with between 20 and 620 individual dense structures present at a given point in time (Sect. 3.3).

4. Major axis lengths of individual clumps range from the resolution limit of the simulation up to more than 30 kpc. Larger clumps have more complex, filamentary morphologies than smaller objects, which tend to be rounder and compacter (Sect. 3.3).

5. We find evidence for uplifting of existing dense gas by the AGN, with dense gas defined to have a maximum temperature of 10⁶ K. During a strong feedback episode, larger, infalling structures fragment into smaller clumps under the influence of the hot outflows driven by the AGN. In the process, they lose up to 75% of their gas mass and become entrained and ejected from the cluster centre (Sect. 3.4).

6. Despite these high mass loss rates for individual clumps, 75.7% of the total dense gas is turned into stars during the course of the simulations. Despite individual clumps loosing up to 75% of their mass during interactions with AGN jets, destruction of dense gas via AGN feedback is therefore not an efficient channel to regulate star formation in clusters (Sect. 4.1).

7. A second round of fragmentation into even smaller clumps occurs at the top of the ballistic orbit, before surviving dense clumps fall back onto the cluster centre where they re-coalesce into larger objects (Sect. 3.4).

8. Condensation takes place preferentially when t_cool/t_ff <  20, which occurs primarily during cooling dominated phases of the cluster, and in the radial range of 5 − 15 kpc. Heating dominated phases see more disturbed profiles of t_cool/t_ff without a clear minimum as the ICM is heating by the AGN feedback (Sect. 3.5).

9. Dense gas continues to be observable even during heating-dominated phases, and is preferentially found at smaller radii than condensation, i.e. at r <  5 kpc, but can be found as far out as 30 kpc due to uplifting. The presence of dense gas is therefore not a reliable tracer for condensation (Sect. 3.5).

10. While 82.1% of condensation of gas from the hot ICM onto dense clumps occurs on infall, there is also evidence for continued condensation for outflowing gas, with outflowing dense clumps having an average condensation rate of 1.06 M_⊙ yr⁻¹, compared to 2.56 M_⊙ yr⁻¹ for infalling clumps (Sect. 3.5).

11. Both direct uplifting of dense gas and condensation of gas from the hot, diffuse to the dense phase in outflowing gas has been invoked to explain the unstructured velocity maps observed in nearby clusters. While we find evidence for both mechanisms, and confirm a general lack of rotation in the dense gas, we also caution that the observed line of sight velocities fail to show coherent radial flow patterns even when they are present in the dense gas (Sect. 3.4).

12. Neither the naive cooling rate Ṁ_cool,naive nor the SFR are reliable observational tracers of the cold gas formation rate Ṁ_cond (Sect. 4.1).

Acknowledgments

The authors thank Maxime Trebitsch and Marta Volonteri for useful discussion, and the referee, Mark Voit, for his excellent comments and useful input. This work was supported by the ANR grant LYRICS (ANR-16-CE31-001 1) and was granted access to the HPC resources of CINES under the allocation A0040406955 made by GENCI. This work has made use of the Horizon Cluster hosted by Institut d’Astrophysique de Paris. We thank Stéphane Rouberol for smoothly running this cluster for us. Visualisations in this paper were produced using the YT PROJECT (Turk et al. 2011).

References

All Figures

Fig. 1. Projections of (from left to right) density, temperature, radial velocity and the time since a cell has been affected by AGN feedback, tAGN, at five different points in time. Radial velocity and temperature are weighted by tAGN. Radial velocity is measured in 3D space with the SMBH at the origin, and negative velocities are inflowing. Contours are based on the plot of tAGN, and are drawn at 10 (solid), 50 (dashed) and 200 (dotted) Myr. The location of the SMBH is marked by a white cross, and black contours in the right hand column denote the outline of dense gas structures. The white dotted line lies along the instantaneous jet axis, which is plotted to be exactly 10 kpc long in 3D space. The shorter it appears, the more it is aligned with the line of sight of the image (here taken to be the z-axis of the box at all times).

In the text

Fig. 2. Top panel: time evolution of cluster properties including stellar mass Mstar, SMBH mass MSMBH and gas mass Mgas. Middle panel: AGN luminosity, X-ray luminosity of hot gas within 50 kpc of the cluster centre, and the dense mass again for comparison. Bottom panel: SFR, as well as the dense gas mass again for comparison, for both the fiducial simulation and for a companion simulation without metal cooling for gas with T >  104 K (see text for details). Dense gas is defined to be gas with a temperature at or below Tdense = 106 K, hot, diffuse gas with a temperature above that. White and grey background colours show the heating and cooling dominated regimes of the fiducial simulation.

In the text

Fig. 3. Spin evolution of the SMBH, showing the spin magnitude (top panel), and the two angles defining the SMBH axis (second and third panel) and the angular momentum of the accreted gas at that particular timestep (bottom two panels). The angles are measured in the box frame, and are defined to be the same as in polar coordinates, where θjet is measured in the x − y plane of the box (shown in Fig. 1) and ϕjet is the angle with the z-axis (the line of sight in Fig. 1). Angles are measured in the range 0 ≤ θ <  360° and 0 ≤ ϕ <  180°. Discontinuous jumps from just below the upper end of the range, to just above the lower end of the range, or vice versa, are due to the cyclic nature of the coordinate system. The top three panels show both the fiducial simulation, and a second, identical simulation initiated with a higher spin value. White and grey background colours show the heating and cooling dominated regimes of the fiducial simulation.

In the text

Fig. 4. Projection plots at t = 874.5 Myr, showing the central gas disc in the cluster: Top row: density projections of the cluster centre along two different lines of sight, bottom left: composite X-ray image, using the same X-ray bins as in Fig. 5, bottom right: jetscalar weighted temperature projection. Contours mark tAGN = 10 and 50 Myr (solid, dashed). The SMBH location is marked by a black cross, and the jet direction is shown by a white dotted line in the right hand panels only.

In the text

	Fig. 5. Synthetic composite X-ray images of the cluster centre, with 0.3−1.2 keV in red, 1.2−2 keV in green and 2−7 keV in blue, to match the image of Perseus in Fabian et al. (2006), towards the beginning, middle and end of the simulation. Each channel is scaled to highlight fainter features. Each image is 100 kpc across.
In the text

	Fig. 6. Example projections of the decomposition of structures into small clumps, big clumps and filamentary structures at three different points in time. Structures are considered to be distinct when not connected in 3D space. From left to right, the plots contain 5, 5 and 9 distinct filaments respectively, partially superimposed due to projection effects.
In the text

Fig. 7. Clump properties for the whole sample (bottom row) and split into the three structure categories (top two rows). From left to right: clump gas mass Mgas, volume ratio fV, distance between the clump centre of mass and the cluster centre rcentre, bulk velocity vr and gas velocity dispersion within the clump σgas, radial. The probability distributions ϕ in the top row is mass weighted, while the one in the row below is unweighted.

In the text

	Fig. 8. Time evolution of (from top to bottom) the number, total dense gas mass, mean gas mass and mean distance to the cluster centre for the three structure categories. Bottom two panels: solid lines show the mean and shaded regions the range from the 10th to the 90th percentile of the distribution. The dashed grey line in all plots shows the AGN luminosity for comparison.
In the text

Fig. 9. a: visual time evolution of one episode of AGN feedback that starts around t = 350 Myr. Only the dense gas is plotted. The colourmap shows the radial velocity of the gas, with negative values denoting infall, with the background colour set to grey for clarity. The location of the SMBH is marked by a cross, and the contours show the extent of the AGN feedback bubbles produced by the feedback event that starts at t = 323 Myr. b: time evolution of dense gas mass contained in the three categories over the same period of time. Vertical grey lines mark the outputs shown in the top panel.

In the text

	Fig. 10. Time evolution of the total number of inflowing and outflowing clumps. The AGN luminosity is shown as a dotted line for comparison. The solid background highlights the event shown in Fig. 9.
In the text

	Fig. 11. Density weighted velocity projections of the dense gas at three different points in time. Top row: radial velocity for each snapshot, bottom row: corresponding line of sight velocity (here chosen to be the z-axis of the simulation box).
In the text

	Fig. 12. Distribution of resolution elements in radial velocity, and line of sight velocity along the x-axis, y-axis and z-axis of the simulation box respectively, for the three snapshots in time shown in Fig. 11.
In the text

Fig. 13. Cluster profiles of the cooling time (tcool) to free fall time (tff) ratio at different snapshots of the simulation. The cluster profiles are sampled each 25 Myr across the full time evolution of the simulation. tcool is calculated for each cell in the simulation, using its instantaneous density, temperature and cooling function as computed by RAMSES. tff is calculated using all mass (DM, gas, stars and the SMBH) contained within a given radius. Profiles are colour-coded by condensation rate (left) or dense gas mass (right), based on the condensation rate and dense gas mass onto clumps at that radius. Snapshots during cooling dominated (top panel) and heating dominated intervals (bottom panel) are plotted separately for clarity.

In the text

	Fig. 14. Time evolution for the gas condensation rate onto the dense structures in the simulation. See text for how the condensation rate is calculated.
In the text

	Fig. 15. Phase plot of the total condensation rate (left) and total dense gas mass (right) over a range of radial positions and radial velocities of the clumps. Data shown here is stacked over all clumps at all snapshots of the simulation.
In the text

	Fig. 16. Probability distribution of clump radius (left) and clump radial velocity (right) weighted by condensation rate and total dense mass respectively.
In the text

	Fig. 17. Cooling flow rates and SFRs (top panel), as well as the resulting star formation efficiencies (bottom panel) using two different measures for the cooling rate, Ṁcond and Ṁcool,naive. Solid lines show time evolution, while dashed lines show time averages. Dotted and dashed-dotted lines show observational values for Perseus from Mittal et al. (2015).
In the text