© The Authors 2026

1. Introduction

Over cosmic time, small primordial fluctuations have grown into large-scale structures, hosting most baryons in the form of diffuse, hot plasma – the circumgalactic, intragroup, and intracluster medium (CGM, IGrM, ICM; White & Frenk 1991; Kravtsov & Borgani 2012; Tumlinson et al. 2017). These gaseous halos act as stratified atmospheres shaped by turbulence, gravity, cooling, and feedback, giving rise to a ‘cosmic weather’ (Gaspari et al. 2020). Multi-wavelength observations increasingly reveal that this hot gas undergoes multiphase condensation, forming warm, neutral, and cold components (e.g. Fabian et al. 2008; Tremblay et al. 2016; Maccagni et al. 2018, 2021; Morganti et al. 2023), yet the physical mechanisms behind this process and its role in galaxy evolution remain open questions.

Most galaxies host a supermassive black hole (SMBH, Kormendy & Ho 2013), which plays a key role in their baryon cycle. Through powerful outflows and jets, SMBHs inject energy into their host halos, regulating gas cooling and heating, star formation, and the galaxy luminosity function (e.g. King & Pounds 2015; Fiore et al. 2017; Iyer et al. 2025). The ensemble of phenomena tied to a SMBH is typically known as an active galactic nucleus (AGN; Lynden-Bell 1969; McNamara & Nulsen 2012). Central open problems are how gas is transported from galactic scales down to the SMBH, how angular momentum is removed, and how AGN feedback couples back to the multiphase medium.

Bridging physical processes across more than nine orders of magnitude, from the megaparsec scale cosmic environment down to the black hole gravitational radius, remains a major theoretical and numerical challenge (Naab & Ostriker 2017; Vogelsberger et al. 2020). Accretion is governed by the interplay of turbulence, cooling, heating, angular momentum transport, and magnetic fields (Balbus & Hawley 1998; Yuan & Narayan 2014). Connecting theory and observation therefore requires one to follow gas across this entire range of scales, with multiphase structures in hot, stratified halos providing key constraints on how cooling, turbulence, and feedback regulate inflow (e.g. Fabian 2012).

In the absence of turbulence and feedback, accretion would proceed as a steady, nearly spherical inflow, concentrating cold gas and star formation at the very centre (Fabian 1994). In reality, ultraviolet observations of brightest cluster galaxies reveal a far more complex picture: star formation occurs in diverse morphologies – often clumpy, in knots, or filamentary – consistent with star formation triggered or enhanced by AGN jet activity (e.g. Donahue et al. 2015; Fogarty et al. 2015). A similar connection between star-forming structures and radio jets is also supported by recent multi-wavelength observations (e.g. Reefe et al. 2025). These observations suggest that gas can cool and form stars at large radii (≈10 − 30 kpc) from the centre of massive halos, along extended, filamentary, and multiphase structures. Millimetre and sub-millimetre observations of molecular gas (e.g. Russell et al. 2019; Olivares et al. 2019; Castignani et al. 2025; Toni et al. 2026) reveal cold gas reservoirs and, in several systems, cold filaments whose morphology closely follows that of the radio jets. Such structures are interpreted as filaments of accreting cold gas, hosting in situ star formation and likely regulated by AGN feedback. Consistent jet – cold gas alignments are also observed in high-redshift systems, where ALMA observations reveal extended cold molecular reservoirs in the CGM aligned with radio jets, supporting a close coupling between AGN activity and cold gas condensation (e.g. Li et al. 2021; Walter et al. 2025).

In realistic intracluster, intragroup, and circumgalactic environments the hot plasma is, however, turbulent, intermittently heated, and radiatively cooling via bremsstrahlung and metal-line emission. In this stratified medium, density perturbations can enter the nonlinear thermal instability regime: overdense fluctuations cool faster than they are mixed and restored, so they grow, decouple from the hot phase, and condense into warm and cold structures. Classical Bondi accretion (Bondi 1952) offers a useful idealised baseline for steady, spherical, and adiabatic inflow, but fails to capture the complexity of realistic accretion in turbulent, multiphase environments. A promising new theory is chaotic cold accretion (CCA; Gaspari et al. 2013, 2015, 2017; Gaspari & Sądowski 2017), in which turbulence and thermal instability promote the condensation of cold clouds that rain stochastically onto the SMBH. This inflow naturally produces variability and stochastic fueling, allowing the black hole to respond dynamically to the state of the surrounding atmosphere. CCA offers a coherent explanation for the observed variability, obscuration, and multi-wavelength signatures of AGN, supported by observations with ATCA, ALMA, MUSE, Chandra, XMM-Newton, JWST, and HST (Maccagni et al. 2014; Voit & Donahue 2015; McDonald et al. 2018; Tremblay et al. 2018; Waters & Proga 2019; Juráňová et al. 2020; Maccagni et al. 2021; McKinley et al. 2022; Temi et al. 2022; Ubertosi et al. 2023, 2025; Olivares et al. 2022, 2023, 2025; O’Sullivan et al. 2024; Eskenasy et al. 2024). Most simulations, however, resolve inflows only down to ∼1–100 pc – several orders of magnitude above the black hole event horizon – and have primarily focused on massive cluster-scale halos. A fully consistent, multi-scale model of SMBH fueling from halo to sub-parsec scales is therefore still lacking and remains a major challenge for the field.

Direct evidence of turbulence has been observed with Hitomi in the core of the Perseus cluster, where spectral line broadening revealed a line-of-sight velocity dispersion of σ_v = 164 ± 10 km s⁻¹ in gas with a temperature kT ≃ 4 keV (Hitomi Collaboration 2016). More recently, XRISM Resolve observations have enabled a similar detection in the Abell 2029 cluster, measuring σ_v = 169 ± 10 km s⁻¹ in a hotter atmosphere with kT ≃ 8 keV (XRISM Collaboration 2025). These results build on earlier XMM–Newton/RGS constraints on line broadening in cool-core systems (e.g. Sanders & Fabian 2013; Pinto et al. 2015). Assuming isotropic turbulence, these measurements correspond to three-dimensional Mach numbers of ℳ ≈ 0.2 − 0.3¹, placing both systems in a weakly turbulent regime and indicating that the ICM hosts mainly subsonic turbulent motions. These direct measurements provide robust constraints on the level of subsonic turbulence in cool-core clusters, complementing indirect estimates derived from surface brightness fluctuations in high-resolution X-ray images, such as those from Chandra, which can be used to infer turbulence in the ICM (Gaspari & Churazov 2013; Hofmann et al. 2016; Khatri & Gaspari 2016; Dupourqué et al. 2024; Romero et al. 2025). By analysing the power spectrum of surface-brightness and thermodynamic fluctuations, one can estimate the amplitude and injection scale of predominantly subsonic turbulence and relate these motions to the underlying density and pressure perturbations. Such measurements commonly indicate ℳ ∼ 0.1 − 0.5 turbulence, spanning clusters to lower-mass groups. As future missions like NewAthena (Nandra et al. 2013) and AXIS (Russell et al. 2024; Koss et al. 2025) are expected to extend these measurements to larger cluster samples and broader radial ranges, a more complete picture of ICM turbulence will emerge.

Numerical hydrodynamical simulations have become an indispensable tool in modern astrophysics, allowing us to model and interpret a wide range of nonlinear, multi-scale processes that are inaccessible to direct experimentation or analytical treatment. They have been successfully applied to study galaxy formation and large-scale structure (e.g. Vogelsberger et al. 2014; Dubois et al. 2014; Pillepich et al. 2019; Hopkins et al. 2023), the evolution of the interstellar and circumgalactic medium (e.g. Springel et al. 2005; Marinacci et al. 2017; Hopkins et al. 2018; Barbani et al. 2023, 2025; Zhang et al. 2025), as well as BH accretion and AGN outflows (e.g. Gaspari et al. 2012; Chen et al. 2024; Fournier et al. 2024, 2025; Grete et al. 2025; Sotira et al. 2025; Jennings et al. 2025). In the AGN context, high-resolution simulations are particularly effective at capturing the nonlinear interplay between cooling, turbulence, and feedback that governs the formation of multiphase gas and the stochastic fueling of SMBHs across a wide range of scales.

This work is part of the BLACKHOLEWEATHER project (PI: Gaspari), which aims to build a unified multi-physics framework for SMBH feeding and feedback across cosmic environments, from galaxies to groups and clusters. BLACKHOLEWEATHER combines high-resolution numerical experiments, synthetic multi-wavelength observables, and theory-guided diagnostics to connect multiphase halo weather, accretion variability, and AGN self-regulation within a single consistent picture. Building on the controlled stratified-halo setups that established the CCA framework (Gaspari et al. 2013, 2017), we present the first of two complementary papers investigating the turbulence-driven condensation and accretion cascade down to sub-pc scales in a hot, stratified atmosphere representative of group-scale halos. We focus on the transitional meso-scale (parsecs to kiloparsecs) that links halo rain to the inner inflow. In this paper, we quantify the emergent multiphase morphology, phase structure and thermodynamics as a function of the turbulence regime. The companion paper (Barbani et al. 2026, B26b hereafter) addresses feeding variability and kinematics in the same set of simulations, characterising the temporal statistics and power spectra of accretion and the kinematical imprints via key CCA diagnostics (k-plots and 𝒞-ratios).

To this end, we carry out a suite of idealised accretion simulations onto SMBHs, designed to isolate the competition between radiative cooling and turbulent mixing. We use ATHENAPK, an adaptive mesh refinement (AMR) hydrodynamical code designed for GPU architectures, to follow the gas from halo scales down to sub-parsec radii. By varying the strength of turbulence driving, we assess how this key process controls the multiphase condensation cascade from the hot to the molecular phase and the resulting mode of SMBH feeding, resolving cold, clumpy inflows near the SMBH and their role in regulating its growth. Within BLACKHOLEWEATHER, activities are organised into thematic work packages (WPs) spanning multiscale feeding and feedback (see main diagram in Gaspari et al. 2020). This paper is part of WP2 (meso- to micro-scale feeding), focused on turbulence-driven condensation and accretion. Complementary efforts currently underway include simulations of micro-scale AGN jet feedback (Cammelli et al. 2026a,b, C26a,b), simulations exploring the impact of BH spin on macro-scale feedback (Piana et al. 2026a,b, P26a,b), and simulations investigating the role of dust physics in the multiphase feeding rain (Barbani et al., in prep.).

This paper is structured as follows. In Section 2, we describe the simulation setup, including the physical models, initial conditions (ICs), and numerical methods. In Section 3 we present the main results, focusing on the formation and evolution of the multiphase structures and their role in feeding the SMBH. In Section 4 we interpret these outcomes in terms of BH weather states and propose an evolutionary link between regimes. In Section 5, we compare our findings with previous observational and numerical studies. Finally, Section 6 summarises our conclusions and outlines the directions for future work.

2. Numerical methods

The simulations performed in this work are made using the AMR code ATHENAPK, an open source, performance-portable magneto-hydrodynamical code, directly descending from ATHENA++ (Stone et al. 2020), based on the block-structured adaptive mesh refinement framework PARTHENON (Grete et al. 2023) and the performance portability programming model KOKKOS (Edwards et al. 2014; Trott et al. 2021). In this way, ATHENAPK is compatible with a wide range of architectures and with both CPUs and GPUs, demonstrating high efficiency and scalability up to 73 728 GPUs.

The equations of hydrodynamics that are integrated by the code in conservative and Eulerian form are as follows:

$$\begin{array}{cc}\hfill \frac{\partial \rho }{\partial t}+\mathrm{\nabla }·(\rho \mathit{v})& =0,\hfill \end{array}$$(1)

$$\begin{array}{cc}\hfill \frac{\partial (\rho \mathit{v})}{\partial t}+\mathrm{\nabla }·(\rho \mathit{v}\otimes \mathit{v}+P\mathbb{I})& =-\rho \mathrm{\nabla }\mathrm{\Phi }+\rho {\mathit{f}}_{\mathrm{turb}},\hfill \end{array}$$(2)

$$\begin{array}{cc}\hfill \frac{\partial E}{\partial t}+\mathrm{\nabla }·(E+P)\mathbf{v}& =-\rho \mathbf{v}·\mathrm{\nabla }\mathrm{\Phi }-\mathcal{C}+{\mathcal{S}}_{\mathrm{turb}},\hfill \end{array}$$(3)

where ρ is the density, v is the velocity, P is the pressure, 𝕀 is the identity tensor, Φ is the total gravitational potential, f_turb is the turbulent driving acceleration (see Section 2.5) and E is the total energy density of the gas, defined as E = e + ρv²/2, where e is the internal energy density. 𝒞 = n²Λ(T) is the cooling rate, which depends on the cooling function Λ(T) (see Section 2.4) and the number density of the gas n, and 𝒮_turb = ρv ⋅ f_turb.

The code employs a finite-volume Godunov (1959) scheme with second-order accuracy in both space and time. Given the highly multiphase chaotic medium, for stability, we use the piecewise linear method (van Leer 1979) for spatial reconstruction and a second-order Runge–Kutta scheme (Butcher 2008) for time integration. Fluxes at cell interfaces are calculated by solving the Riemann problem using the approximate Harten-Lax-van Leer contact (HLLC) solver (Toro et al. 1994), which accurately captures contact and shear discontinuities. Given the complexity of astrophysical phenomena, we applied a first-order flux correction to ensure numerical stability and accuracy near discontinuities (Brüggen et al. 2023). In gas cells where non-physical values, such as negative densities or temperatures, arise, fluxes are recalculated using a forward Euler time integration, first-order (piecewise constant) reconstruction, and the Local Lax–Friedrichs Riemann solver. This fallback approach ensures physically meaningful solutions without the need to impose artificial density or temperature floors.

2.1. Grid structure

The simulation is performed using a static mesh refinement (SMR) setup. The computational domain is organised as a hierarchy of nested Cartesian grid blocks with increasing resolution towards the centre of the domain, where the highest spatial accuracy is required. Each grid block contains 32³ cells, and refinement proceeds by subdividing parent blocks into child blocks with twice the resolution in each spatial direction. The refinement hierarchy is static, defined at initialisation, and remains fixed throughout the simulation. It is centred on the origin and extends radially outwards, with higher refinement levels covering progressively smaller volumes, while coarser grids are employed in the outer halo to reduce computational cost and still capture the large-scale gas dynamics (see Figure 1). The root grid spans the entire simulation volume and contains 128³ cells, with a total of ℓ = 12 refinement levels. The cell size at refinement level ℓ is given by

$$\begin{array}{c}\hfill \mathrm{\Delta }{x}_{\ell }=\frac{L}{N_0·2^{\ell }},\end{array}$$(4)

Fig. 1. Simulation grid and SMBH sink scheme. The domain is decomposed into blocks (shown as a green square), which are recursively refined towards the centre of the domain (each block contains 323 cells). The decreasing block size increases the effective resolution in the central region, where the black hole sink is located. The black hole sink (shown in black) removes mass within a radius of four finest-level cells at each time step (an example of cell is shown as a red square).

where L = 50 kpc is the box size and N₀ = 128 is the number of cells per side of the root grid. This configuration yields a maximum spatial resolution of Δx ≃ 0.1 pc within the central 3.2 pc region. Consistency across refinement levels is ensured through appropriate prolongation and restriction operations at all refinement boundaries. This SMR configuration enables accurate modelling of multiscale processes such as gas cooling, inflow, and feedback-driven turbulence in the central region of the halo.

2.2. SMBH sink

Following Gaspari et al. (2013) setup, to represent the presence of a central SMBH, we introduce a sink region at the centre of the computational domain. The sink is implemented as a fixed spherical region of radius r_sink = 4Δx_min = 0.4 pc, centred at the origin, and acts as an effective inner boundary condition for gas accretion.

At each timestep, gas within r_sink is removed from the simulation by resetting its thermodynamic and kinematic properties to ρ_sink = 10⁻³⁰ g cm⁻³, T_sink = 1 K, and v_sink = 0 km s⁻¹. This prescription avoids nonphysical gas pile-up at the resolution limit of the simulation and maintains numerical stability. Thanks to the sub-parsec spatial resolution achieved in the central region, the simulations resolve gas dynamics well within the Bondi radius of a group-scale SMBH (i.e. r_B ∼ 1–30 pc, see Section 2.3). However, the physical processes occurring at event-horizon scales remain unresolved. The sink therefore provides a physically motivated approximation for the unresolved accretion flow, absorbing gas that would otherwise be expected to feed the SMBH. The gravitational softening length ϵ (Plummer softening, i.e. $\mathrm{\Phi }=-GM/\sqrt{r^2+{ϵ}^2}$) is chosen to be smaller than r_sink, ensuring that gravitational forces are well resolved within the sink region. The SMBH mass is fixed and does not evolve during the simulation. The total mass that would be accreted over the simulated time span is negligible compared to the assumed black hole mass, and does not affect the gas dynamics or the gravitational potential in the central region.

2.3. Initial conditions

The ICs consist of a gaseous halo in hydrostatic equilibrium inserted in a static gravitational potential.

Since galaxy groups are the fundamental building blocks of the cosmic web and host the bulk of present-day baryons in halos, we focus on a representative, intermediate-mass system rather than an extreme, massive cluster. This also helps the computational costs where the cooling time is significantly shorter (i.e. t_cool ≈ 10 − 20 Myr in the central kiloparsec of our simulations).

Gravity is modelled using a static gravitational potential that includes contributions from a dark matter halo, a central dominant group galaxy, and a SMBH (see also Fournier et al. 2024). The gravitational acceleration is defined as g = −∇Φ. For the dark matter halo, we assume a Navarro–Frenk–White (NFW) density profile (Navarro et al. 1997), for which the radial component of the gravitational acceleration is

$$\begin{array}{c}\hfill g_{\mathrm{NFW}}(r)=\frac{G}{r^2}\frac{{\mathrm{M}}_{\mathrm{NFW}}[ln(1+\frac{r}{R_{\mathrm{NFW}}})-\frac{r}{r+R_{\mathrm{NFW}}}]}{ln(1+c_{\mathrm{NFW}})-\frac{c_{\mathrm{NFW}}}{1+c_{\mathrm{NFW}}}},\end{array}$$(5)

where M_NFW is the virial mass of the halo, R_NFW is the scale radius, and c_NFW is the halo concentration parameter.

The scale radius R_NFW is related to the total halo mass and the characteristic density ρ_s, and is defined by the following:

$$\begin{array}{c}\hfill R_{\mathrm{NFW}}={\left(\frac{{\mathrm{M}}_{\mathrm{NFW}}}{4\pi {\rho }_s[ln(1+c_{\mathrm{NFW}})-c_{\mathrm{NFW}}/(1+c_{\mathrm{NFW}})]}\right)}^{1/3},\end{array}$$(6)

which ensures that the mass enclosed within R₂₀₀ = c_NFWR_NFW matches M_NFW. The characteristic density ρ_s determines the amplitude of the NFW density profile and is calculated as

$$\begin{array}{c}\hfill {\rho }_s=\frac{200}{3}{\rho }_{\mathrm{crit}}\frac{c_{\mathrm{NFW}}^3}{ln(1+c_{\mathrm{NFW}})-c_{\mathrm{NFW}}/(1+c_{\mathrm{NFW}})},\end{array}$$(7)

where the critical density of the Universe, ρ_crit, sets the normalisation, and is defined as

$$\begin{array}{c}\hfill {\rho }_{\mathrm{crit}}=\frac{3{\mathrm{H}}_0^2}{8\pi G},\end{array}$$

with H₀ being the Hubble parameter. This formulation allows the NFW profile to be fully specified by the two parameters M_NFW and c_NFW, which are typically derived from cosmological simulations or observational constraints.

The central dominant elliptical galaxy (cD) is modelled with a Hernquist profile (Hernquist 1990), which gives a gravitational acceleration g(r) of

$$\begin{array}{c}\hfill g_{\mathrm{cD}}(r)=G\frac{{\mathrm{M}}_{\mathrm{cD}}}{R_{\mathrm{cD}}^2}\frac{1}{{(1+\frac{r}{R_{\mathrm{cD}}})}^2},\end{array}$$(8)

where R_cD = 10 kpc and M_cD = 1.4 × 10¹¹ M_⊙, similar to NGC 5044, where recurrent evidence for CCA has been found (e.g. Temi et al. 2022; Rajpurohit et al. 2025).

The gravity of the SMBH is derived from its mass as a point-like source as

$$\begin{array}{c}\hfill g_{\mathrm{SMBH}}(r)=\frac{G{M}_{•}}{r^2},\end{array}$$(9)

where the SMBH mass is set to M_• = 2.8 × 10⁸ M_⊙, representative of a central black hole in a low-redshift galaxy group and consistent with the seeding prescriptions adopted in semi-analytic and cosmological models of black hole formation (e.g. Piana et al. 2021; Cammelli et al. 2025; Piana & Pu 2025). Therefore, the total gravitational acceleration is given by g_tot = g_NFW + g_cD + g_SMBH.

The entropy profile is defined following a power law as in the ACCEPT cluster database (Cavagnolo et al. 2009), with values rescaled to represent a galaxy group,

$$\begin{array}{c}\hfill K(r)=K_0+K_{100}{(r/100\mathrm{kpc})}^{{\alpha }_K},\end{array}$$(10)

where K₀ is the central entropy, K₁₀₀ is the entropy at r = 100 kpc, and α_K is the slope of the power law.

Knowing the gravitational acceleration g(r) and the entropy K(r) profiles, we can derive the density profile ρ(r), assuming hydrostatic equilibrium,

$$\begin{array}{c}\hfill \frac{1}{\rho }\frac{\mathrm{d}P}{\mathrm{d}r}=-g(r),\end{array}$$(11)

which can be rewritten in terms of entropy as

$$\begin{array}{c}\hfill \frac{\mathrm{d}}{\mathrm{d}r}[K(r)\rho {(r)}^{\gamma }]=-\rho (r)g(r).\end{array}$$(12)

To derive the density, we first need to close the system defining the electron density n_e(r_fix) at a fixed radius r_fix = 10 kpc. The derived density and temperature profiles are shown in Figure 2. The pink circles correspond to observations of electron density in galaxy groups taken from the ACCEPT catalogue, while the solid black lines represent the profiles derived from our simulation. The simulated density profile lies within the observed scatter and broadly reproduces the characteristic slope and normalisation of the intragroup and intracluster medium over more than two orders of magnitude in radius. In the right panel, the temperature profile of our ICs also agrees well with observational trends, following the universal shape identified by Kim et al. (2020) from Chandra observations of early-type galaxies, when both the temperature and radius are normalised by the peak temperature (T_peak) and the virial radius (R_vir), respectively. All parameters relevant for the definition of our ICs are listed in Table 1.

Fig. 2. Thermodynamic properties of the simulated galaxy group at t = 0 Myr (black lines). Left: radial gas number density profile compared with observed electron densities from the ACCEPT database (Cavagnolo et al. 2009, pink circles). Right: normalised temperature profile (T/Tmax) as a function of r/Rvir, compared with early-type galaxies from the Chandra Galaxy Atlas (Kim et al. 2020, pink circles). The ICs are consistent with observed systems, and the temperature profile follows a nearly universal shape with a broad peak around r ≃ 0.04 Rvir. The vertical dashed line marks the maximum radius contained in the computational domain (rmax = 43.3 kpc, corresponding to half of the box diagonal).

2.4. Radiative cooling

Radiative processes are the main cause of gas cooling in the Universe and are a fundamental process for the formation and evolution of galaxies. The cooling function Λ is shown in Figure 3. For T ≥ 10^4.2 K, we adopt the high-temperature cooling curve of Schure et al. (2009) with solar metallicity. Their model assumes an optically thin plasma in collisional ionisation equilibrium and includes radiative losses from a very extensive set of metal-line transitions in the ultraviolet and soft X-ray bands, together with free–free and recombination continuum emission. The cooling function accounts for thousands of transitions from the main astrophysical elements (H, He, C, N, O, Ne, Mg, Si, S, Ar, Ca, Fe, Ni), providing a physically consistent treatment of radiative losses in warm and hot gas.

Fig. 3. Radiative cooling function at solar metallicity used in our simulations. The curve shows Λ(T) as a function of temperature, with tabulated values from Schure et al. (2009) for T ≥ 104.2 K and the analytic fit given by Eq. (13) for T < 104.2 K. Cooling function includes metal-line, recombination, bremsstrahlung and low temperature cooling; it is tabulated on a uniform grid in log T and interpolated during the simulation to compute radiative losses.

For T < 10^4.2 K, we use the prescription from Gaspari et al. (2017) with a volumetric cooling rate defined as n_H²Λ_cold

$$\begin{array}{cc}\hfill {\mathrm{\Lambda }}_{\mathrm{cold}}& =2\times {10}^{-19}exp[-1.184\times {10}^5/(T+{10}^3)]\hfill \\ \hfill & +2.8\times {10}^{-27}\sqrt{T}exp[-92/T],\hfill \end{array}$$(13)

which incorporates crucial low temperature processes such as atomic line cooling, rotovibrational line cooling, and molecular collisions with dust grains. This implementation is crucial to realistically generate cold (10–100 K) gas.

We integrated the cooling source term, ℒ, with an exact implicit scheme (Townsend 2009; Gaspari et al. 2012), which provides an analytically exact update of the thermal energy over each timestep. This method is unconditionally stable, preserves positive internal energies, and eliminates the stringent timestep constraints normally imposed by stiff cooling functions, making it well suited for simulations where gas can cool rapidly over many orders of magnitude in temperature.

2.5. Turbulence driving

No explicit AGN feedback through jets and/or winds is included in this setup, as our primary focus is on the dynamics of the pure gas inflow transitioning from macro- to micro-scales. Nevertheless, we include gentle, large-scale AGN-motivated turbulence to capture the dynamical impact of SMBH activity on the surrounding medium. This approach provides a controlled framework for isolating inflow and cooling processes in the central regions of group-scale halos.

To sustain turbulence, we drive the gas with a stochastic acceleration field constructed from a superposition of random Fourier modes and evolved as an Ornstein–Uhlenbeck (OU) process (Uhlenbeck & Ornstein 1930; Schmidt et al. 2009; Grete et al. 2018, 2025), following the original CCA setup (Gaspari et al. 2013, 2017), in which continuous, predominantly subsonic stirring of a stratified atmosphere seeds a turbulent cascade and fosters nonlinear multiphase condensation.

The forcing peaks at a characteristic dimensionless mode number n_peak, corresponding to an injection scale L_inj = L_box/n_peak. The associated physical wavenumber is k_peak = 2πn_peak/L_box. The OU process is characterised by a correlation time t_corr, while the solenoidal versus compressive composition of the forcing is controlled by a weighting parameter ζ, such that ζ = 1 yields purely solenoidal forcing and ζ = 0 purely compressive forcing.

In Fourier space, the OU update of each individual acceleration mode reads

$$\begin{array}{c}\hfill {\widehat{a}}_i(\mathit{k},t+\mathrm{\Delta }t)=c_{\mathrm{drift}}{\widehat{a}}_i(\mathit{k},t)+\sqrt{1-c_{\mathrm{drift}}^2}{P}_{\mathrm{a}}(k){\mathcal{P}}_{\mathit{ij}}{\mathcal{N}}_j,\end{array}$$(14)

where k is the Fourier wavevector, c_drift = exp(−Δt/t_corr) sets the temporal correlation, 𝒩_j are complex random numbers with zero mean and modulus smaller than unity, and the target spectrum peaks at k_peak with

$$\begin{array}{c}\hfill P_{\mathrm{a}}(k)={\stackrel{\sim }{k}}^2(2-{\stackrel{\sim }{k}}^2)\mathrm{\Theta }(2-{\stackrel{\sim }{k}}^2),\end{array}$$(15)

where $\stackrel{\sim }{k}\equiv k/k_{\mathrm{peak}}$ and Θ denotes the Heaviside step function, which confines the forcing to a narrow shell around k_peak.

The solenoidal or compressive nature of the forcing is imposed through the projection tensor which implements the Helmholtz decomposition

$$\begin{array}{c}\hfill {\mathcal{P}}_{\mathit{ij}}=\zeta {\delta }_{\mathit{ij}}+(1-2\zeta )\frac{k_i{k}_j}{{|\mathit{k}|}^2},\end{array}$$(16)

where δ_ij is the Kronecker delta. Eq. (14) describes the OU evolution of each individual Fourier coefficient. The real-space turbulent acceleration field is then constructed by summing a finite number N_modes of such modes, whose wave-vectors k_m are randomly drawn within a narrow shell around k_peak:

$$\begin{array}{c}\hfill \mathit{a}(\mathit{x},t)=\mathfrak{R}\left[{\displaystyle \sum _{m=1}^{N_{\mathrm{modes}}}}\widehat{a}({\mathit{k}}_m,t)e^{i{\mathit{k}}_m·\mathit{x}}\right],\end{array}$$(17)

where ℜ denotes the real part of the complex Fourier sum. The turbulent forcing term entering the hydrodynamic Equations (2 and 3) is defined as f_turb(x, t)≡a(x, t), where a is the real-space stochastic acceleration field obtained from the Fourier-space OU driving described above. The overall forcing strength is controlled by the prescribed root-mean-square (RMS) acceleration amplitude,

$$\begin{array}{c}\hfill a_{\mathrm{rms}}\equiv {\langle {|\mathit{a}(\mathit{x},t)|}^2\rangle }^{1/2},\end{array}$$(18)

where angle brackets denote a volume average over the computational domain. At each timestep, the acceleration field is renormalised to match the prescribed value of a_rms and to ensure zero net momentum injection.

This method injects kinetic energy at large scales in a controlled and statistically stationary way. The resulting turbulent velocity dispersion and energy injection rate are determined by the combination of a_rms, L_inj, and t_corr. The parameters used for our simulations are listed in Table 3. Energy is injected solenoidally at an intermediate spatial scale, L_inj = 12.5 kpc, which is neither purely small-scale nor fully large-scale. As such, this choice is not uniquely associated with a single physical mechanism, but is broadly consistent with turbulence driven by AGN feedback, sloshing and/or minor mergers (e.g. Gaspari et al. 2012; Vazza et al. 2012). This intermediate-scale forcing initiates a turbulent cascade in which kinetic energy is transferred non-linearly from the injection scale to smaller scales, where it is eventually dissipated. This approach allows the system to self-consistently develop a turbulent spectrum characteristic of galaxy groups and clusters.

2.6. Simulation suite

All simulations were performed with the ATHENAPK code using SMR. The simulation box has a size of L_box = 50 kpc, with 12 nested levels of refinement reaching a finest resolution of ≃0.1 pc within the inner region. All the runs start from identical initial density and temperature profiles (see Figure 2) representing a static hot gaseous halo in hydrostatic equilibrium.

The simulation suite is summarised in Table 2. We progressively add physical processes to isolate their effects on the gas dynamics. The bondi run is obtained with an adiabatic simulation, with both driven turbulence and radiative cooling not active, and we use it as a benchmark to compare to the other simulations. The turb_low and turb_high runs include stochastic solenoidal driven turbulence with root-mean-square accelerations of a_rms = 6.2 × 10⁻⁹ cm s⁻² and a_rms ≃ 1.6 × 10⁻⁸ cm s⁻², but without radiative cooling. These two cases isolate the dynamical impact of different levels of subsonic turbulent stirring in a stratified hot atmosphere. Furthermore, we carry out two main simulations to investigate how the strength of turbulence influences the development of CCA. Both runs include radiative cooling but differ in their level of turbulence during the simulation time. We test two extremal weather regimes: a weak turbulence case (cca_low) and a strong turbulence case (cca_high), characterised by typical velocity dispersions of σ_v ≃ 60 − 90 km s⁻¹ and σ_v ≃ 210 − 230 km s⁻¹, respectively. Radiative cooling is activated after an initial evolution of 35 Myr, allowing turbulence to evolve towards our target velocity dispersions before condensation begins.

In our simulations, the sound speed of the IGrM in the central region is c_s ≈ 500 − 600 km s⁻¹, corresponding to average turbulence (3D) Mach numbers of ℳ ∼ 0.15 in cca_low and ℳ ∼ 0.4 in cca_high. Both cases therefore exhibit subsonic turbulence but differ by a factor of ∼3, bracketing direct and indirect observational constraints (e.g. XRISM and fluctuations power spectra; see Section 1). Note that galaxy groups can reach higher Mach numbers due to the lower sound speed or X-ray halo temperature (e.g. Hofmann et al. 2016). These bracketing setups provide a controlled framework to study how turbulence and cooling jointly regulate multiphase condensation and AGN feeding under significantly different weather conditions. Throughout the paper, we adopt a ‘BH weather’² terminology to qualitatively label different thermodynamical and dynamical states of the IGrM, using ‘sunny’, ‘rainy’, and ‘stormy’ conditions (see Figure 7 for a quantitative illustration of these regimes and Figure 13 for the conceptual framework) to denote, respectively, hot turbulence-dominated atmospheres and two different states of CCA.

3. Results

We now present the key results of our analysis, focusing on how turbulence affects gas condensation and accretion onto the SMBH. Our aim is to investigate how multiphase structure and inflow dynamics differ between the simulations with strong (cca_high) and weak (cca_low) turbulence. Although we treat them as distinct setups, cca_high and cca_low can also be interpreted as two possible evolutionary stages of the same system (see Section 4). Because stronger stirring enhances mixing and provides additional non-thermal support, the hot atmosphere in cca_high requires a longer time to cool and condense than in cca_low, so the two simulations evolve on intrinsically different physical timescales. For this reason, we normalise the time by the raining time t_rain, defined as the moment when the first cold (T ≈ 10⁴ K) raining gas appears. This normalisation allows for a direct comparison of their evolution relative to their respective raining cycles. t_rain is measured starting from the moment when radiative cooling is switched on. Each simulation is evolved up to t/t_rain = 3, corresponding to t_rain = 30 Myr for cca_high and t_rain = 7 Myr for cca_low (tied to the very inner cooling times).

For each analysis, we divided the gas into five thermal phases corresponding to distinct temperature ranges and observational bands: cold molecular (radio), cold (optical), warm (UV), and hot gas further separated into soft X-ray and hard X-ray components. The temperature boundaries for each phase are listed in Table 4. To investigate how the multiphase components form and evolve across scale (cf. the BLACKHOLEWEATHER diagram in Gaspari et al. 2020), we partitioned the computational domain into four radial regimes: micro (r ≤ 0.1 kpc), meso (0.1 < r ≤ 1 kpc), inner macro (1 < r ≤ 10 kpc), and outer macro (r > 10 kpc). In this context, the intermediate meso-scale (parsecs to kiloparsecs) remains one of the least explored regimes in current numerical studies, despite being the critical bridge between halo-scale condensation in cosmological simulations and the inner accretion flows resolved in general-relativistic magneto-hydrodynamical (GR-MHD) black hole simulations. This definition also aligns naturally with our radially concentric SMR design, allowing scale-by-scale diagnostics at a nearly uniform effective resolution within each regime.

The analysis of the results is organised as follows. In Section 3.1, we examine the morphology of cold clumps and filamentary structures. In Section 3.2, we study the SMBH accretion rate variability in the two turbulence regimes. In Section 3.3, we present radial profiles of the key thermodynamic variables and discuss the implied thermal and dynamical balance of the atmosphere. We then characterise the emergence of multiphase gas using density-temperature phase diagrams (Section 3.4) and phase mass evolution (Section 3.5). Finally, in Section 3.6, we analyse the probability density functions of the condensed multiphase gas.

3.1. Filamentary structure

Turbulence governs the chaotic motions of the hot gas and seeds the density and temperature fluctuations from which the multiphase medium condenses via nonlinear thermal instability (Gaspari et al. 2013, 2017). Its strength plays a key role in regulating how gas cools, fragments, and is ultimately accreted onto the SMBH. We begin our analysis by comparing the morphological structure of the filaments and clumps that emerge in the two turbulence regimes cca_high and cca_low.

Figure 4 shows gas density slices at different times and physical scales in the cca_high run, illustrating how the gas transitions from an initially smooth medium to a fully developed multiphase structure. The first panel in each row provides a large-scale view of the hot halo showing the full 50 kpc simulation box. At the beginning of the simulation, the hot gas forms a coherent and thermally stable structure, with only mild fluctuations induced by turbulence. As turbulence develops, it stirs the hot atmosphere, generating elongated density fluctuations aligned with the large-scale shear of the flow. Nonlinear compressions induced by turbulence trigger localised catastrophic cooling, generating cold gas filaments and clumps around t ≈ 30 Myr. By t/t_rain ≈ 1 (see Figure 4, first row), the first filamentary condensations appear along regions of converging motions. The turbulent compression locally enhances the density and cooling rate, allowing cold gas to form despite the hot background. The filaments are often aligned with the turbulent eddies and concentrated in the central ≈10 kpc. At later times (t/t_rain ≃ 2.5 − 3), the system shows fragmentation of the cold and dense gas driven by the turbulent cascade, which transfers kinetic energy to smaller scales and increases the velocity dispersion in the core to σ_v ≈ 230 km s⁻¹. Multiple filaments interact, merge, and break apart under the action of turbulence, cooling, and gravity. The dense structures develop internal denser sub-filaments and knots, forming a tangled network surrounding the central region.

Fig. 4. Gas density slices of the central 50 kpc region in the cca_high simulation. Rows show the time evolution from t/train = 1.5 (top) to t/train = 3 (bottom), while within each row the panels from left to right present a progressive zoom-in from the halo scale (macro) to the filamentary condensation region (meso) and finally to the innermost clumps (micro). The sequence reveals the characteristic morphology of CCA: a complex network of cold, dense filaments condensing out of the turbulent hot halo and converging towards the galactic centre. Panels at different spatial scales adopt independent logarithmic colour ranges to optimise contrast and highlight the gas substructures; colour scales are thus not directly comparable across panels. The top right corner of each panel in the last column shows a zoom-in of the innermost rotating structure.

From left to right (Figure 4), the columns progressively zoom in from the halo scale (tens of kpc, i.e. macro-scale) to the filamentary condensation region (≈1 kpc, i.e. meso-scale) and finally to the innermost clumps (central parsec, i.e. micro-scale), capturing the end state of the condensation cascade. Towards the galactic centre, this clumpy multiphase complex extends all the way to the SMBH, where thin cold filaments transport gas inward as clouds rain in from kiloparsec scales.

In the cca_low simulation (Figure 5), the central ∼1 kpc region is shown at t/t_rain = 1.5 and t/t_rain = 3, zooming in towards the centre of the halo, exhibiting distinct morphology and distribution compared to cca_high. The initially smooth inflows condense after t = 7 Myr, earlier than in the cca_high case due to the lower turbulent heating-forming a multiphase structure that remains predominantly filamentary rather than clumpy. The resulting morphology is irregular yet more coherent than in cca_high, with cold gas circulating in a less disordered pattern. Compared to cca_high, the cold phase is more centrally concentrated, being largely confined within the innermost 100 pc.

Fig. 5. Gas density slices of the central region in the cca_low simulation. Rows show the time evolution from t/train = 1.5 (top) to t/train = 3 (bottom), while within each row the panels from left to right present a zoom-in from the inner kiloparsec to 15 pc around the SMBH. At both epochs, the system hosts a dense, clumpy cold core on parsec scales, continuously fed by smoother filamentary inflows from larger radii. The top right corner of each panel in the last column shows a zoom-in of the innermost rotating structure. Cooling dominates over turbulent mixing, allowing cold gas to accumulate and recycle within a long-lived multiphase region where clouds and filaments move along intersecting, colliding, and shearing trajectories. The morphology remains irregular and fragmented, illustrating the CCA mode also in a weakly turbulent galactic atmosphere.

In both simulations, a clumpy rotating structure (see the inset panels in the last panels of Figures 4 and 5) emerges spontaneously at the centre despite the initially non-rotating conditions. Indeed, turbulent eddies induce local vorticity and stochastically broaden the angular momentum distribution in the hot phase (see also Gaspari et al. 2015); as gas cools and condenses, anisotropic inflows transfer a small residual angular momentum to the collapsing material. Since multiphase inflows do not perfectly cancel out, the remaining angular momentum accumulates in the centre, forming a rotating structure of characteristic radius ≈10 pc in cca_high and ≈5 pc in cca_low composed mainly of cold and molecular gas but with a substantial fraction of warm gas. Given our resolution and sink treatment, we do not interpret this feature as a fully resolved, rotationally supported accretion disc; rather, it is best viewed as a torus-like, clumpy rotating complex at the micro-scale, potentially analogous to the commonly observed clumpy AGN torus (and/or the outermost reservoir feeding the unresolved inner disc). This structure is continuously perturbed by infalling filaments; therefore, the inflow remains chaotic and highly time-variable despite the presence of net rotation. These simulations show that the meso-scale (parsecs to kiloparsecs) is not a passive transition region, but an active regime where turbulence and radiative cooling together shape the filamentary network. In this range, cold filaments form, fragment, and interact, establishing the geometry and radial extent of the cold phase before it is accreted towards the central region.

Figure 6 illustrates the contrasts in the thermodynamic properties by showing slices of temperature (left panel), pressure (middle panel), and velocity magnitude (right panel) across the group core in the cca_high (top row) and cca_low (bottom row) simulations. The temperature map reveals the coexistence of hot ambient gas (T ≈ 10⁶ − 10⁷ K) and cooler condensed gas (T ≈ 10 − 100 K), with a clear multiphase structure within the filaments and clumps, where warm gas (≈10⁴ K) surrounds cold molecular cores. The diffuse halo gas instead remains in the 10⁶ − 10⁷ K range. The pressure slice reveals thermal pressure contrasts between the cold structures and the volume-filling hot gas, especially in the cca_high. Cold structures correspond to regions of lower thermal pressure (≈10⁻¹³ − 10⁻¹² dyne cm⁻²) compared to the surrounding hot medium (≈10⁻¹⁰ − 10⁻⁹ dyne cm⁻²), showing a non-isobaric condensation. Finally, the velocity field shows the strong, chaotic motions that continuously reshape the filamentary network. The velocity maps reveal complex kinematics driven by turbulent motions, spanning velocities from ∼10 to 10³ km s⁻¹, with most cold structures being transported and accreted towards the central region. The slices emphasise the contrasting filament distributions, with cca_high showing extended structures across the central kpc and cca_low retaining cold gas confined to the innermost tens of parsecs. An in-depth quantitative analysis of cold clouds and filaments, including their properties and spatial distribution, will be presented in a separate, upcoming work.

Fig. 6. Slices of the simulated galaxy group core for the cca_high (top row) and cca_low (bottom row) simulations at t/train = 3, showing temperature (left), pressure (middle), and velocity magnitude (right). The maps highlight the multiphase structure of the IGrM and the role of turbulence in driving thermodynamic fluctuations in the hot atmosphere, seeding local thermal instability and the condensation of cold gas. Each panel in the top row shows an 8 kpc × 8 kpc region, while each panel in the bottom row shows the central 0.2 kpc × 0.2 kpc region.

3.2. Accretion rates

To highlight the differences between these different CCA regimes we analyse the time evolution of the accretion rate in our simulations. Here we provide a brief overview, while a more detailed analysis of gas accretion is presented in the companion paper (B26b). Ṁ is computed as the mass accreted within the sink region (r < 0.4 pc) at each simulation timestep. We express the accretion rate in units of the Bondi rate Ṁ_B, which provides a natural reference scale for hot mode accretion. Figure 7 shows the normalised accretion rates Ṁ/Ṁ_B as a function of normalised time t/t_rain for the cca_high (blue, stormy weather) and cca_low (cyan, rainy weather) simulations, compared with the adiabatic Bondi run (bondi, black), the purely turbulent control runs (sunny weather) turb_high (red) and turb_low (yellow), and the cooling-only simulation cool (purple).

Fig. 7. SMBH accretion rate Ṁ normalised to the Bondi accretion rate ṀB as a function of normalised time t/train for the cca_high (blue line, stormy weather) and the cca_low (cyan line, rainy weather) simulations, compared with the turbulence simulations (sunny weather) turb_high (red line) and turb_low (yellow line), with the radiative cooling simulation cool (purple line) and with the idealised adiabatic simulation bondi (black line). The CCA phenomenon produces a highly variable accretion rate due to its chaotic nature.

The classical Bondi accretion rate (Bondi 1952) is given by

$$\begin{array}{c}\hfill {\dot{M}}_{\mathrm{B}}=\lambda 4\pi {(G{M}_{•})}^2\frac{{\rho }_{\infty }}{c_{\mathrm{s},\infty }^3},\end{array}$$(19)

where λ = (1/2)^{(γ + 1)/2(γ − 1)}[(5 − 3γ)/4]^{−(5 − 3γ)/2(γ − 1)} is a normalisation factor of order unity, γ is the gas adiabatic index, M_• is the black hole mass, ρ is the gas density and c_s is the gas sound speed (the infinity symbol denotes very large radii). Although this expression is widely employed in numerical and observational studies, often even when the Bondi radius is unresolved (e.g. Di Matteo et al. 2005; Cattaneo & Teyssier 2007; Booth & Schaye 2009; Yang et al. 2012), it assumes a steady, homogeneous, adiabatic, and spherically symmetric inflow, which is not found in realistic galactic environments. Therefore, the adiabatic bondi run provides a stable reference for the hot-mode flow but cannot be used to compute a realistic accretion rate. The BH accretion rate has a nearly constant value of Ṁ_• ≃ 6 × 10⁻⁴ M_⊙ yr⁻¹, in close agreement with the analytic Bondi prediction based on the ambient hot-gas density and temperature (Eq. (19)).

The impact of turbulence in the absence of cooling is isolated in the purely turbulent runs. Both turb_low and turb_high remain firmly in the hot accretion regime, but exhibit systematically sub-Bondi accretion rates. In the turb_low run, characterised by moderate subsonic stirring (ℳ ≃ 0.1 − 0.2), the accretion rate fluctuates around Ṁ_• ≈ 0.6 Ṁ_B. Increasing the turbulent velocity dispersion in turb_high (ℳ ≃ 0.4) further suppresses the inflow, yielding Ṁ_• ≈ 0.5 Ṁ_B and a gradual decline at late times. In these runs, turbulent motions generate vorticity and small rotating eddies that support gas against the accretion, decreasing direct inflow towards the SMBH. In addition, bulk turbulent motions introduce a non-zero relative velocity between the gas and the black hole, reducing the effective accretion rate in a Bondi–Hoyle–like manner (Bondi & Hoyle 1944).

The cool simulation shows a strong increase in the accretion rate, from ≈ 1 Ṁ_B up to ≈ 10⁴ Ṁ_B. This occurs because gas in the central region cools efficiently and accumulates towards the centre, producing accretion rates several orders of magnitude above the Bondi value. This configuration is not realistic, as it neglects perturbations in the medium, which are expected to be present in real astrophysical systems.

Including radiative cooling in a turbulent medium fundamentally changes the accretion regime. Once multiphase condensation sets in, the SMBH transitions from a turbulence-dominated, sunny atmosphere to a CCA mode, corresponding to rainy and stormy weather conditions. The formation of cold filaments and clumps, followed by their inelastic interactions and loss of angular momentum, leads to accretion rates that are strongly super-Bondi. In both cca_high and cca_low, the accretion rate rises sharply after t/t_rain ∼ 1, reaching values one to two orders of magnitude above the Bondi rate. Typical levels of Ṁ_• ∼ 10−100 Ṁ_B are sustained for extended periods, accompanied by strong variability. The accretion becomes highly intermittent, with rapid fluctuations driven by chaotic cloud–cloud collisions and the sporadic infall of cold filaments. Despite their markedly different filamentary morphologies, the cca_high and cca_low runs display qualitatively similar accretion histories once CCA is established, indicating that the presence of multiphase gas dominates the regulation of SMBH fueling. On average, both runs remain strongly super-Bondi, highlighting the inefficiency of Bondi-based prescriptions in capturing accretion in cooling, multiphase halos.

3.3. Thermodynamic variables radial profiles

To investigate the thermodynamic stratification and time evolution of the multiphase IGrM, we computed mass-weighted radial profiles of density, temperature, pressure, and velocity magnitude for the two fiducial simulations, cca_high and cca_low (Figure 8). Gas is separated into five temperature phases (see Table 4): hot hard X-ray (red), hot soft X-ray (orange), warm (violet), cold (blue) and molecular (cyan) phases. Each profile represents a mass-weighted average within 3D spherical shells centred on the centre of the simulation, extending out to ≈25 kpc. The colour intensity encodes the temporal evolution, from early (faint, t/t_rain = 0) to late (bright, t/t_rain = 3) stages. The overlaid ACCEPT profiles (pink circles) show that the simulated hot X-ray phase is broadly consistent with observed IGrM density, temperature, and pressure profiles on kpc scales³.

Fig. 8. Mass-weighted radial profiles of gas density, temperature, thermal pressure, and velocity magnitude for simulations cca_high (top row) and cca_low (bottom row). The gas is divided into five temperature phases: hot, hard X-ray (red) and soft X-ray (orange), warm (violet), cold (blue) and molecular (cyan, see Table 4). Increasing colour intensity indicates later stages of evolution. Profiles are computed in spherical shells around the group centre out to 25 kpc, showing how the thermodynamic and kinematic structure of the multiphase IGrM evolves in the two CCA regimes. Pink circles show observational density, projected temperature, and pressure profiles from the ACCEPT database (Cavagnolo et al. 2009), where applicable.

In the cca_high run, the hot plasma fills the entire simulation volume, with densities of 10⁻²⁵ − 10⁻²⁴ g cm⁻³ at r ∼ 1 kpc and temperatures of several 10⁷ K. The profiles also reveal a non-negligible amount of gas below 10⁶ K, confined within the innermost ∼10 kpc. As expected, the gas density increases towards colder phases, from the hard X-ray component down to the molecular phase, which reaches 10⁻¹⁹ − 10⁻¹⁸ g cm⁻³ in the central 100 pc. This central region exhibits the strongest variability in the density profiles, as clumps and filaments can be disrupted from their mutual interactions or turbulence effects, while the outer halo remains comparatively stable. The temperature of the hard X-ray phase remains nearly flat beyond the central 100 pc and slightly increases towards the centre. The hot soft X-ray, warm, and cold components show mildly increasing temperatures with radius. The molecular phase instead exhibits small-scale fluctuations with radius but no clear systematic trend, remaining roughly constant on average.

The pressure peaks in the centre and decreases with radius in all gas phases. The hot hard and soft X-ray components show a pronounced central excess of ≈10⁻⁶ dyne cm⁻². The colder phases reach central pressures of 10⁻⁷ − 10⁻⁸ dyne cm⁻² and start to separate beyond ∼10 pc: the warm gas maintains the highest pressure at ≈10⁻¹¹ dyne cm⁻², followed by the cold phase, while the molecular cold component drops to nearly 10⁻¹² dyne cm⁻². This radial decline reflects gravitational stratification and is consistent with approximate local pressure balance between the different phases. Differences between the temperature phases pressure profiles are expected, since the warm, cold, and molecular gas tend to occupy different parts of the multiphase structures, such as mixing layers or filaments centres.

Although the molecular, cold and warm components seem to be thermally under-pressured relative to the ambient hot medium, the hot and cold phases are in approximate pressure equilibrium once non-thermal turbulent pressure (e.g. Khatri & Gaspari 2016) is included. We defined a total pressure,

$$\begin{array}{c}\hfill P_{\mathrm{tot}}\equiv P_{\mathrm{th}}+P_{\mathrm{nt}}\simeq n{k}_{\mathrm{B}}T+\frac{1}{3}\rho {\sigma }_v^2,\end{array}$$(20)

where k_B is the Boltzmann constant and σ_v is the gas velocity dispersion. We can simply derive this pressure equilibrium for the molecular component in cca_high, as it is the final outcome of the condensation cascade. For the hot phase, adopting typical values of the gas in the central kiloparsec (n_hot ≈ 0.1 cm⁻³ and T_hot ≈ 10⁷ K, see Section 3.4) yields

$$\begin{array}{c}\hfill P_{\mathrm{th},\mathrm{hot}}=n_{\mathrm{hot}}{k}_{\mathrm{B}}{T}_{\mathrm{hot}}\approx 1.4\times {10}^{-10}{\mathrm{dyne}\mathrm{cm}}^{-2}.\end{array}$$(21)

The corresponding non-thermal term for σ_hot ≈ 200 km s⁻¹ is

$$\begin{array}{c}\hfill P_{\mathrm{nt},\mathrm{hot}}=\frac{1}{3}{\rho }_{\mathrm{hot}}{\sigma }_{\mathrm{hot}}^2\approx 1.3\times {10}^{-11}{\mathrm{dyne}\mathrm{cm}}^{-2},\end{array}$$(22)

so that P_tot, hot ≈ 1.5 × 10⁻¹⁰ dyne cm⁻². In contrast, the cold molecular gas has a much smaller thermal pressure. Using n ≈ 100 cm⁻³ and T ≈ 100 K, we obtain P_th, cold ≈ 10⁻¹² − 10⁻¹³ dyne cm⁻² (see also the second panel in Figure 6). The high densities of the molecular gas imply that even modest internal motions provide substantial non-thermal support. Adopting μ ≃ 2.3 for molecular gas and σ_cold ≈ 10 km s⁻¹, we obtain

$$\begin{array}{c}\hfill P_{\mathrm{nt},\mathrm{cold}}=\frac{1}{3}{\rho }_{\mathrm{cold}}{\sigma }_{\mathrm{cold}}^2\approx 1.3\times {10}^{-10}{\mathrm{dyne}\mathrm{cm}}^{-2}.\end{array}$$(23)

Therefore, despite being thermally under-pressured, the cold molecular phase satisfies P_tot, cold ≈ P_tot, hot, implying that the multiphase medium is close to pressure balance when considering the non-thermal pressure component (magnetic fields could further increase this contribution). More broadly, recent SZ observations of AGN-inflated X-ray cavities also suggest that non-thermal pressure components can provide an important contribution to the pressure budget in hot halo atmospheres (e.g. Abdulla et al. 2019; Orlowski-Scherer et al. 2022).

Importantly, the degree of thermal under-pressure depends on turbulence strength. In the cca_high run, strong non-thermal support allows cold clouds to remain dynamically confined while being significantly under-pressured in thermal terms. In contrast, in the cca_low case the reduced non-thermal support requires a higher thermal compression of the cold gas, bringing the system closer to an almost isobaric configuration (Figure 6). The weaker turbulent motions in the hot phase also limit the spatial coherence of cooling compressions, leading to the formation of smaller and more compact cold clouds compared to cca_high.

The velocity field is dominated by subsonic turbulent motions on kiloparsec scales, with characteristic speeds of ≈300 km s⁻¹ in the hard and soft X-ray phases and ≈200 km s⁻¹ in the cold and warm gas, with a mild decline towards colder phases, reflecting the decoupling of cooler gas from the hot turbulent eddies and the partial dissipation of turbulent motions during condensation. Towards the innermost parsecs, the flow accelerates in the deep gravitational potential of the SMBH, and the typical velocities rise to ≈10³ km s⁻¹ within the central sub-parsec region.

The cca_low run shows broadly similar radial trends but with a more centrally concentrated cold and warm component. The molecular phase reaches high densities in the inner few hundred parsecs, but its density declines more rapidly beyond r ∼ 0.1 kpc and it is effectively confined within ∼1 kpc. Its temperature remains in the tens to ≈100 K range and the corresponding pressure profile lies systematically below that of cca_high beyond r ∼ 0.01 kpc, reflecting the smaller radial extent of the coldest gas. The warm and cold phases are likewise more compact, typically extending only out to ∼2 − 3 kpc, while the hot atmosphere dominates at larger radii. Overall, cca_low still develops a multiphase core, but the cold and warm gas are less extended and more centrally concentrated than in cca_high, consistent with weaker turbulent transport.

3.4. Phase diagrams

To characterise the multiphase nature of the IGrM (as recurrently found in observations, e.g. Sun et al. 2009; Temi et al. 2018; Eckert et al. 2021; Olivares et al. 2022; Ubertosi et al. 2025; Temi et al. 2026), we start by comparing the phase diagrams obtained from the different simulations. Figure 9 shows the phase diagrams of the benchmark runs in which radiative cooling is deactivated (turb_high and turb_low) and of the run in which turbulent driving is switched off (cool). The oblique grey dashed lines indicate reference thermodynamic paths. Lines with negative slope represent isobaric evolution (T ∝ n⁻¹), corresponding to gas cooling or heating while remaining in pressure equilibrium with its surroundings. Lines with positive slope show adiabatic evolution (T ∝ n^2/3), associated with compressions and expansions driven by turbulent motions. These guides help illustrate the dominant thermodynamic behaviour of the gas.

Fig. 9. Mass-weighted density-temperature phase diagrams for the cool (blue), turb_high (red) and turb_low (orange) simulations. Dashed oblique lines indicate adiabatic (T ∝ nγ − 1 = n2/3) and isobaric (T ∝ n−1) relations.

In the turbulence-only simulations, the ICs form a stripe at low density. With time, a high-density tail emerges as gas is adiabatically compressed towards the centre, similarly to Bondi accretion. Meanwhile, turbulence causes the initial stripe to widen through repeated adiabatic compressions and expansions. turb_high maintains gas at systematically higher temperatures, while reaching lower maximum densities with respect to turb_low. The stronger turbulent motions disrupt the gas in the centre, preventing it from accumulating. The gas remains confined to a hot, low-density branch. As a result, these runs do not produce a realistic multiphase medium, but instead maintain a nearly single-phase atmosphere.

In the cool simulation, the opposite behaviour is observed. Gas undergoes radiative cooling, forming a cooling sphere expanding from the centre in which material rapidly transitions from the hot phase to very cold temperatures, from the stable T ∼ 10⁴ phase down to 10 K. As a result, cold gas accumulates near the centre and fails to reproduce the extended multiphase structure observed in galaxy groups. A realistic multiphase IGrM does not arise from turbulence or cooling alone, but from their combined action.

Figure 10 shows the temporal and spatial evolution of the gas density–temperature phase diagram in the cca_high and cca_low fiducial simulations. In the first two rows, the panels illustrate how the phase structure develops in the full simulation volume at increasing times, expressed in units of t_rain (t/t_rain = 0,  1,  2, and 3). The last two rows show the gas phase structure at t/t_rain = 3, decomposed into the four radial ranges: micro-scale (r ≤ 0.1 kpc), meso-scale (0.1 < r ≤ 1 kpc), an inner macro-scale (1 < r ≤ 10 kpc), and an outer macro-scale (r > 10 kpc). In each panel we show dashed oblique lines indicating adiabatic (T ∝ n^2/3) and isobaric (T ∝ n⁻¹) evolution.

Fig. 10. Mass-weighted density–temperature phase diagrams for the cca_high (rows 1 and 3) and cca_low (rows 2 and 4) simulations. The first two rows show the time evolution at t/train = 0.1,  1,  2,  3, while the last two rows show the corresponding distributions at t/train = 3 in four radial ranges: 0 < r ≤ 0.1 kpc (micro-scale), 0.1 < r ≤ 1 kpc (meso-scale), 1 < r ≤ 10 kpc (inner macro-scale), and r > 10 kpc (outer macro-scale). Colours denote the total gas mass in each (n, T) bin, from dark brown (low mass – 10−5 M⊙) to yellow (high mass – 105 M⊙). Phases (a-f) are classified by temperature (hot: T ≳ 106 K; warm: 104 < T < 106 K; cold: T < 104 K) and density (diffuse vs. dense, with thresholds at n = 1 cm−3 for hot/warm and n = 100 cm−3 for cold). Dashed oblique lines indicate adiabatic (T ∝ n2/3) and isobaric (T ∝ n−1) relations in each panel.

At the beginning of both simulations (t/t_rain = 0.1), the gas lies almost entirely in the hot, volume-filling phase with T ≳ 10⁶ K and densities 10⁻² < n < 1 cm⁻³, as seen in the ICs profile (Figure 2). The atmosphere is still monophase and in approximate hydrostatic equilibrium. At this stage, the cooling time is much longer than the dynamical time; condensation is negligible and most of the mass remains in the hot turbulent halo. By t/t_rain = 1, the picture starts to change. The locus in the phase diagram broadens significantly: while a large fraction of the gas remains on the hot branch, a tail of material extends towards higher densities (up to n ≈ 10² cm⁻³) and lower temperatures. In cca_high, the gas spans a broader range of thermodynamical conditions, reaching temperatures as low as 10² K, whereas in cca_low, the cold component is confined to the 10⁴ K phase, indicating a slower progression towards the cold molecular regime. The phase plot at this stage shows gas populating mainly intermediate densities (n ∼ 1 − 10² cm⁻³) and temperatures (T ∼ 10^4 − 5 K), forming a continuous bridge between the hot ambient medium and the emerging cold phase. This is the first signature of thermal instability triggered by the interplay of radiative cooling and turbulent fluctuations. Such overdensities can grow non-linearly when the cooling time becomes sufficiently short relative to the free-fall time, commonly expressed by t_cool/t_ff ≲ 10 (e.g. Field 1965; Pizzolato & Soker 2005; Gaspari et al. 2012; McCourt et al. 2012; Sharma et al. 2012), thereby forming the chaotic cold rain. Overdense regions, compressed by turbulent motions, cool more efficiently and move away from the hot phase locus.

At t/t_rain = 2, a distinct multiphase structure has developed. In addition to the hot background, there is now a well-populated branch extending down to T ≲ 10 K reaching densities up to n ∼ 10⁵ − 10⁶ cm⁻³. Dense clumps and filamentary structures condense at n ≳ 1 cm⁻³, becoming thermodynamically decoupled from the surrounding medium.

At t/t_rain = 3, a wider multiphase picture appears. In the cca_high run, the stronger turbulence produces a highly fragmented and extended phase-space distribution: the gas spans over ten orders of magnitude in density and almost eight in temperature. The cold component reaches down to T ≲ 10 K and n ≳ 10⁶ cm⁻³. Although the hot phase is still dominant, the IGrM is multiphase, signaling a fully developed regime of CCA. The broader extent of the distribution in the cca_high simulation is a direct consequence of the stronger turbulent motions. Turbulence continuously compresses and expands the gas, driving its evolution along the oblique adiabatic trajectories. These rapid compressions and rarefactions widen the distribution in both density and temperature, preventing the gas from remaining confined to a narrow distribution. As a result, the cca_high simulation exhibits a significantly broader phase-space structure than cca_low, where weaker turbulence produces more limited density contrasts and a more compact distribution. In contrast, the cca_low run remains more compact in phase space. Weaker turbulence limits both density contrasts and the nonlinear efficiency of cooling. Cold gas forms only in the centre (see Figure 5) partially, reaching minimum temperatures of T ∼ 10 K and densities up to n ≲ 10⁶ cm⁻³. Most of the mass stays in the hot or warm regime, with a smoother transition between phases and without the sharp multiphase separation seen in the high-turbulence case. In both simulations, and more prominently in cca_high, part of the gas evolves along isobaric tracks towards higher temperatures (T ≳ 10⁷ K), consistent with heating driven by turbulent dissipation.

The diagram reveals the presence of a rich multiphase medium, we highlight 6 different coexisting components in the forth panels: (a) The diffuse hot phase (T ≳ 10⁶ K, n < 1 cm⁻³) represents the volume-filling hot plasma in quasi-hydrostatic equilibrium. (b) The dense hot phase (T ≳ 10⁶ K, n > 1 cm⁻³) consists of gas locally compressed by turbulent stirring or shocks, where cooling is about to become effective. (c) The diffuse warm phase (10⁴ < T < 10⁶ K, n < 1 cm⁻³) traces the mixing interfaces between the hot atmosphere and condensed filaments. (d) The dense warm phase (10⁴ < T < 10⁶ K, n > 1 cm⁻³) marks thermally unstable gas rapidly cooling towards the cold regime. (e) The diffuse cold phase (T < 10⁴ K, n < 100 cm⁻³) forms the extended envelopes of filaments and clumps. (f) The dense cold phase (T < 10⁴ K, n > 100 cm⁻³) corresponds to the molecular gas that can eventually form stars (Ferrière 2001). The presence of ongoing star formation in active AGN environments has been detected also by recent JWST observations (e.g. Reefe et al. 2025). In our simulations, cold dense gas is present and could potentially form stars before reaching the central SMBH, provided that the local free-fall time is shorter than the Jeans time (e.g. t_ff ≈ 10 Myr > t_jeans ≈ 1 Myr). This scenario is supported by observations revealing star formation within cold filaments located several kiloparsec from the galaxy centre (Tremblay et al. 2016). While star formation is not explicitly included in the present work, we plan to investigate its role in future studies. The emergence of the cold and molecular phases reflects the thermodynamic transition driven by density fluctuations induced by turbulence and radiative cooling, leading to the formation of cold clumps and filaments that decouple from the hot halo and sink towards the SMBH enhancing accretion.

The spatial distribution of the gas changes markedly with radius (Figure 10, third and fourth rows). At the micro-scale, both simulations show the majority of the multiphase gas. This highlights the crucial importance of resolving these small physical scales to correctly study the accretion of gas at the centre of galaxies. Both very hot and very cold gas phases are present in this compact region.

At the meso-scale region the two simulations diverge. In the cca_high simulation, the enhanced turbulent motions preserve a complex network of cold filaments and clumps embedded within the hot atmosphere. In contrast, cca_low still hosts both warm and cold gas, but the IGrM is dominated by more discrete phases – a warm component with T ≈ 10⁴ K and n ≈ 1 − 10² cm⁻³, and a colder, denser component with n ≈ 10² cm⁻³. This suggests that weaker turbulence leads to less efficient mixing and fragmentation.

At the inner macro-scale, the cca_high simulation retains significant amounts of intermediate-temperature gas (10^4.5 − 10⁶ K) and cold gas down to 10 K, whereas the cca_low case is mainly composed by the hot phase with a small fraction of warm and 10⁴ K gas. Finally, at the outer macro-scale, both runs are dominated by a tenuous, hot medium (T ≳ 10⁶ K), but the cca_high simulation shows a non negligible amount of cold gas at T ≈ 10⁴ K, indicating that turbulence promotes the persistence of multiphase structures even in the outskirts, whereas cca_low shows no presence of cold gas. Overall, at the same normalised time, both simulations show a rich multiphase environment but with very different spatial distribution. cca_high develops a more violent fragmentation and a fully established multiphase structure until ≳10 kpc, while cca_low remains more dominated by the hot medium, with smaller and more stable condensations in the central kiloparsec.

3.5. Mass evolution

The differences between the two turbulence regimes become clear when examining the amount of gas in the different thermal phases. Figure 11 shows the gas mass of each phase as a function of time at the micro (r ≤ 0.1 kpc), meso (0.1 < r ≤ 1 kpc), and inner macro (1 < r ≤ 10 kpc) scales. The total average mass of each phase for the two simulations is shown in Table 5 together with the mass ratio. Overall the cca_high simulation contains larger masses of warm, cold, and molecular gas at all radii.

Fig. 11. Gas mass for cca_high (solid) and cca_low (dashed) simulations divided in molecular (cyan), cold (blue), warm (violet), soft (orange), and hard (red) X-ray gas and in different scales: macro, meso, inner macro and outer macro.

At the micro-scale (top panel), all the phases masses in cca_high vary strongly with time due to the stronger turbulence and the interaction between different clumps and filaments that can heat up or disrupt the cold gas. At t < t_rain the medium is dominated by the hot gas in both simulations (M_hot ≈ 10⁴ M_⊙). After t/t_rain = 1 warm, cold and molecular gas start to appear, with the molecular gas mass reaching M_mol ≈ 10⁶ M_⊙ at t/t_rain = 3. At the same time, the cca_low run shows a much smaller molecular mass, M_mol ≈ 5 × 10⁴ M_⊙, with a smoother and more gradual increase. A similar behaviour is found for the cold phase, although with lower masses of M_cold ≈ 10⁵ M_⊙ in cca_high and M_cold ≈ 10⁴ M_⊙ in cca_low. In contrast, the warm gas mass is comparable in the two runs, with values around M_warm ≈ 10³ M_⊙. Even at these small scales, differences between the two regimes are already present.

The meso-scale region (middle panel) marks the first appearance of warm and cold gas in cca_high, at t/t_rain = 1. By contrast, in cca_low the cold phase first forms at the micro-scale. This highlights a difference in where condensation initially develops, although in both simulations a large fraction of M_cold is later found in this region. The cold gas mass in cca_high (M_cold ≈ 10⁶ − 10⁷ M_⊙) remains systematically higher than the warm phase (M_warm ≈ 10⁵ − 10⁶ M_⊙) throughout the simulation, without a steady monotonic increase but again with an oscillatory trend showing the strong turbulence and variability of the medium. The molecular gas forms later, after the appearance of the warm and cold phases. This happens at all scales for both simulations and provides insight into the condensation cascade: gas first cools into the warm phase and then rapidly transitions into the colder, more stable phase. As a result, most of the warm gas is confined to the interfaces between hot and cold material within clumps and filaments, as also seen in Figure 6. The molecular gas mass in cca_high remains high, although slightly below the cold gas mass, almost reaching M_mol ≈ 10⁷ M_⊙ at t/t_rain = 3. In contrast, in cca_low the molecular component is much smaller, appearing only after t/t_rain ≈ 1.5 and reaching M_mol ≈ 10⁵ M_⊙ at t/t_rain = 3. The hot phase in cca_high is irregular over time (M_hot ≈ 10⁵ − 10⁷ M_⊙), whereas in cca_low it is almost constant (M_hot ≈ 10⁷ M_⊙).

At the inner macro-scale (bottom panel), the cca_low simulation shows almost no warm or cold gas during the early evolution, with these phases starting to grow only after t/t_rain ≈ 2. In cca_high, at larger distances from the centre, the phase masses display a smoother, less oscillatory evolution and increase steadily with time. As a consequence, the hot gas mass slowly decreases from ≈5 × 10⁹ M_⊙ to ≈10⁹ M_⊙ as gas cools and transitions into lower-temperature phases. The mass of cold gas in cca_high strongly increases over time passing from M_cold ≈ 10⁶ M_⊙ to 10⁹ M_⊙, the same happens to the warm phase which increases reaching M_warm ≈ 10⁸ M_⊙. The molecular phase is formed soon after also steadily increasing with slight oscillations arriving to 2 − 3 × 10⁷ M_⊙ by the end of the simulation.

Despite the cca_high simulation developing a substantially larger cold and molecular gas reservoir (Table 5), the SMBH accretion rate remains comparable, and in some phases slightly lower, to the one measured in the cca_low run (see Figure 7). This indicates that BH feeding is not directly set by the total mass of cold gas present, but rather by the efficiency with which cold structures can lose angular momentum and be transported towards the centre. This is consistent with high-resolution ALMA observations showing no clear correlation between ≲100 pc molecular gas mass and AGN activity tracers (Elford et al. 2024). In the stormy regime, a significant fraction of the cold gas remains stored in extended, rotationally supported and dynamically stirred structures, while only a limited portion effectively participates in accretion. Conversely, in the rainy regime, cold gas forms less abundantly but is accreted more efficiently once it condenses, leading to comparable inflow rates despite the much smaller cold gas reservoir.

3.6. Multiphase density structure across scales

In Figure 12 we present the time evolution of the mass-weighted gas density probability distribution functions (PDFs) for different thermal phases, measured in three radial shells spanning the micro, meso, and inner macro regimes. For each shell we select all gas cells within the specified radial range and bin them uniformly in log n. The PDFs are mass-weighted: in each density bin i we compute the mass fraction $p_i\equiv \frac{1}{M_{\mathrm{shell}}}{\displaystyle \sum _{j\in i}}{m}_j$, where m_j is the mass of cell j and M_shell = ∑_jm_j is the total gas mass in the shell. The corresponding differential PDF is d𝒫/dlog n ≃ p_i/Δlog n, with Δlog n the bin width. By construction, ∑_ip_i = 1 (equivalently ∫(d𝒫/dlog n) dlog n = 1), so the integral over any log n interval yields the corresponding mass fraction. We adopt 150 bins spanning log n = −3 to 6 and quantify each phase via the mass-weighted mean μ_{log n} and dispersion σ_{log n} (and the skewness 𝒮_{log n} where quoted below).

Fig. 12. Probability density functions (PDFs) of gas density measured in three radial shells, spanning the micro, meso, and inner macro-scales. The top two rows show the evolution within the inner r ≤ 0.1 kpc region (micro-scale), the middle two rows correspond to 0.1 < r ≤ 1 kpc (meso-scale), and the bottom two rows to 1 < r ≤ 10 kpc (inner macro-scale). For each radial shell, the upper row displays the cca_high (strong turbulence) run, while the lower row shows the cca_low (weak turbulence) case. Different columns represent successive snapshots, with times normalised to the rain timescale t/train. Coloured curves trace the various thermal phases: hot hard X-ray (red), hot soft X-ray (orange), warm (violet), cold (blue), and molecular (cyan).

The innermost 100 pc region (micro-scale), where most of the cooling and condensation occur, is shown in the top two panels of Figure 12. In cca_high, the initially hot and nearly homogeneous core (n ∼ 0.1 cm⁻³) exhibits a narrow, nearly lognormal density PDF, as expected for subsonic, pressure-supported turbulence in the hot phase. After t/t_rain ≈ 1 a cold molecular component with μ_log n ≃ 2.0 dex and a relatively small dispersion (σ_log n ≃ 0.2 dex) has emerged, with strongly negative skewness, indicating that most of the molecular mass sits in a well-defined high-density peak with only a shallow tail towards lower densities. Cold and warm gas populate smaller densities (μ_log n ≃ 1.4 dex and 0.3 dex, respectively) with moderate dispersions (σ_log n ∼ 0.3–0.5 dex) and near-symmetric skewness, tracing filaments and clumps in the intermediate stage of the condensation cascade. At t/t_rain ≃ 2 there is a slight broadening of the PDFs in all phases, while maintaining a similar average and distribution. At later times (t/t_rain ≃ 3) the molecular PDF broadens a lot, the mean shifts to μ_log n ≃ 4.3 dex with σ_log n ≃ 1.6 dex, and skewness becomes mildly positive, which could signal the recurrent formation and disruption of dense clumps. The hot X-ray phases, in turn, become more dilute (soft X-ray μ_log n drops from ≃ −0.4 dex to ≃ −1.2 dex) and develop positive skewness, which is consistent with cooling interfaces embedded in this volume-filling medium, which indicates the hot gas from which the condensation cascade is generated. Also the distributions of the warm and cold phases have broadened over time, spanning intermediate densities between the hot and molecular gas. Taken together, the broadening and skewness evolution likely reflect intermittent compressions and mixing layers at hot–warm–cold interfaces; in the companion paper (B26b) we will corroborate this interpretation using scale-dependent velocity statistics and kinematic CCA diagnostics (such as k-plots; Gaspari et al. 2018).

Also in the cca_low simulation the gas becomes strongly multiphase. By t/t_rain ∼ 1, a molecular component appears with μ_log n ≃ 1.2 dex and moderate dispersion, while cold and warm phases occupy intermediate densities (μ_log n ≃ 1.6 dex and 0.4 dex). The hot gas remains relatively dense and positively skewed. By t/t_rain ≃ 3, the coldest gas reaches high densities (μ_log n ≃ 5.3 dex for the molecular phase), but its PDF stays narrower than in cca_high and the hot components retain higher mean densities and very positive skewness, pointing to rare cooling sites embedded in a largely smooth hot core. In cca_high, the hot gas PDFs on average extend to lower densities (≲ 10⁻² cm⁻³), reflecting the stronger compression–rarefaction cycle driven by turbulence in this regime.

At meso-scales (middle panels, 0.1 < r ≤ 1 kpc), both simulations still exhibit multiphase structures, but with different intensities. In cca_high, the hot gas becomes progressively more dilute and structured as condensation proceeds. From t/t_rain ≃ 1 to 3, the soft X-ray phase shifts from μ_{log n} ≃ −0.6 dex to ≃ − 1.1 dex, with dispersions increasing from σ_{log n} ∼ 0.1 dex to ∼0.2 dex and skewness evolving from mildly positive to mildly negative values. The hard X-ray component follows a similar trend, maintaining μ_{log n} ≃ −0.8 to −1.2 dex with small dispersions (σ_{log n} ∼ 0.04–0.11 dex) and fluctuating skewness, indicative of intermittent turbulent compressions and rarefactions. This evolution reflects the continuous depletion of the densest hot gas as it cools and feeds the colder phases. The warm and cold components populate intermediate densities throughout this radial range, with typical means μ_{log n} ≃ −0.4–0.1 dex and μ_{log n} ≃ 0.7–1.3 dex respectively, and moderately broad PDFs (σ_{log n} ∼ 0.3–0.5 dex). Their skewness is generally mildly negative or near zero. The molecular phase is the most intermittent: its PDF alternates between relatively narrow and very broad shapes, with μ_{log n} ranging from ∼0.6 to ∼3 dex and dispersions up to σ_{log n} ∼ 0.6 dex, tracing episodic formation and destruction of dense clumps driven by CCA.

In cca_low, the meso-scale atmosphere remains significantly more stable. Both the soft and hard X-ray phases retain nearly constant mean densities over time, with μ_{log n} ≃ −0.63 to −0.70 dex for the soft X-ray gas and μ_{log n} ≃ −0.83 to −0.85 dex for the hard X-ray component. Dispersions remain small (σ_{log n} ≲ 0.1 dex), and skewness stays systematically positive, particularly for the soft X-ray phase (𝒮_{log n} ∼ 1–1.6), indicating localised, non-runaway cooling perturbations embedded in an otherwise smooth hot medium. The hard X-ray PDFs remain close to lognormal at all times. Cold and warm phases are present but evolve weakly: their mean densities and dispersions change little with time, and their PDFs remain relatively narrow compared to cca_high. The molecular component appears later than in the strong-turbulence case, becoming appreciable only after t/t_rain ≈ 2, with μ_{log n} ∼ 2 dex and modest dispersions. Overall, between 0.1 and 1 kpc, cca_high develops increasingly dilute, negatively skewed hot gas and a broad, intermittent cold component, whereas cca_low maintains denser, positively skewed hot phases and a more weakly evolving multiphase structure, reflecting the reduced efficiency of turbulence-driven condensation at meso-scales.

At the inner macro-scale (bottom panels, 1 < r ≤ 10 kpc) the thermodynamic evolution of the two simulations is very different. In cca_high, the hot soft and hard X phases have an evolution similar to the other scales. Cold phases are already present from t/t_rain ≃ 1: a diffuse warm component with μ_log n ∼ −0.3 dex to − 0.1 dex and an optical component with μ_log n ∼ 0.2–0.9 dex (dispersions σ_log n ∼ 0.3–0.5 dex) trace low-level condensation and mixing at meso–macro transition scales. A molecular component appears intermittently with μ_log n ∼ 0.6–1.6 dex and broad PDFs, signaling sporadic transport of dense clumps out to a few kiloparsecs.

In cca_low, the same radial range remains much closer to a single-phase hot atmosphere for most of the evolution. From t/t_rain = 0 to 2, only the hot phases are present. Only by t/t_rain ≃ 3 a weak multiphase component emerge: faint molecular, cold, and warm phases appear with μ_log n ∼ 1.2 dex, 1.4 dex, and − 0.1 dex, respectively, and relatively narrow dispersions, indicating that the condensed gas remains a minor, centrally sourced component rather than an extended filament network.

Overall, between 1 and 10 kpc, cca_high sustains intermittent cooling and partial multiphase mixing, with cold and warm gas present – albeit in small amounts – throughout the inner halo. In cca_low, the hot phase stays denser and more stable, and multiphase structure at these radii only appears late and remains weak, confirming that most condensation is confined to the central kiloparsec.

In the outer macro-scale halo (panel not shown, 10 < r ≤ 50 kpc), both simulations are trivially dominated by diffuse hot plasma. In cca_high, only at late times (t/t_rain ≃ 2 − 3) faint warm and cold components appear, with ⟨log n⟩ ≲ 0.

In cca_low, the outer region remains almost perfectly single-phase. At all epochs, only the hard X-ray gas is present, with ⟨log n⟩ ≃ −1.54 to −1.49 and small dispersions. No warm, cold, or molecular phases develop in this region, indicating that condensation and turbulent mixing are effectively confined to smaller radii.

The density PDFs provide a compact, quantitative summary of the multiphase condensation cascade across scales. In the strong-turbulence case (cca_high), the hot phases progressively develop broader, increasingly asymmetric distributions, while the cold and molecular components become highly intermittent, with extended high-density tails that appear early and persist from the micro- and meso-scale out to the inner macro region. This behaviour is consistent with a turbulent condensation cascade in which mixing layers and repeated compressions continuously seed nonlinear overdensities, leading to recurrent formation and disruption of dense clumps. In contrast, the weak-turbulence case (cca_low) shows systematically narrower PDFs outside the central region, with hot phases remaining closer to a quasi-stationary atmosphere and cold gas largely confined to the inner scales. The resulting multiphase medium is therefore more centrally concentrated and less stochastic at meso–macro radii, indicating a reduced efficiency of turbulence-driven fragmentation and transport. Overall, the PDF evolution links the morphological picture (filaments and clumps) and the phase-diagram view into a single scale-resolved statistical framework: stronger stirring broadens the thermodynamic distributions and sustains multiphase structure to larger radii, while weaker stirring favours a more stable hot halo with localised condensation.

4. An evolutionary link between the two turbulence regimes

Although cca_high and cca_low are designed as two fixed-turbulence experiments, they do not represent intrinsically different classes of group atmospheres. Instead, they can be interpreted as two limiting manifestations of the same underlying CCA process along a feedback–feeding cycle, in which the level of stirring and the efficacy of radiative condensation vary in time. In this context, the ‘BH weather’ terminology provides a compact way to connect the phenomenology seen in the simulations to physical regimes: (i) a sunny state, in which turbulence and heating dominate over cooling so that the atmosphere remains largely hot and condensation is suppressed or strongly reduced, (ii) a stormy CCA state, in which multiphase condensation is active over an extended region and produces a filament-rich network at meso to macro radii, and (iii) a rainy CCA state, in which condensation is still present but is more centrally concentrated, leading to a compact, clumpy cold core at micro to inner-meso radii.

Our results naturally map onto this picture. The two fiducial runs bracket the outcome of CCA under different stirring levels, while preserving the same stratified halo and cooling physics. The strongly stirred case (cca_high) delays the onset of raining (with t_rain = 30 Myr versus 7 Myr in cca_low) and sustains a broader, more radially extended multiphase medium: cold and warm structures persist to larger radii, the density PDFs broaden markedly beyond the central region, and the morphology is dominated by a tangled filamentary network that reaches kiloparsec scales. Conversely, the weak-turbulence case (cca_low) exhibits faster condensation but a more localised precipitation pattern: the hot halo outside the central region remains closer to a quasi-stationary configuration with comparatively narrower PDFs, while cold gas is largely confined to the inner tens of parsecs and organises into a compact, clumpy multiphase core. These differences show that the same halo can host qualitatively distinct ‘weather’ outcomes depending on how turbulence redistributes, mixes, and compresses the gas during the condensation cascade.

A simple evolutionary interpretation, schematically summarised in Figure 13, is as follows. In a low-stirring phase (rainy-like), condensation is concentrated towards the nucleus and the cold phase is spatially confined, which can in principle promote stronger central inflow and trigger an AGN response (winds or jets) in a full feedback model. The injected energy and momentum would then heat the atmosphere and raise the turbulent velocity dispersion, potentially shifting the system towards a hotter, turbulence-dominated state in which multiphase condensation is temporarily reduced (sunny-like). As the driving weakens and turbulence decays, radiative cooling progressively regains importance and condensation can restart over a broader radial range, producing a filament-rich, extended precipitation configuration (stormy-like) reminiscent of cca_high. With further decay of stirring and continued cooling, the condensation region can contract again and cold structures can collapse back towards the nucleus, returning the system towards a more centrally concentrated, rainy-like configuration closer to cca_low. In this narrative, stormy and rainy represent two CCA flavours regulated primarily by how turbulence distributes multiphase gas in radius, rather than two distinct accretion modes.

Fig. 13. Cartoon scheme illustrating a possible evolutionary sequence driven by the interplay between turbulence and gas condensation. Rainy phase (cca_low): in a low-turbulence state, cold clumps and filaments are concentrated near the nucleus, leading to coherent accretion and a centrally confined cold phase. Sunny phase (turb_low/turb_high): increased turbulent stirring redistributes the cold gas and temporarily suppresses further condensation, producing a more diffuse, extended atmosphere. Stormy phase (cca_high): in a strongly stirred cooling atmosphere, multiphase filaments reform over a broader radial region, before progressively collapsing towards the nucleus as turbulence weakens, ultimately returning the system to the initial, centrally concentrated state.

The characteristic timescales in our simulations support the plausibility of such a cycle. For example, the onset of raining differs by a factor ∼4 between the two endpoints (t_rain = 30 Myr in cca_high versus 7 Myr in cca_low): stronger stirring increases mixing and non-thermal support, leading to a longer raining time and a more extended multiphase structure once condensation develops, whereas weaker stirring produces earlier raining and a compact multiphase core. While our setups impose constant turbulent forcing and therefore cannot demonstrate temporal transitions within a single run, they show that modest changes in turbulence amplitude are sufficient to move the same group-scale halo between an extended, filament-rich CCA configuration and a centrally concentrated, compact multiphase state. The companion paper (B26b) will further connect the weather classification to kinematic CCA diagnostics and time-domain variability, providing an independent corroboration of the physical picture proposed here.

Differences in weather states may also arise from differences in halo mass. Lower-mass groups, with shallower gravitational potentials, may be more susceptible to AGN-driven uplift and turbulent stirring, allowing cold gas to be displaced more efficiently and favouring extended stormy configurations. Conversely, more massive groups may confine the condensed phase more effectively, promoting faster recycling towards the nucleus and producing compact, centrally concentrated rainy-like reservoirs. Thus, the distinction between stormy and rainy weather may reflect both temporal variability within individual systems and a halo-mass dependence in the likelihood of different CCA states.

The companion jet-regulated BLACKHOLEWEATHER simulations including self-consistent AGN feedback (cf. C26a,b) suggest that, in those simulations, a fully sunny state is rarely achieved. Even when AGN-driven turbulence and heating suppress condensation in the central region, multiphase gas can continue to survive and condense at larger radii, where jet-driven uplift, mixing, and compression around the outflow cones can locally trigger thermal instability and promote cold gas formation. This naturally gives rise to an intermediate ‘cloudy’ regime, in which the inner atmosphere is hotter and more turbulence-dominated while cold clouds persist at the meso-scale. Complementary BLACKHOLEWEATHER simulations including SMBH spin evolution and jet precession (cf. P26a,b) instead show transitions from rainy or stormy configurations towards a sunny phase at the meso-scale, with sub-Bondi accretion rates (10⁻⁴ − 10⁻⁵ M_⊙ yr⁻¹). A possible interpretation is that jet precession distributes the feedback energy more isotropically across the central atmosphere, increasing the efficiency of volumetric heating and suppressing condensation more globally. This suggests that the long-term weather cycle may depend sensitively not only on the feedback power, but also on the geometry and coupling of AGN energy injection.

5. Comparison with previous works

Our setup is intentionally based on the stratified-halo CCA framework introduced in Gaspari et al. (2013, 2017) which showed that in a turbulently stirred atmosphere radiative cooling can trigger a nonlinear top-down multiphase condensation cascade: warm filaments and colder clouds form out of the hot phase, inherit its turbulent kinematics, and rain towards the centre where chaotic cloud collisions promote angular-momentum cancellation and strongly time-variable SMBH feeding. The phenomenology we recover is largely consistent with these key CCA expectations: extended multiphase structures arise once cooling becomes competitive with turbulent stirring with correlated thermo-kinematics, the accretion rate is recurrently boosted by orders of magnitude relative to the classical hot-mode baseline, and the fueling remains highly bursty. At the same time, our results advance the CCA framework in a regime that is particularly relevant to link the meso-micro scales under different BH weather conditions. By explicitly contrasting two subsonic turbulence levels, we show that changes in the turbulent weather (stormy vs. rainy) are sufficient to move the same group-scale halo between an extended, filament-rich CCA precipitation state and a more centrally concentrated CCA configuration with a compact, clumpy inner structure. Advancing here the dynamic range in one calculation further allows us to follow how large-scale filaments progressively fragment into smaller clumps and streams across the three main phases with improved spatial and temporal resolution. This provides an important controlled reference baseline for the forthcoming BLACKHOLEWEATHER extensions that include jets, spin, and dust, and for more direct confrontation with multi-wavelength constraints.

Other numerical works studied SMBH accretion with high-resolution simulations, albeit with significantly different conditions. Cho et al. (2024) developed a multizone GR-MHD framework that self-consistently bridges the Bondi scale down to the event horizon. While their simulations achieve a global steady-state solution over several decades in radius, the absence of radiative cooling restricts the accretion flow to a purely hot mode. As a result, despite realistic outer boundary conditions and strong magnetic fields, the mass accretion rate remains strongly suppressed (∼1% of the Bondi rate), with no contribution from cold or multiphase accretion. Guo et al. (2023) simulated an M87-like elliptical galaxy with radiative cooling and distributed heating down to the SMBH micro-scales (with an aggressive first-order flux correction), finding that despite cold discs and chaotic phases seeded by small initial perturbations, the innermost accretion flow remains hot-mode dominated with suppressed accretion rates. Guo et al. (2024) extended this to MHD, showing that magnetic fields create filamentary cold inflows, enhance angular momentum transport, and boost accretion by an order of magnitude while launching powerful polar outflows. At variance, we use 3D hydrodynamic simulations of a hot intragroup halo with driven subsonic turbulence that models realistic astrophysical turbulence injection and weather (e.g. Wittor & Gaspari 2020), as found by multi-wavelength observations. Rather than seeding the cold phase through small decaying initial perturbations, our turbulence continuously drives multiphase condensation, allowing us to track multiple realistic raining episodes and isolate how turbulence shapes the multiphase structure.

6. Summary and conclusions

In this work, we analysed the process of gas accretion onto SMBHs in group-scale halos, focusing on how turbulence regulates the condensation, structure, and feeding of multiphase gas. We performed a set of high-resolution (Δx_min ≃ 0.1 pc) 3D hydrodynamical simulations with ATHENAPK, including radiative cooling and driven, subsonic turbulence in a hot, stratified atmosphere representative of the IGrM. Our two fiducial runs explore contrasting turbulence regimes, a strongly stirred halo (cca_high) and a weakly stirred one (cca_low), and are evolved through recurrent raining episodes. Thanks to the SMR, we follow the condensation cascade from tens of kpc across the crucial meso-scale (parsecs to kiloparsecs) down to the inner 0.1 kpc, resolving well within the Bondi radius and explicitly tracking how cold clouds and filaments feed the SMBH. Our main results can be summarised as follows.

• (i)

  In both turbulence regimes, radiative cooling in a stratified, turbulently stirred intragroup medium triggers nonlinear multiphase condensation and rain, giving rise to CCA. Turbulence primarily regulates the radial distribution and morphology of the condensed phase: strong stirring (cca_high) delays the onset of raining and sustains extended, filamentary structures reaching kpc scales (stormy weather), whereas weaker stirring (cca_low) yields a more centrally concentrated, clumpy multiphase core, with most cold gas confined within the inner tens of parsecs (rainy weather).

• (ii)

  The meso-scale (parsecs to kiloparsecs) is an active, controlling regime rather than a passive bridge. In this range, filaments repeatedly form out of the hot phase, fragment into clumps and streams, and interact through collisions and shearing layers before being funneled towards the inner region. This sustained, multiscale condensation cascade couples halo precipitation to micro-scale inflow, in part mediated by a clumpy, rotating nuclear structure reminiscent of an AGN torus.

• (iii)

  We directly follow cold clumps and filaments down to sub-pc radii and measure the SMBH accretion rate at the sink. In both cca_high and cca_low, the accretion rate is recurrently boosted by ∼10 − 100 relative to the classical hot-mode (Bondi-like) baseline, with stochastic excursions of up to ∼2 dex, reflecting the intrinsically chaotic nature of CCA.

• (iv)

  Although the strongly stirred case develops a much larger cold and molecular reservoir than the weakly stirred case, the accretion rates remain comparable. This demonstrates that SMBH feeding is not set by the total condensed mass, but by how efficiently the condensed structures couple to the central inflow once they form.

• (v)

  Cooling in a turbulent halo produces a strongly multiphase medium spanning ∼10 orders of magnitude in density and ∼8 in temperature. In cca_high, stronger stirring broadens thermodynamic distributions and sustains warm and cold gas to larger radii, yielding markedly broader density PDFs beyond the nucleus. In cca_low, density contrasts outside the centre remain smaller, the outer halo stays closer to a classical hot atmosphere, and the density PDFs beyond ∼100 pc are narrower and more time-stationary. Across phases, incorporating turbulent support reduces apparent thermal pressure imbalances, bringing the multiphase medium closer to local pressure balance in terms of P_tot = P_th + P_nt.

Taken together, these results show that the explored turbulence regimes are best interpreted as different ‘BH weather’ realisations within the same CCA process. Stronger stirring delays condensation and favours an extended, filament-rich precipitation pattern (a longer, hot turbulence-dominated sunny stage in a full cycle, followed by stormy CCA), whereas weaker stirring promotes earlier condensation and a compact, centrally concentrated multiphase core (rainy CCA). In this first paper we therefore focus on the thermodynamic and morphological manifestation of multiphase condensation across scales, quantified through phase structure, PDFs, and mass budgets. The companion paper (Barbani et al. 2026) builds directly on the same simulations to quantify the condensation criteria (via CCA diagnostics such as 𝒞-ratios and k-plots) and to connect the weather states to time-domain inflow variability and its statistical signatures.

Overall, the simulations presented here provide a new step towards a unified, multiscale theory of AGN feeding in group atmospheres. By combining stratified intragroup ICs, driven subsonic turbulence, radiative cooling, and GPU-enabled dynamic range, we follow in a single calculation the full CCA cascade from halo rain at tens of kpc down to ∼0.1 pc. In doing so, we corroborate the core CCA framework established by Gaspari et al. (2013, 2017, 2020) and extend it by quantifying how distinct ‘BH weather’ regimes imprint on multiphase morphology and thermodynamics across the meso-scale. This controlled reference baseline sets the stage for forthcoming BLACKHOLEWEATHER studies that will incorporate additional physics and confront the predicted diagnostics with multi-wavelength observations across a broad range of environments and scales.

Acknowledgments

We thank the anonymous referee for constructive comments and suggestions. The BHW authors acknowledge funding support from the European Research Council (ERC) under the European Union’s Horizon Europe research and innovation programme (Consolidator Grant BlackHoleWeather, No. 101086804; PI: Gaspari). Views and opinions expressed are, however, those of the author(s) only and do not necessarily reflect those of the European Union or the European Research Council Executive Agency; neither the European Union nor the granting authority can be held responsible for them. We acknowledge ISCRA for awarding this project access to the LEONARDO supercomputer, owned by the EuroHPC Joint Undertaking, hosted by CINECA (Italy). The numerical work was in part supported by the NASA High-End Computing (HEC) Program through the NASA Advanced Supercomputing (NAS) Division at Ames Research Center. VO acknowledges support from the DICYT ESO-Chile Comite Mixto PS 1757, and Fondecyt Regular 1251702. FMM acknowledges support from the Next Generation EU funds within the PNRR, Mission 4 - Education and Research, Component 2 – From Research to Business (M4C2), Investment Line 3.1 – Strengthening and creation of Research Infrastructures (project IR0000034 “STILES”). MF acknowledges funding by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy – EXC 2121 “Quantum Universe” – 390833306. PT acknowledges support from NASA NNH22ZDA001N Astrophysics Data and Analysis Program under award 24-ADAP24- 0011. RS acknowledges funding from the CAS-ANID grant No. 220016. The authors thank Dong-Woo Kim for providing the observational data used in Figure 2, Sara Sandoni for providing the illustration used in Figure 13, and Philipp Grete for support with the ATHENAPK code. We thank the organisers and participants of the following conferences for the stimulating discussions that helped improve this work: ‘BlackHoleWeather I’ (Sexten, ITA), ‘Modelling of Multiphase Astrophysical Media’ (Ringberg, GER), ‘Multi-phase, Multi-temperature, and Complex’ (Olbia, ITA).

References

All Tables

All Figures

Fig. 1. Simulation grid and SMBH sink scheme. The domain is decomposed into blocks (shown as a green square), which are recursively refined towards the centre of the domain (each block contains 323 cells). The decreasing block size increases the effective resolution in the central region, where the black hole sink is located. The black hole sink (shown in black) removes mass within a radius of four finest-level cells at each time step (an example of cell is shown as a red square).

In the text

Fig. 2. Thermodynamic properties of the simulated galaxy group at t = 0 Myr (black lines). Left: radial gas number density profile compared with observed electron densities from the ACCEPT database (Cavagnolo et al. 2009, pink circles). Right: normalised temperature profile (T/Tmax) as a function of r/Rvir, compared with early-type galaxies from the Chandra Galaxy Atlas (Kim et al. 2020, pink circles). The ICs are consistent with observed systems, and the temperature profile follows a nearly universal shape with a broad peak around r ≃ 0.04 Rvir. The vertical dashed line marks the maximum radius contained in the computational domain (rmax = 43.3 kpc, corresponding to half of the box diagonal).

In the text

Fig. 3. Radiative cooling function at solar metallicity used in our simulations. The curve shows Λ(T) as a function of temperature, with tabulated values from Schure et al. (2009) for T ≥ 104.2 K and the analytic fit given by Eq. (13) for T < 104.2 K. Cooling function includes metal-line, recombination, bremsstrahlung and low temperature cooling; it is tabulated on a uniform grid in log T and interpolated during the simulation to compute radiative losses.

In the text

Fig. 4. Gas density slices of the central 50 kpc region in the cca_high simulation. Rows show the time evolution from t/train = 1.5 (top) to t/train = 3 (bottom), while within each row the panels from left to right present a progressive zoom-in from the halo scale (macro) to the filamentary condensation region (meso) and finally to the innermost clumps (micro). The sequence reveals the characteristic morphology of CCA: a complex network of cold, dense filaments condensing out of the turbulent hot halo and converging towards the galactic centre. Panels at different spatial scales adopt independent logarithmic colour ranges to optimise contrast and highlight the gas substructures; colour scales are thus not directly comparable across panels. The top right corner of each panel in the last column shows a zoom-in of the innermost rotating structure.

In the text

Fig. 5. Gas density slices of the central region in the cca_low simulation. Rows show the time evolution from t/train = 1.5 (top) to t/train = 3 (bottom), while within each row the panels from left to right present a zoom-in from the inner kiloparsec to 15 pc around the SMBH. At both epochs, the system hosts a dense, clumpy cold core on parsec scales, continuously fed by smoother filamentary inflows from larger radii. The top right corner of each panel in the last column shows a zoom-in of the innermost rotating structure. Cooling dominates over turbulent mixing, allowing cold gas to accumulate and recycle within a long-lived multiphase region where clouds and filaments move along intersecting, colliding, and shearing trajectories. The morphology remains irregular and fragmented, illustrating the CCA mode also in a weakly turbulent galactic atmosphere.

In the text

Fig. 6. Slices of the simulated galaxy group core for the cca_high (top row) and cca_low (bottom row) simulations at t/train = 3, showing temperature (left), pressure (middle), and velocity magnitude (right). The maps highlight the multiphase structure of the IGrM and the role of turbulence in driving thermodynamic fluctuations in the hot atmosphere, seeding local thermal instability and the condensation of cold gas. Each panel in the top row shows an 8 kpc × 8 kpc region, while each panel in the bottom row shows the central 0.2 kpc × 0.2 kpc region.

In the text

Fig. 7. SMBH accretion rate Ṁ normalised to the Bondi accretion rate ṀB as a function of normalised time t/train for the cca_high (blue line, stormy weather) and the cca_low (cyan line, rainy weather) simulations, compared with the turbulence simulations (sunny weather) turb_high (red line) and turb_low (yellow line), with the radiative cooling simulation cool (purple line) and with the idealised adiabatic simulation bondi (black line). The CCA phenomenon produces a highly variable accretion rate due to its chaotic nature.

In the text

Fig. 8. Mass-weighted radial profiles of gas density, temperature, thermal pressure, and velocity magnitude for simulations cca_high (top row) and cca_low (bottom row). The gas is divided into five temperature phases: hot, hard X-ray (red) and soft X-ray (orange), warm (violet), cold (blue) and molecular (cyan, see Table 4). Increasing colour intensity indicates later stages of evolution. Profiles are computed in spherical shells around the group centre out to 25 kpc, showing how the thermodynamic and kinematic structure of the multiphase IGrM evolves in the two CCA regimes. Pink circles show observational density, projected temperature, and pressure profiles from the ACCEPT database (Cavagnolo et al. 2009), where applicable.

In the text

	Fig. 9. Mass-weighted density-temperature phase diagrams for the cool (blue), turb_high (red) and turb_low (orange) simulations. Dashed oblique lines indicate adiabatic (T ∝ nγ − 1 = n2/3) and isobaric (T ∝ n−1) relations.
In the text

Fig. 10. Mass-weighted density–temperature phase diagrams for the cca_high (rows 1 and 3) and cca_low (rows 2 and 4) simulations. The first two rows show the time evolution at t/train = 0.1,  1,  2,  3, while the last two rows show the corresponding distributions at t/train = 3 in four radial ranges: 0 < r ≤ 0.1 kpc (micro-scale), 0.1 < r ≤ 1 kpc (meso-scale), 1 < r ≤ 10 kpc (inner macro-scale), and r > 10 kpc (outer macro-scale). Colours denote the total gas mass in each (n, T) bin, from dark brown (low mass – 10−5 M⊙) to yellow (high mass – 105 M⊙). Phases (a-f) are classified by temperature (hot: T ≳ 106 K; warm: 104 < T < 106 K; cold: T < 104 K) and density (diffuse vs. dense, with thresholds at n = 1 cm−3 for hot/warm and n = 100 cm−3 for cold). Dashed oblique lines indicate adiabatic (T ∝ n2/3) and isobaric (T ∝ n−1) relations in each panel.

In the text

	Fig. 11. Gas mass for cca_high (solid) and cca_low (dashed) simulations divided in molecular (cyan), cold (blue), warm (violet), soft (orange), and hard (red) X-ray gas and in different scales: macro, meso, inner macro and outer macro.
In the text

Fig. 12. Probability density functions (PDFs) of gas density measured in three radial shells, spanning the micro, meso, and inner macro-scales. The top two rows show the evolution within the inner r ≤ 0.1 kpc region (micro-scale), the middle two rows correspond to 0.1 < r ≤ 1 kpc (meso-scale), and the bottom two rows to 1 < r ≤ 10 kpc (inner macro-scale). For each radial shell, the upper row displays the cca_high (strong turbulence) run, while the lower row shows the cca_low (weak turbulence) case. Different columns represent successive snapshots, with times normalised to the rain timescale t/train. Coloured curves trace the various thermal phases: hot hard X-ray (red), hot soft X-ray (orange), warm (violet), cold (blue), and molecular (cyan).

In the text

Fig. 13. Cartoon scheme illustrating a possible evolutionary sequence driven by the interplay between turbulence and gas condensation. Rainy phase (cca_low): in a low-turbulence state, cold clumps and filaments are concentrated near the nucleus, leading to coherent accretion and a centrally confined cold phase. Sunny phase (turb_low/turb_high): increased turbulent stirring redistributes the cold gas and temporarily suppresses further condensation, producing a more diffuse, extended atmosphere. Stormy phase (cca_high): in a strongly stirred cooling atmosphere, multiphase filaments reform over a broader radial region, before progressively collapsing towards the nucleus as turbulence weakens, ultimately returning the system to the initial, centrally concentrated state.

In the text