Constraining hot gas turbulence from the warm and cold phase

# Shaken Snow Globes: Kinematic Tracers of the Multiphase Condensation Cascade in Massive Galaxies, Groups, and Clusters

M. Gaspari, M. McDonald, S. L. Hamer, F. Brighenti, P. Temi, M. Gendron-Marsolais, J. Hlavacek-Larrondo, A. C. Edge, N. Werner, P. Tozzi, M. Sun, J. M. Stone, G. R. Tremblay, M. T. Hogan, D. Eckert, S. Ettori, H. Yu, V. Biffi, S. Planelles Department of Astrophysical Sciences, Princeton University, 4 Ivy Lane, Princeton, NJ 08544-1001, USA
Kavli Institute for Astrophysics and Space Research, MIT, Cambridge, MA 02139, USA
CRAL, Lyon Observatory, CNRS, Université Lyon 1, 9 Avenue Charles André, F-69561 Saint Genis-Laval, France
Astronomy Department, University of Bologna, Via Piero Gobetti, 93/3, 40129 Bologna, Italy
Astrophysics Branch, NASA Ames Research Center, Moffett Field, CA 94035, USA
Department of Physics, University of Montreal, Montréal, QC H3C 3J7, Canada
Department of Physics, Durham University, Durham, DHL 3LE, United Kingdom
MTA-Eötvös University Lendület Hot Universe Research Group, Pázmány Péter sétány 1/A, Budapest, 1117, Hungary
Dep. of Theoretical Physics and Astrophysics, Faculty of Science, Masaryk University, Kotlářská 2, Brno, 611 37, Czech Republic
School of Science, Hiroshima University, 1-3-1 Kagamiyama, Higashi-Hiroshima 739-8526, Japan
INAF, Astronomy Observatory of Florence, Largo Enrico Fermi 5, 50125, Firenze, Italy
Physics Department, University of Alabama in Huntsville, Huntsville, AL 35899, USA
Harvard-Smithsonian Center for Astrophysics, 60 Garden Street, Cambridge, MA 02138, USA
Department of Physics and Astronomy, University of Waterloo, Waterloo, ON, N2L 3G1, Canada
Max Planck Institute for Extraterrestrial Physics, Giessenbachstr., 85741, Garching, Germany
INAF, Astronomy Observatory of Bologna, Via Piero Gobetti, 93/3, 40129 Bologna, Italy
INFN, Sezione di Bologna, viale Berti Pichat 6/2, 40127 Bologna, Italy
Department of Astronomy, Beijing Normal University, Beijing, 100875, China
Department of Physics and Astronomy, University of Trieste, via Tiepolo 11, 34131 Trieste, Italy
INAF, Astronomy Observatory of Trieste – OATs, via Tiepolo 11, 34131 Trieste, Italy
Department of Astronomy and Astrophysics, University of Valencia, C/Dr. Moliner 50, 46100 Valencia, Spain
###### Abstract

We propose a novel method to constrain turbulence and bulk motions in massive galaxies, galaxy groups and clusters, exploring both simulations and observations. As emerged in the recent picture of the top-down multiphase condensation, the hot gaseous halos are tightly linked to all other phases in terms of cospatiality and thermodynamics. While hot halos ( K) are perturbed by subsonic turbulence, warm ( K) ionized and neutral filaments condense out of the turbulent eddies. The peaks condense into cold molecular clouds ( K) raining in the core via chaotic cold accretion (CCA). We show all phases are tightly linked in terms of the ensemble (wide-aperture) velocity dispersion along the line of sight. The correlation arises in complementary long-term AGN feedback simulations and high-resolution CCA runs, and is corroborated by the combined Hitomi and new Integral Field Unit measurements in Perseus cluster. The ensemble multiphase gas distributions (from UV to radio band) are characterized by substantial spectral line broadening ( - ) with mild line shift. On the other hand, pencil-beam detections (as HI absorption against the AGN backlight) sample the small-scale clouds displaying smaller broadening and significant line shift up to several 100  (for those falling toward the AGN), with increased scatter due to the turbulence intermittency. We present new ensemble of the warm H+[NII] gas in 72 observed cluster/group cores: the constraints are consistent with the simulations and can be used as robust proxies for the turbulent velocities, in particular for the challenging hot plasma (otherwise requiring extremely long X-ray exposures). Finally, we show the physically motivated criterion best traces the condensation extent region and presence of multiphase gas in observed clusters and groups. The ensemble method can be applied to many available spectroscopic datasets and can substantially advance our understanding of multiphase halos in light of the next-generation multiwavelength missions.

multiphase ICM, IGrM, CGM – AGN feedback – 3D hydrodynamic simulations – spectroscopic observations – turbulence – X-rays, UV, optical, radio: galaxies, groups, clusters
* * affiliationtext: E-mail: mgaspari@astro.princeton.edu $\dagger$ $\dagger$ affiliationtext: Einstein and Spitzer Fellow

## 1. Introduction

Despite our everyday solid-state experience, baryons populate the universe mostly in a diffuse gaseous form. A new picture has recently emerged – from both the theoretical and observational side – describing the gaseous atmospheres of galaxies, groups, and clusters of galaxies. While initially modeled as hydrostatic monophase systems, the gaseous halos filling the potential well of cosmic systems are complex atmospheres akin to Earth weather, following a top-down multiphase condensation cascade (e.g., Gaspari:2017_cca; Gaspari:2017_uni). After falling at large redshift into the potential wells of dark matter halos, baryons heat up, forming hot plasma halos (intracluster, intragroup, and circumgalactic medium – ICM, IGrM, CGM; McNamara:2012; Sun:2012 for reviews). Such hot halos are the progenitors for other major condensed structures, including warm filaments, cold molecular clouds, and stellar/planetary systems.

During their evolution, the diffuse halos experience cyclical states, akin to the rapid alternation on Earth of sunny, cloudy, and rainy weather. From the thermal point of view, cosmic atmospheres span temperatures from several keV ( K) of hot plasma halos to  K of warm ionized and neutral filaments to tens K of cold molecular clouds (as beautifully detected by ALMA), with particle number density on average anticorrelated with temperature ( - ). At the same time, from the dynamical point of view, cosmic atmospheres experience a continuous competition between chaotic turbulent motions and coherent rotational flows (turbulent Taylor number or , respectively). Hotter, thermal pressure-supported halos often reside in the former chaotic regime due to the multiple drivers acting through cosmic time in a partially uncorrelated way: at larger radial distances (Mpc-scale) mergers and galaxy motions drive subsonic turbulence in the volume-filling phase (e.g., Vazza:2011; Miniati:2014; Khatri:2016), while in the core ( kpc – where the entropy profile slope changes) active galactic nucleus (AGN) feedback recurrently pumps energy via massive outflows and jets (e.g., Lau:2017; Hillel:2017); at the smallest scales, supernovae and stellar winds further preserve a minimum level of (compressive) turbulence (e.g., Kim:2013).

In the turbulent gaseous halos of clusters, groups, and galaxies (particularly massive ones), extended filaments and clouds condense out of the hot plasma in a top-down nonlinear111This nonlinear condensation process has significantly different properties from those of classic linear thermal instability (TI); the latter is mainly concerned with small overdensities overcoming buoyancy oscillations (e.g., Field:1965; Balbus:1989; Burkert:2000; Pizzolato:2005; McCourt:2012 – more in §5). condensation cascade, forming a chaotic multiphase rain. The thermal state and kinematics of the progenitor hot plasma halo drive the formation and evolution of all the condensed structures, which inherit some of the parent properties. Part of the inner condensed gas eventually accretes onto the central supermassive black hole (SMBH), igniting the feedback response and efficiently self-regulating the whole atmosphere over several Gyr (e.g., Gaspari:2011a; Gaspari:2011b; Gaspari:2012b; Gaspari:2012a; Li:2014; Barai:2016; Yang:2016; Soker:2016; Meece:2017; Voit:2017). This feeding process is known as chaotic cold accretion (CCA; Gaspari:2013_cca) and can intermittently boost the accretion rates up to the hot (Bondi) rate. If turbulence is subdominant, the halo tends instead to condense in a disk structure (due to the preservation of angular momentum), reducing feeding and feedback – this regime is more important for low-mass, spiral galaxies222The top-down rain differs from the bottom-up condensation in the disk of spiral galaxies, where the hot/warm phase is created in situ by supernovae which drive compressive, non-solenoidal turbulence (e.g., McKee:1977; Kim:2013). Nevertheless, the two complement each other, producing multiphase gas in the more extended halo and in the disk, respectively. Massive galaxies, groups, and clusters, lacking an extended disk (e.g., Werner:2014), typically reside in the top-down condensation regime.. Finally, if the entropy of the halo (or cooling time) becomes too high, the whole atmosphere may simply prevent condensation and remain hot for an extended period of time, dramatically stifling the feedback response. Overall, assessing the dynamical state of the multiphase halos is crucial to understand the past and predict the future evolution of cosmic structures.

Although the thermal properties of gaseous halos are fairly well-constrained thanks to the last-generation X-ray, optical/IR, and radio telescopes (e.g., Chandra, XMM, Hubble, Herschel, and IRAM; Combes:2007; McDonald:2010; McDonald:2011a; McNamara:2012; Canning:2013; Werner:2014; Tremblay:2015; Hamer:2016; Russell:2016; David:2017), constraining the kinematics of the hot phase has proven to be very challenging, mainly due to the limited spectral resolution at high energies. The kinematics of the gas can be directly retrieved from the spectral line width (which is tied to the turbulent velocity dispersion) or the line centroid offset (which traces bulk motions). Recently, Hitomi gave us a sneak peek into the complexity of hot halos, finding line-of-sight (LOS) velocity dispersions in the Perseus cluster core (Hitomi:2016). Turbulent motions can also be roughly estimated via relative plasma density fluctuations, which are related to the turbulent Mach number (e.g., Gaspari:2013_coma; Hofmann:2016; Eckert:2017_PS; Zhuravleva:2017), finding subsonic Mach numbers in the ICM, although substructures contamination can introduce a significant noise. The subsonic turbulence is corroborated by the linewidth upper limits in combination with resonant scattering set by XMM-RGS (Werner:2009; dePlaa:2012; Sanders:2013; Pinto:2015; Ogorzalek:2017). Such a level of turbulence is also required to substantially suppress the emission measure in the soft X-ray spectrum (Gaspari:2015_xspec).

This paper continues our systematic investigation of the multiphase condensation and CCA mechanism (e.g., Gaspari:2012a; Gaspari:2013_cca; Gaspari:2015_cca; Gaspari:2017_cca), focusing on the gas kinematics. By using state-of-the-art high-resolution hydrodynamic 3D simulations, complemented by new observations, we present a novel method to constrain the gas motions, taking advantage of the ensemble333Theoretically, meaning the global large-volume statistics of all of the condensed elements for a given phase; practically, referring to the use of (spectroscopic) observations with wide projected radial aperture ( arcmin or exceeding several kpc). kinematics of the condensed multiphase filaments and clouds – one of the most robust properties of the low-energy phases. We will show that, singularly, each structure can take a large range of values of the random and bulk velocity components. Globally, and with enough statistics, the condensed structures can be considered as quasi-linear tracers of the turbulent eddies and cascade – reminiscent of shaken snow globes. Vice versa, we can apply the same method to infer the kinematics of the cooler phase from the warm phase, or any different multiwavelength combination. As shown by the new observational constraints in §4, this can be easily and efficiently leveraged by the Integral Field Unit (IFU) spectroscopy, which is advancing at a remarkable pace (e.g., MUSE, VIMOS, SITELLE). In the other direction, small-aperture/‘pencil-beam’ (e.g., a few 100 pc or  arcsec) detections – such as HI absorption against the AGN backlight or CO emission – can shed light on the mode of accretion onto SMBHs (e.g., CCA versus hot mode accretion; David:2014; Hogan:2014; Tremblay:2016) and on the properties of the small-scale clouds (e.g., Maccagni:2017). As we live in an era of new exciting telescopes covering the radio and IR spectrum (e.g., ALMA, JWST, SKA, CARMA2) the proposed multiwavelength kinematics methods can be tested and used to advance our understanding of cosmic halos in galaxies, groups, and clusters.

Retrieving the velocity dispersion of the hot halo opens up a simple and direct way to assess the presence of multiphase gas or ensuing condensation. A major debate in the recent literature concerns which is the best (and minimal) criterion to assess the condensation state of the hot halo as a function of characteristic timescales (e.g., McCourt:2012; Sharma:2012; Gaspari:2012a; Voit:2015_nat; McNamara:2016; Hogan:2017), including  Gyr,  - 30, or , where , , and are the cooling, gravitational, and conduction timescales, respectively. We will show that the crossing locus of the cooling time and the turbulent eddy time (which is a function of predominantly the ensemble gas velocity dispersion, , and directly accessible from observations) provides a robust criterion for the physical state of the hot halo, separating multiphase and monophase systems. Although the method can be applied to a large number of available datasets, it can also augment the next-generation X-ray missions (e.g., Athena, XARM, and Lynx) by providing precise and testable observables.

This work is structured as follows. In §2, we review the high-resolution 3D hydrodynamic simulations used in this study. In §3, we dissect the resulting correlations between all of the different phases (soft X-ray to UV/optical band to radio/21 cm), in both long-term AGN feedback simulations (§3.1) and CCA feeding simulations with pc-scale resolution (§3.2). The main correlation is tested with the Hitomi and new SITELLE IFU direct measurements in the Perseus cluster. In §3.3, we discuss the current limitations of the models and future improvements. In §4, we present new observational constraints – together with available literature data – obtained via the proposed ensemble (§4.1) and pencil-beam (§4.2) methods for the warm and cold gas in massive galaxies (many of which are central brightest cluster galaxies – BCGs) and we compare them with the above numerical predictions. In §5, we discuss a key application of the ensemble measurement, presenting a new condensation criterion tied to the turbulence eddy turnover time for the presence and radial extension of the multiphase gas structures in clusters and groups. In §6, we summarize the main results of the study and provide concluding remarks.

## 2. Simulations

The core of the theoretical study (§3) is based on 3D hydrodynamic simulations (carried out with the Eulerian adaptive-mesh-refinement – AMR – code FLASH4), combining them with new and recent multiwavelength observations (§4). We use as reference two complementary simulations, one covering the large-scale and long-term evolution, and the other covering the high-resolution and full multiphase cascade from the hot plasma to the molecular phase. As we privileged high accuracy in space and time, the total computational cost is substantial, over 4 million CPU-hours. Since the simulations are unchanged compared with our previous investigations, we refer the interested reader to Gaspari:2012a (Gaspari:2012a – G12 hereafter) and Gaspari:2017_cca (Gaspari:2013_cca; Gaspari:2015_cca; Gaspari:2017_cca – G17 hereafter) for the details and nuances related to the modules and numerics adopted. Here we review the key features and relevant physics.

### 2.1. G12 simulation: self-regulated AGN feedback

The goal of the G12 suite of simulations is to study the evolution and properties of the long-term self-regulated kinetic AGN feedback affecting the X-ray plasma halo. The simulation models a typical cool-core galaxy cluster with central plasma444The plasma average particle weight is , with adiabatic index . The metal abundances are  -  for the cluster and central galaxy, respectively. entropy  keV cm (minimum ). The plasma halo is initially perturbed by random fluctuations in density and temperature with 0.3 relative amplitude to model the presence of cosmic weather (G12; Sec. 2.3). The maximum AMR resolution reaches 300 pc, so that it is possible to evolve the system for several Gyr in a large  Mpc domain. The static gravitational potential is dominated by the dark matter component with Navarro-Frenk-White profile; the cluster virial mass is with gas fraction .

Besides hydrodynamics555In all runs, we employ the third-order accurate piecewise parabolic method (PPM) to solve the Euler hydrodynamics equations. Boundary conditions have all outflow permitted and inflow prohibited., the two key competing physics are radiative cooling and AGN feedback. The plasma radiative cooling induces the gas to lose thermal energy, and thus pressure support, forming extended warm filaments via nonlinear condensation. The plasma halo emits radiation mainly via Bremsstrahlung above 1 keV and line recombination below such soft X-ray regime. The emissivity is , where is the Sutherland:1993 plasma cooling function in collisional ionization equilibrium. The plasma cooling curve incorporates calculations for H, He, C, N, O, Fe, Ne, Na, Si, Mg, Al, Ar, S, Cl, Ca, and Ni, and all stages of ionization. Due to the limited resolution in this run, condensation is halted at the warm phase regime around  K. The cooling source term is integrated with an exact solver with conservative time-step limiter (G12, Sec. 2.1).

Radiative cooling is counterbalanced by AGN feedback, in the form of massive subrelativistic outflows. The bipolar outflow mass, momentum, and energy are injected through an internal boundary nozzle in the innermost resolved region (G12, Sec. 2.2) – locus of the SMBH. The injected velocity is , which is typical of observed entrained FRI jets or ionized ultrafast outflows (e.g., Tombesi:2013). The injected kinetic power is self-regulated by the central inflow rate, , where is the mechanical efficiency able to avoid both overcooling and overheating. The triggering nuclear mass inflow ( pc) is dominated by the condensed gas phase (thus linked to the cooling rate of the hot halo) in the form of raining filaments and clouds (a.k.a. CCA – §2.2).

The self-regulated, bipolar AGN outflows propagate outwards and gently dissipate the mechanical energy via bubble mixing, turbulence, and weak shocks, thus restoring most of the previously radiated internal energy (i.e., a global quasi thermal equilibrium), with a duty cycle of the order of the central (Gaspari:2017_uni for a review). This simulation has been shown to be consistent with several spectroscopical X-ray constraints, including the suppression of the cooling flow by over 20 fold and the Gyr survival of the cool core (preserving the positive gradient; Gaspari:2013_rev).

### 2.2. G17 simulation: ultra high-resolution AGN feeding

While G12 runs model a more realistic AGN feedback injection, the G17 (and previously related) simulations aim to resolve with maximally feasible resolution ( pc) the top-down multiphase condensation cascade, from the keV plasma phase to the warm gas ( -  K) and to the cold (20 - 200 K) molecular gas. The simulation zooms in, with static mesh refinement, on the central massive galaxy (akin to NGC 5044) in the inner gaseous cool core (52 kpc domain; again, minimum plasma ), reaching a 100 Myr evolution. The static potential over such a small domain is dominated by the central galaxy stellar mass (with effective radius  kpc). The central SMBH () is modeled with the pseudo-relativistic Paczynski:1980 potential and has a characteristic Bondi radius of 85 pc (G17, Sec. 2.2).

The plasma radiative emissivity follows the same prescription as described in §2.1 () though now the cooling curve is extended down to the neutral and molecular regime (Fig. 1 in G17) following an analytic prescription analogous to reference ISM studies (e.g., Dalgarno:1972; Kim:2013; Sec. 2.5 in G17). The included processes are atomic line cooling from hydrogen Ly, CII, OI, FeII, SiII; rovibrational line cooling from H and CO; and atomic and molecular collisions with grains (Koyama:2000). The typical ionization fraction in the warm/cold phase is of the order of 1% (mimicking the influence of post-AGB stars, AGN, and cosmic rays).

As we are here interested in the feeding stage and due to the limited integration time, we model only the gentle 4 deposition stage of the mechanical AGN feedback, which prevents the cooling flow catastrophe and related monolithic collapse of the hot atmosphere (G12, Fig. 9). The injected plasma halo heating rate thus balances the average cooling rate in coarse radial shells (G17, Sec. 2.4), i.e., the hot atmosphere is globally stable but locally unstable. We further include heating of the cold/warm phase (G17, Sec. 2.5), which is mainly due to the photoelectric effect; however, since massive/early-type galaxies have minor star formation ( yr), this term is subdominant. We note that we are fully aware of the complexity of the heating and cooling microprocesses in the multiphase ICM/IGrM; these simulations are part of our ongoing numerical campaign to dissect each physics step by step (§3.3).

The additional key physics, alongside hydrodynamics, cooling, and heating, is turbulence. As probed naturally in the previous self-regulated AGN outflow runs (and independent studies; §1), the hot X-ray halo is continuously stirred by subsonic turbulence (due to the buoyant bubbles, Kelvin-Helmholtz instabilities, and shocks). Being these runs controlled experiments, the diffuse hard X-ray halo is continuously perturbed by subsonic turbulent motions () via a spectral forcing scheme based on an Ornstein-Uhlenbeck random process (G17, Sec. 2.3), which also reproduces experimental high-order structure functions (Fisher:2008). The gas is only stirred at low- Fourier modes, allowing the development of a self-consistent and natural turbulence cascade. Such a subsonic, mainly solenoidal turbulence mimics well that produced by AGN feedback during the deposition stage (e.g., Gaspari:2012b). Since we focus in this work on massive galaxies, the driven chaotic gas motions are stronger than the gas rotational velocity (i.e., ).

In this heated and turbulent atmosphere, warm filaments and cold clouds condense out of the hot plasma, some of which rain on the SMBH. In the nuclear region inelastic collisions allow the angular momentum to mix and cancel via chaotic cold accretion. As shown in our previous series of investigations (cf. Gaspari:2013_cca; Gaspari:2015_cca; Gaspari:2017_cca), CCA displays important properties that can explain diverse observed phenomena, including the rapid flickering of AGN, their obscuration properties (the broad-/narrow-line region and clumpy torus), the shallow X-ray temperature profiles, and the cospatiality of the inflowing/outflowing multiphase gas in the soft X-ray, optical, and radio bands. In this work, we focus on the multiphase CCA kinematics.

## 3. Linking the multiphase gas kinematics

While in previous works we focused on the thermodynamics and accretion process, here we analyze the LOS luminosity-weighted (LW) kinematics of the multiphase gas and possible correlations. Numerically, the mean LW LOS velocity is computed as

 ¯vlos=ΣlosvkΔLkΣlosΔLk, (1)

and the related velocity dispersion as

 σ2v,los=Σlos(vk−¯vlos)2ΔLkΣlosΔLk, (2)

where the summation is computed over a given cylindrical aperture along the full line of sight (with the SMBH as center). The luminosity for each cell with volume is , where is the radiative cooling curve (§2.2) in the temperature band of the given gas phase. In the following subsections, we mainly analyze the numerical results, while we dedicate §4 to an in-depth comparison with recent and new observations of multiphase gas in massive galaxies.

### 3.1. Long-term self-regulated kinetic AGN feedback

We start from the long-term simulation, which can follow condensation only down to the warm phase but in a long-term AGN outflow feedback evolution (G12; §2). Fig. 1 addresses a key question: what is the degree of correlation between the velocity dispersions (spectral line width) of the condensed ensemble warm gas and the volume-filling X-ray plasma in cluster cores?

During the Gyr evolution, the hot halo is continuously perturbed by the cosmic weather and the AGN outflows at large and small radii, respectively. The turbulent motions promote the nonlinear condensation of extended warm ionized filaments ( K), which are mainly observed in H+[NII] emission. Fig. 1 shows that during the top-down multiphase condensation and recurrent AGN feedback cycles (blue points), the ensemble warm phase behaves as a quasi linear tracer of the turbulent eddy evolution (see also §3.2.1). The best-fit linear relation has the following slope and normalization:

 σv,los,warm=0.97+0.01−0.02σv,los,hot+8.3+3.5−5.1 kms−1. (3)

The linear correlation emerges within a wide ensemble extraction region of size 1-2 condensation radii (for our massive cluster  45 kpc), over which the warm gas forms out of the progenitor X-ray (0.3 - 8 keV) plasma. At small scales, the single clouds show instead a larger variance driven by the local eddies and cloud-collisional kinematics (see §3.2).

The turbulence eddy turnover timescale tied to this coupling region (which is also comparable to the AGN bubble injection scale) is the effective dynamical timescale of the top-down multiphase condensation process (§5). For well-resolved objects, it is preferable – but not essential – to cut the very inner region (here 4 kpc) in the presence of jets activity666As long as the wide aperture captures the bulk of the condensed gas and related emission, the ensemble detection is not sensitive to changing the inner/outer extraction radius by a few kpc, remaining within the retrieved scatter.. The simulated velocity dispersion distribution has mean , reaching values up to 250 during the stronger AGN feedback phases. The logarithmic scatter of the entire warm/hot gas distribution is 0.13 dex. Focusing instead on the deviation from the best-fit line, the RMS is 14% (with maximum residuals up to %) which is mainly generated by the AGN duty cycle and related time hysteresis between driving perturbations and recurrent residual condensation.

In §4, we compare the simulated distribution with new ensemble warm gas constraints for 76 clusters, resulting to be consistent with the simulated range predicted here. For one cluster – Perseus – we can directly probe the correlation here, as direct LOS velocity dispersion detections for both the warm gas and X-ray plasma are available. Specifically, we combine the Fe XXV-XXVI linewidth fiducial detection, (Hitomi:2016777Recently, a few more uncertain regions have been added to the analysis, which nevertheless resulted in a similar single-spectrum value (Hitomi:2017, Sec. 3.4).) , with a new wide-field SITELLE888A wide-field imaging Fourier transform spectrometer (IFTS) with IFU capabilities in the visible (350 - 900 nm) for the Canada-France-Hawaii telescope (CFHT; Drissen:2010): http://www.cfht.hawaii.edu/Instruments/Sitelle. IFTS observation999Data taken in January 2016 with the SN3 filter (651-685 nm) for 2.14h (308 exposures of 25 s), with a spectral resolution of 1800. The data reduction and calibration were conducted by using ORBS and the analysis tools from ORCS (Martin:2015). of the ensemble H+[NII] linewidth (Gendron-Marsolais et al. in prep.; see §4.1 for the analysis method). By fitting the H+[NII] lines of the single spectrum integrated over the same wide extraction region as above, we find . Selecting the hard X-ray band as for the Hitomi iron-lines measurements ( 5 - 9 keV), we find simulated plasma velocity dispersions that are on average 20% higher than those in the entire X-ray band, since the hard X-ray, less dense gas is more easily accelerated by feedback. Nevertheless, whether or not this correction is applied (orange circle versus triangle in Fig. 1), the warm and hot gas velocity dispersions are consistent with having comparable turbulent kinematics within uncertainties, in agreement with the predicted correlation.

### 3.2. High-resolution chaotic cold accretion feeding

We move on from the long-term evolution to the detailed ultra high-resolution kinematics of the top-down multiphase condensation (0.8 pc – ensuring convergence of the main properties), which tracks all the phases down to the molecular regime (G17 and §2). While the previous simulation focuses on the realistic feedback process, the current run focuses on the detailed feeding process in a central massive galaxy for a shorter time, 100 Myr, which is still 10 times the central (kpc-scale) cooling time. In this turbulent and heated halo, extended warm filaments ( -  K) condense out of the hot keV plasma atmosphere and produce a condensation rain. The thin outer layer of the filament is ionized and strongly emitting in optical and UV. The spine of the filaments is mostly neutral gas ( K), containing most of the warm gas mass. The denser peaks further condense into molecular clouds ( K) with radii spanning several pc to 100 pc for the giant molecular associations. Total molecular masses can reach up to several , consistent with recent ALMA data (e.g., David:2014; massive clusters can show even , e.g., Vantyghem:2016; Pulido:2017). While temperature radial profiles are fairly flat for all condensed phases, density profiles have logarithmic slope -1, with inner densities up to  g cm for the molecular phase.

The CCA process is also responsible for efficiently boosting SMBH accretion rates with rapid intermittency up to two orders of magnitude (e.g., Gaspari:2016 for a brief review). Given that subsonic turbulence is a common state of hot halos (§4), CCA is also a recurrent state of observed systems (e.g., McDonald:2011a; McDonald:2012_Ha; Werner:2014; Tremblay:2016; David:2017), although overheated halos can experience a pure hot low-accretion mode, and rotation-dominated halos can be associated with a decoupled thin disk. We refer the interested reader to G17 for the detailed thermodynamic properties of each phase and in-depth discussions; here, are interested in the statistical kinematic properties related to the CCA rain, in particular confronting the ensemble versus local variance and the mean of the velocity field (i.e., the broadening and shift of the spectral lines) for all the gas phases. The limitations of the current simulations and future improvements are discussed in §3.3.

Fig. 2 shows the ensemble velocity dispersion in six major temperature bins normalized to the subsonic turbulent velocity, which is stably driven in the hot plasma ( K). The different phases correspond to key observational bands, covering the radio, optical/IR, and UV/soft X-ray regimes, as highlighted by different colors in the top panel. Turbulent LOS velocity dispersions are detected through the broadening of the observed spectral lines by measuring the line’s full width at half maximum, . The lower the temperature, the smaller the contribution of thermal broadening101010The relative turbulent and thermal Doppler broadening are respectively given by and , where is the line center frequency and is the mass of the given ion., which is  1 - 8  for molecular and warm gas, respectively. The H, CO, HI, [CII], and H+[NII] lines are all excellent probes of the gas kinematics.

The top and bottom panels of Fig. 2 show the same velocity dispersion diagnostics for the ensemble-beam (excluding the collisional nuclear region) and for a pencil-beam (aperture pc) detection, respectively. The blue bars indicate the logarithmic mean and 1-2 standard deviation111111 The average fractional difference (as absolute value) between the mean/1 and median/34.1% interval for the ensemble-beam points in log space is 7.8/8.6%, respectively. Gaussian fits are thus a good representation of the distributions and are accurate enough for the scope of the present work. of the underlying points (not shown), which are tracked every 1 Myr for 100 Myr. The ensemble measurement substantially reduces the turbulence intermittency noise and shows again a tight linear correlation throughout the phases, corroborating the result in §3.1. The RMS deviation from the hot gas turbulent velocity is 13% (brown), which is analogous to the long-term simulation deviation from the best-fit line in Fig. 1. The ratio is not unity as condensed structures do not fill the entire halo at every moment in time. This demonstrates that we can use the ensemble warm or cold gas as tracers of the kinematics of the turbulent hot gas, and vice versa we can predict the kinematics of the multiphase CCA cascade from the turbulent plasma halo. The optical/NUV phase near the stable  K regime has one of the lowest scatters and better equivalence, making nebular H+[NII] emission an excellent tool to study turbulence (§4.1). The FUV phase shows larger mean velocity dispersion, experiencing the most rapid and unstable condensation transition due to the strong line cooling, while tracing the low-mass filament skin (cf. G17 for the multiwavelength synthetic imaging). The molecular clouds, having the lowest volume filling, display the largest scatter. Globally, the condensed gas cannot be treated as ballistic or free-falling objects, as all phases participate in the hydrodynamical layer-within-layer cascade. Note that although some of the condensed gas can be accreted by the SMBH, the phases are continuously replenished by the turbulent condensation rain.

The bottom panel of Fig. 2 shows that, in the synthetic observations with a pencil beam (small aperture through the center), the velocity dispersion decreases significantly, down to a few 10% of the hot gas value. Therefore, we expect to detect commonly narrow lines with this approach, with FWHM down to a few 10 . At the same time, the scatter increases substantially (the distribution is lognormal with dispersion over all the phases of 0.41 dex), thus a broad component can also be present in pencil-beam measurements (e.g., by using absorption lines against the AGN backlight; §4.2). The broad component is typically associated with structures at large radial distance having small line shift. The narrower component tends instead to track the inner denser clouds (especially for the colder phases), which have experienced inelastic collisions in the nuclear region and are being funneled toward the SMBH. Such clouds can be better probed via major blue-/redshifted lines (§3.2.2).

The increased scatter (bottom versus top panel) reflects the chaotic and intermittent nature of turbulence, since the single warm/cold structures fill smaller volumes while condensing down the turbulence cascade. The pencil-beam approach indeed tends to sample a few or single clouds (e.g., for the molecular phase, the inner volume filling is 2%, gradually decreasing beyond kpc). Because of the turbulence inertial cascade, the velocity dispersion decreases121212By using the power spectrum analysis tool developed in Eckert:2017_PS, we checked that in projection the Kolmogorov cascade retains the same power index, in particular at small scales. as . From characteristic 2 kpc to 20 pc scales, this implies a factor of 0.2 decrease in velocity dispersion, as retrieved for the molecular phase in Fig. 2. Warmer phases are less compressed, thus having a lower decrement, as shown by the positive best-fit line slope.

The scatter related to multiple off-center pencil-beam measurements of the condensed gas line broadening can be used as a new way to quantify the level of small-scale intermittency in the turbulent medium. It is worth noting that the cold molecular phase suffers the largest scatter in line broadening. While the numerous cold clouds trace the condensation out of the peaks of the filamentary warm gas (in turn formed out of the turbulent hot halo), they also experience chaotic collisions and drag with all other phases, in particular at small distances from the SMBH. A fraction of the clouds may turn into young star clusters, decoupling via the collisionless dynamics (likely retaining the progenitor cold gas velocity dispersion). However, a significant in all the condensed gas phases implies that turbulent pressure is a key component (dominating over thermal pressure) which prevents most clouds from major collapse (cf. G17), in agreement with the low mean star formation rates and large cloud virial parameters () observed in early-type galaxies (e.g., David:2014; Temi:2017).

#### 3.2.2 Bulk/inflow motions: line shift

In Fig. 3, we analyze the bulk motions during the same CCA run via the mean velocity along the line of sight, or analogously, via the spectral line blue-/redshift (as offset from the systemic velocity). We note that our galaxy (stellar) systemic velocity is always null, as the simulation box is centered on the (static) gravitational potential. The ensemble gas detection (top panel) typically shows a small line shift with fairly contained scatter, slightly increasing toward the cooler phase. The logarithmic mean and dispersion of its magnitude over all the phases are . The line centroid would sometimes appear consistent with the galactic dynamics, given the typical measurement uncertainties. The pencil-beam detections (bottom panel) show instead a substantially larger line shift with global logarithmic mean and dispersion, . The fastest structures are typically inner clouds which have collided, canceling angular momentum, and are often falling toward the SMBH, within the Bondi capture radius. In both cases, we expect on average a similar fraction of blue- and redshifted lines (at least in emission), as clouds can drift in front of or behind the accretor.131313In absorption, a powerful outflow (on which some AGN surveys are selected) may skew the line absorption distribution toward the strongly blueshifted () side, swamping the infalling clouds. Conversely, during major angular momentum cancellation episodes, absorbed redshifted lines may be more frequent. In observations, separating inflow from outflow is non-trivial and multiple constraints are essential (e.g., reverberation mapping).

Interestingly, the condensed gas kinematics distributions are best fitted by lognormal distributions.141414 The average fractional difference (as absolute value) between the mean/1 and median/34.1% interval for the pencil-beam points in log space is just 2.6/7.5%, respectively. This is particularly evident when the range increases to several dex as for the velocity shift, since the linear approximation is no longer valid, and the high-end tail of the distribution becomes prominent. The reason is that turbulence continuously drives nonlinear perturbations in all thermodynamic properties with a characteristic lognormal shape (e.g., Gaspari:2014_coma2). The turbulence cascade (and related multiphase condensation) is indeed a multiplicative process, with smaller and smaller eddies generated within larger vortices, which can also be seen as a power-law inertial range in Fourier space.

Overall, as shown for the line-broadening kinematics, the adoption of a pencil beam leads to sampling smaller structures, thus tracing a lower velocity dispersion and a larger velocity shift of the infalling multiphase clouds and filaments. This method of probing the inner CCA via observations of narrow and significantly shifted lines is a simple and promising method, which is particularly efficient in absorption against a strong AGN backlight emitting in the band of the targeted gas phase (e.g., GHz radio for CO gas). An excellent case study is A2597 (Tremblay:2016), where ALMA detected 3 central infalling narrow-line clouds with redshifted velocities up to 4.2 for a large dataset comparison).

The previous simulations are part of our numerical campaign to dissect the multiphase physics of clusters, groups, and massive galaxies, as we test and disentangle each physical process step by step. We discuss below the main limitations of the current runs. In general, such extra – typically subdominant – physics tends to alter the details of the condensed gas (morphology, ionization layers, etc.), while the statistical thermodynamic and kinematic properties are expected to remain similar, i.e., the gas condenses via the turbulent top-down cascade, feeding the SMBH via accretion of clouds and filaments that are disordered on large scales.

Magnetic fields and cosmic rays (CRs) can provide non-thermal pressure, in addition to turbulence, and further support the condensation collapse, thus altering the size of the clouds and filaments (e.g., Sharma:2010). The draping provided by the large-scale -fields around the warm and cold structures is also expected to mitigate the heat and mass exchange between the different phases. On the other hand, the typical strength of magnetic fields is a few G, thus in the hot phase, implying that they are dynamically unimportant over most of the volume. Regarding CRs, their transport properties (as diffusion and streaming) out of the Galaxy are highly uncertain. Fermi telescope has also put severe upper limits on the amount of gamma-ray emission in the ICM, with no significant signal even in stacked analyses (e.g., Huber:2013), implying CR-to-thermal energy ratios of less than a few percent.

Speaking of diffusion processes, anisotropic conduction and viscosity will likely promote the formation of more elongated filamentary structures, with more equilibrated temperature and momentum along the magnetic lines. However, ram-pressure stripping observations (e.g., DeGrandi:2016; Eckert:2017_ramP), analytic studies (e.g., Burkert:2000), and plasma particle-in-cell simulations (e.g., Komarov:2016; Kumar:2017) suggest that the transport Spitzer/Braginskii coefficients are suppressed by at least one to two orders of magnitude due to plasma micro-instabilities (e.g., firehose and mirror) and/or highly tangled magnetic fields below the Coulomb mean free path.

Hot stars and AGN radiation can all contribute to ionize the external layers of the warm filaments. While the ionized skin depth varies widely depending on the clump density and location (cf. Valentini:2015), the typical heat input is modest in massive/early-type galaxies (see Loewenstein:1990), which have both very low star formation rates and Eddington ratios. Connected to this, stellar feedback is several orders of magnitude lower compared with AGN heating (e.g., Gaspari:2012b). Self-gravity is also negligible as most of the clouds have a large virial parameter because of non-thermal pressure, as found in ALMA data (e.g., David:2014; Temi:2017). Despite its simplicity, our feeding runs showed that radiative emission from line recombination during the condensation cascade down to  K, together with mixing with the hot plasma, can reproduce H maps similar to those retrieved by SOAR (Gaspari:2015_cca, Sec. 8). Regardless of the ionization level and microphysics, the turbulence cascade is a top-down process that develops in an analogous way, generating extended filaments and cloud associations that trace the self-similar turbulent eddies.

Overall, although more physics is likely at play in the ICM/IGrM – in addition to gas dynamics, gravity, radiative multiphase cooling, AGN heating, and turbulence implemented here – the current simulations already provide a robust framework and laboratory to assess the main statistical properties of the multiphase condensation cascade. Passing different multifrequency observational tests (e.g., radial profiles, surface brightness maps, cooling rates, emission measures, AGN variability; cf. G12-G17) gives confidence that the simulations are already capturing the main features of the real condensation process. Our forthcoming works will be important to dissect in depth the above physics and assess any deviation from the current results, thus improving the theoretical understanding of the multiphase gas precipitation.

\capstartfalse \capstarttrue

## 4. Observational Data and Comparison with Simulations

We present here new observational constraints – together with the available literature data – on the line broadening and line shift of the warm and cold gas kinematics in massive galaxies mainly within the cores of clusters and groups (spanning the mass range  - ). We compare them with the above numerical results, discussing the key differences between the ensemble and pencil-beam method, as well as the related limitations. At the same time, the following data points provide an estimate for turbulence velocities and local cloud kinematics, which can be used in subsequent analytical, numerical, or observational works. For instance, this can be used to calibrate and remove systematics of indirect methods estimating the gas kinematics, or to model non-thermal pressure support in semianalytic studies and subgrid models for large-scale simulations. We remark that the goal is to show the potential for such a kind of global analysis and not to delve into the details of each object, which is left to future work.

### 4.1. Ensemble-beam detections

As shown in §3.1-3.2.1, the primary method which allows us to retrieve the volume-filling turbulent velocity is to analyze the ensemble LOS velocity dispersion. In observations, this can be achieved by extracting a single, integrated spectrum over the maximally feasible radial aperture covering the entire warm ( -  K) or cold ( K) gas emission region. To show such capability, we applied the ensemble method to the Hamer:2016 (Hamer:2016 – H16) sample including 68 well-resolved H+[NII] objects – mostly galaxy cluster cores, plus 12 massive groups. An analogous method is applied to other 8 objects not included in the H16 sample and displaying cold/warm gas emission (e.g., McDonald:2012_Ha; Werner:2014; and the Perseus cluster, Gendron-Marsolais et al. in prep.). Table 1 lists the retrieved constraints and sample details for all the 76 objects.

The observational analysis to retrieve the gas FWHM () and line offset from the systemic velocity (shift) was carried out as follows – by taking the H16 sample as reference. Initially, the total spectra of the continuum and line-emitting regions are extracted from the data cube (e.g., VIMOS IFU for the H16 sample). For the line emission, this is done by taking the H flux maps (emission detected at S/N ). A masked cube is then created by discarding all spaxels in each wavelength channel with no emission. The remaining spaxels are summed to give a total flux value for that channel, producing a total spectrum for the line-emitting regions. Table 1 lists the extraction radius151515Although representative of the warm gas bulk emission with high S/N, in a few objects, may not match the full extent of the nebula due to the limited field of view and/or H absorption., with a median  kpc for clusters and 3 kpc for groups.

The total continuum spectrum is determined in the same way, using a collapsed white-light image in place of the H map and discarding spaxels where the continuum flux is less than 1/10 of the peak flux from the galaxy center161616This empirical threshold ensures that sufficient sky pixels are recovered to calculate the sky emission, which is then subtracted. We tested different thresholds (1/5 and 1/20) and found no significant difference in the retrieved velocity offsets.. The H masking is then inverted (discarding only spaxels with H detection) and the standard deviation of the remaining spaxels in each channel is calculated to provide an estimate of the error related to the two total spectra. The kinematics of the ionized gas and stellar components of the galaxy are then determined by fitting Gaussian profiles to the relevant total spectrum (a triplet fit to the H+[NII] complex for the line-emitting gas and a negative doublet fit to the NaD absorption171717We note NaD can have both stellar and gaseous origins; however, the NaD features we fit are all substantially broader than the emission lines, indicating they do not originate from the ionized gas in the galaxy and thus are most likely of stellar origin. for the stellar component) using a minimization procedure (see H16 for more details). Finally, the FWHM of the line-emitting gas is extracted directly from the fit181818 We tested an alternative approach, computing the ensemble velocity dispersion as the RMS of the projected velocity shifts from several patches within . This method introduces substantial noise due to patches with low signal. Moreover, each velocity shift has experienced positive/negative summation along LOS, so this RMS is typically a lower limit to the actual line broadening . We recommend to use the more robust single-spectrum FWHM. . The velocity difference between the H emission and the NaD absorption gives the velocity offset. It is worth noting that the total uncertainty is dominated by the systematic error (spectral sampling and instrumental line spread function), while random noise is drastically reduced due to the aggregation of several 100 spectra for each source (, with the number of spectra).

As a diagnostic tool to understand the kinematic properties with different methods, we propose a novel diagram confronting the line broadening versus the line shift magnitude. Fig. 4 shows the observed data points for warm and cold gas, which are compared with the simulation predictions (shaded contours). The shaded contours show the simulation bivariate distributions with mean and 1-3 standard deviations found in the previous run (§3.1), which are best fitted by lognormal distributions. As the §3.1 run probes a large dynamical range and varying cluster regimes, we use as reference line-broadening normalization the mean of its entire hot gas distribution, (which is comparable to that of the warm gas; Fig. 1), and show the logarithmic standard deviation from such a mean. The simulated global velocity offsets are discussed in §3.2.2.

The top panel in Fig. 4 shows that the ensemble detection – both in simulations and observations – cover a specific section of the diagram, namely the top-left region, which is the locus of substantial line broadening and relatively low velocity offsets. The scatter in velocity dispersion is mild due to the ensemble/single-spectrum approach, which decreases the statistical noise of single clouds and filaments. The log mean of the two distributions differs by 3% and 6% along the broadening and shift axis, respectively. The overlap within 2 is evident to the naked eye, with essentially a null correlation angle and analogous symmetry. Along the broadening/shift direction, the observed data have a mildly larger standard deviation (5%/6%). Given the limitations of the models and observational biases, we deem the simulations and observations to be in good agreement. E.g., the observed velocity offsets can be very sensitive to the redshift measurement of the host galaxy, with systematics not always captured. Further, coherent warm gas structures dominating the field of view have the tendency to increase the LW velocity offsets. Running a 2D Kolmogorov-Smirnov (KS)191919 Following NRec:1992 based on Fasano:1987. The synthetic data sample is randomly extracted from the simulation distribution with the same number of data points as that of the observational sample (excluding the few upper limits) and bootstrapped 1000 times to give a mean -value. As with any single likelihood value, this should be taken with a grain of salt. First, the KS statistic in 2D is only approximated (as cumulative distribution functions are not well-defined in more than 1D). Second, the number of current observed points is limited. Perhaps more important, a single value attached to a comparison is reductive of the complexity included in either the simulation or observational measurement. two-sample test results in a -value of 0.02, implying that the null hypothesis cannot be rejected at the customary confidence levels of 99% and 99.9%. This corroborates the visual inspection that the two datasets do not deviate by a large amount but are perhaps not drawn from an identical distribution. Nevertheless, perfect equivalence shall not be expected considering the detailed differences between the synthetic and real value measurements.

The above analysis benefiting from a large (76) sample size shows that turbulence in galaxy cluster/group cores is contained within a relatively narrow window,  - 250, which implies subsonic Mach numbers  - 0.3. This corroborates indirect X-ray hot gas estimates (§1); e.g., Hofmann:2016 and Zhuravleva:2017 find a similar range via plasma density fluctuations in the core, with within 50  from our values (interesting cases are A2052, A85, and A1795), although their uncertainties remain substantial because of the density contamination tied to substructures. The upper and lower limits set via XMM-RGS (e.g., Ogorzalek:2017) and cosmological simulations (e.g., Lau:2009; Lau:2017; Nagai:2013) further support such range of subsonic ICM/IGrM turbulence. Last but not least, it is remarkable that the SITELLE warm gas constraint matches the direct high-resolution observed by Hitomi in the archetypal galaxy cluster Perseus (§3.1). While waiting for the next-generation X-ray instruments (e.g., XARM – the successor of Hitomi – and Athena; Ettori:2013) to give us high-precision turbulence detections, this analysis allows us to use warm and cold gas velocity dispersions as robust proxies for the hot gas turbulent velocities in large samples (note that future X-ray instruments will still require several days of exposure for just one object).

For five objects (Tab. 1) we report literature detections in more than one band (besides Perseus, which was tackled in §3.1). Except for NGC 5846, the four other galaxies have ensemble broadening comparable between the cold and warm gas within 50 , with NGC 5044 and NGC 6868 two exemplary cases. The velocity offsets are also broadly aligned. On the other hand, the currently available extraction regions are often dissimilar; indeed, we observe larger values, all associated with larger extraction regions. Upcoming systematic multiwavelength investigations – which our team is currently undertaking with Chandra, XMM, ALMA, VLT, SOFIA, HST, Magellan, SOAR, SITELLE, and IRAM – will be key to calibrate the multiphase kinematics over the same extraction region and with similar depth. At the present, the ensemble H[NII] nebulae appear to be the best channel to test the volume-filling turbulence and related velocity dispersion.

\capstartfalse \capstarttrue

### 4.2. Pencil-beam detections

For the pencil-beam (small aperture, below a few arcsec) detections in massive galaxies, we report the published value from the literature and a few new detections. A larger sample requires new observational programs (e.g., one recently approved in ALMA Cycle 5 – PI: A. Edge). The most used pencil-beam approach is to detect absorption lines (e.g., due to HI and CO clouds having significant optical depth) against the (radio) AGN, which acts as a backlight source. The BCG stellar light can also be used as backlight, in combination with NaD absorption. The analysis procedure is then analogous to that above, extracting the systemic velocity offset and FWHM from Gaussian fitting of the absorption lines in the continuum-subtracted spectrum. Hogan:2014 discusses in detail the observational reduction with different instruments such as VLA and WRST (see Tremblay:2016 for ALMA data). Although more challenging due to the low S/N in massive galaxies, small-aperture observations can be used to track small-scale clouds via emission features such as CO(2-1). Table 2 lists the 47 fiducial detections for the 19 available objects, with the addition of the mean properties of the Maccagni:2017 (Maccagni:2017 – M17) radio galaxy sample and of the Perseus IRAM regions.

Fig. 4 (bottom panel) shows the observational detections as red points (identical non-circle symbols denote the same host galaxy), compared with the simulation results for all the condensed gas with 1–3 confidence intervals. For the pencil-beam approach, the §3.2 run covers a meaningful range to statistically test the broadening distribution; the reference mean is the same as in the previous section, but now the condensed gas has a ratio lower than 1:1 due to the pencil-beam sampling (Fig. 2, bottom). The simulated velocity offsets over all of the condensed gas are discussed in §3.2.2.

As anticipated by the numerical analysis, it is clear that this method substantially boosts the logarithmic scatter in both the broadening and shift dimensions up to nearly 0.5 dex. The mean line broadening () is significantly lower than the ensemble measurement, as predicted by the simulations. The mean velocity offset has instead larger values . This is because the pencil beam is sampling smaller and typically fewer condensed elements. About 1/3202020Multicomponent systems should be much easier to observe with the more modern receivers (e.g., CABB on ATCA), which combine a wide bandwidth with high resolution. of the central galaxies observed in absorption display a dual component: either with a large line broadening and small shift or with a large shift and fairly contained broadening. A mild anticorrelation appears to be present, with the faster (and denser) clouds associated with the nuclear inflow toward the SMBH sink region. It is interesting to observe that the broad component in absorption is often a reasonable proxy for the ensemble velocity dispersion, as it tends to sample multiple clouds along the LOS. For instance, the pencil-beam broad components of A2390, A1795, A2597, and N1275 all reside within from the actual ensemble value, although there are exceptions, such as Hydra-A (which has an abnormal 5 kpc disk).

In Fig. 4 (bottom), we also plot the mean and deviation of the HI absorption detections in the M17 sample of 66 radio galaxies (yellow bars). Although comprising mostly non-central radio/elliptical galaxies and requiring deeper follow-up observations for each target, the HI absorption broadening and shift are consistent with the above central galaxies and simulation properties, with a slightly lower average line shift. An interesting case study is the early-type galaxy PKS 1718, where Maccagni et al. (in prep.) find redshifted clouds in HI, H, and CO clouds infalling within the Bondi radius with velocities up to , as found very similarly in A2597 central galaxy (Tremblay:2016). Other interesting cases are PKS 1740 (Allison:2015), PKS 1657 (Moss:2017), and NGC 3998 (Devereux:2018). This suggests that the top-down condensation and CCA are also likely common phenomena in low-mass/non-central galaxies, given that plasma atmospheres are expected to be ubiquitous (e.g., Anderson:2015) and subsonically perturbed by any AGN type (radio, quasar, etc.) and/or cosmic flow.

The pencil-beam method can be further used in emission and off center. For instance, taking advantage of ALMA high angular resolution and CO sensitivity, we probed several giant molecular associations in the massive galaxies NGC 4636, NGC 5846, and NGC 5044 (Temi:2017). Remarkably, the emission features show a similar mean and scatter to the above HI absorption features. Commonly, small-scale clouds with low broadening () are associated with large velocity offsets above 100 . At the same time, the inspected masses, radii, and cospatiality of the giant molecular associations are consistent with the G17 simulation, corroborating the incidence of in-situ cooling via the multiphase cascade. Line absorption against the stellar light of a background galaxy is another promising way to retrieve the kinematics of gas at varying clustercentric distances (e.g., Smith:2017). Notice the beam must be small to achieve a proper pencil-beam detection (below the kpc or a few arcsec scale), otherwise the measurement will tend toward the ensemble approach; the IRAM data of Perseus (orange) is an example of intermediate regime with a beam of 4 kpc (12 arcsec).

Focusing on the comparison between the observed and simulated distributions in the diagram, the mean differs by 2% along both axes. The fairly good overlap within 2 is again evident to the naked eye. The observed data have mildly larger/lower standard deviation (4%/10%) along the broadening/shift direction, respectively. Running a KS two-sample test (keeping in mind the limitations discussed in §4.1) results in a -value of 0.14 (0.18 including the M17 and IRAM points). The absence of low -values implies that the null hypothesis cannot be rejected, even at the less significant 95% level. The point to take away is that the two distributions are similar, though not necessarily identical. In particular, the simulation shows a milder anticorrelation (0.54 rad difference in angle rotation), although the scarcity of observed points in the bottom left section of the diagram may be attributed to the difficulty in detecting both small shift and narrow lines in the spectra.

Overall, granted that the sample requires larger statistics, both simulations and data agree well on the same picture: the pencil-beam method is useful to track the inner infalling clouds (selecting the narrow features with a large velocity shift) or to have a preliminary estimate of the large-scale chaotic motions (selecting the broad features). The wide range of velocity shifts is key to allowing a clean separation of such components in the energy spectrum (e.g., in the radio band). The ensemble measurement instead provides a robust and direct constraint on the volume-filling turbulent motions, as the ensemble condensed elements actively participate in the large-scale kinematics via the top-down multiphase condensation cascade. The combination of the two approaches provides a powerful complementary diagnostic of the global and local gas kinematics.

## 5. The nonlinear multiphase condensation criterion

In the previous sections, we showed how the kinematics of the different phases is tightly related during the top-down condensation cascade and how we can convert between the ensemble velocity dispersions, in particular to that associated with the hot plasma. Here we discuss a key application of the retrieved turbulent velocities aimed to further understand the condensation process and the above multiwavelength observations.

### 5.1. Previous condensation criteria

A majorly debated topic in recent literature concerns the dominant criterion governing the formation of the condensed structures in the ICM/IGrM. Previous studies proposed that the ratio of the plasma cooling time to free-fall time212121Defined as , where is the gravitational acceleration due to the total mass within a radius . falling below a few tens is the triggering threshold of thermal instability (a.k.a. TI ratio). The mean value of the TI ratio is highly uncertain. McCourt:2012 and Sharma:2012 propose a value of 10. Gaspari:2012a simulations show a TI ratio between 8 and 25 (see their Fig. 3, bottom-right panel), which is similarly adopted by Voit:2015_nat. Valentini:2015 find TI-ratio thresholds that can reach a value of 70 in early-type galaxies. Taken at face value and over different studies, the uncertainty on the TI ratio can be an order of magnitude. It is important to note that this theoretical uncertainty accumulates on top of the intrinsic scatter due to system-to-system variations mainly tied to the cooling time (more in §5.2).

Four key issues arise with the TI ratio. First, the unity problem: why is the condensation occurring at such an elevated threshold, well above 1, which should be instead the physically sound transition222222For example, in a steady spherical cooling flow, the Lagrangian entropy equation can be written as (where ), i.e., as the r.h.s. ratio crosses below 1, condensation ensues (similar relations apply for the other conservation equations).? Second, the large 1 dex theoretical uncertainty in the actual threshold – even for similar systems – hints at the fact that the free-fall time is not the primary timescale of the condensation process. Third, observed condensed clouds do not follow ballistic orbits but are drifting with subvirial velocities (e.g., McNamara:2014; Russell:2016). Last but not least, there is the major observational hurdle of retrieving the free-fall times, i.e., the total masses, including the galactic and cluster baryonic and dark matter masses, which are notoriously challenging to constrain.

Another suggested criterion (Hogan:2017) is to consider the cooling time below a fixed threshold, e.g., 1 Gyr. Although empirically effective and less noisy, this has the limitation of being applicable only to some classes of objects, e.g., massive clusters, but not groups or galaxies (e.g., OSullivan:2017). Finally, the uplift action of the AGN cavity may significantly increase the gas total inflow time (), well above the free-fall timescale, and bring it to altitudes where (e.g., McNamara:2016). However, it appears too restrictive to consider just one specific directional (bipolar) and short stage of the feedback cycle, namely the trailing phase. The full process includes bubble inflation, uplift, and cocoon shocks, with isotropic turbulence being generated as a common by-product (this is true even if the driver is a merger or sloshing event). The turbulence time described below can be considered as a more universal concept of the reduced inflow time (with a clear quantitative observable), i.e., gas drifts in the halo regardless of the radial direction (inflow or outflow) and location (behind or around cavities).

### 5.2. The new C-ratio criterion

Given the results in §3, namely the cold and warm gas following the progenitor chaotic kinematics, we show in Fig. 5 that a robust and physically motivated condensation criterion is the ratio of the cooling time to the turbulence eddy turnover time , which marks the multiphase state of the cluster/group core. The condensation criterion can be alternatively written and interpreted as the velocity ratio , similar to a new dimensionless ‘Mach’ number, where is the cooling velocity.

Figure 5 shows the two characteristic timescales and the locus of extended H+[NII] filaments, for clusters (left) and groups (right). In the left panel, the blue curves represent the plasma cooling time of observed non-cool-core (NCC) clusters (Cavagnolo:2008) and cool-core (CC) clusters (following the recently updated constraints by Panagoulia:2014 and Hogan:2017). The CC profiles are best fitted by a broken power law: at large radii following the cosmic adiabatic baseline (Tozzi:2001) with (for a typical 5 keV cluster,  Mpc and ; Sun:2009a), while at small radii . The (X-ray) plasma cooling time is

 tcool=3kbTneΛ, (4)

where and is the plasma cooling function (Sutherland:1993; for clusters). The observed profiles are typically contained within 0.25 dex (blue bands; % confidence level). The 0.5 - 1 Gyr dotted black band crudely separates the observed clusters hosting (below) or not hosting (above) extended warm filaments traced via H+[NII] emission (e.g., Hogan:2017).

The eddy turnover timescale, i.e., the time which a turbulent vortex requires to gyrate – thus producing density fluctuations in a stratified halo – is

 teddy=2πr2/3L1/3σv,L, (5)

where is the velocity dispersion at the injection scale and using the Kolmogorov cascade appropriate for subsonic turbulence (Gaspari:2013_coma; we use here the approximation that smaller radii contain eddies of smaller size, ).232323The eddy turnover time differs from the turbulence dissipation time, which is 10 - 30 longer for subsonic turbulence, . At this level, the fluctuation generation (and related condensation) outweighs the slow turbulent heating. Turbulence is mainly injected via AGN feedback in the core. While it was previously difficult to assess the plasma kinematics, the ensemble measurement makes it now trivial to apply the turbulence velocity dispersion derived from one of the condensed phases. We use the mean three-dimensional () ensemble velocity dispersion, , found in the AGN feedback simulation and similarly in the observational sample (Fig. 4). We remark the ensemble dispersion is comparable in all phases, cold, warm, or hot. The injection scale can be traced by the diameter of the bubble inflated by the AGN outflows/jets, which start to decelerate as they deposit the kinetic energy (or, alternatively, from the size of the warm gas nebula; more below).242424 We note for the initial, brief cavity trailing stage and partial volume-filling turbulence (e.g., Brighenti:2015), the criterion should be applied in sectors (of size ) rather than azimuthally. The criterion can also be applied for the rarer merger or sloshing events, adopting their maximum diameter as the injection scale. From the large high-quality sample of Shin:2016, observed 5 keV cluster bubbles have a typical diameter  kpc. Given the weak dependence on the injection scale (), the eddy time scatter for a given system is driven by , which has a 90% confidence level of 0.2 dex (Fig. 4).

The key result from Fig. 5 is that, despite its simplicity, the crossing of and traces the condensation region. For our analyzed galaxy clusters (Tab. 1), the logarithmic median and dispersion of the extent radius () are , i.e., a 2 range of 3.5 - 52.5 kpc (which is close to that found by McDonald:2010 in a smaller sample of CC clusters with detected warm filaments; dashed ellipse). This range is tracked by our proposed dimensionless cooling number (see the overlapping blue and red bands). Around such a threshold, turbulence in a stratified halo drives significant density/entropy perturbations in the hot phase, (cf. Gaspari:2014_coma2); the overdense hot gas rapidly cools down, promoting the multiphase cascade down to warm filaments and molecular clouds. The cospatiality between different phases, in particular the soft X-ray and radial extent of H gas, corroborates this scenario (e.g., Hogan:2017b). We remark that this process differs from linear thermal instability models, which assume tiny perturbations growing exponentially against the restoring buoyancy force252525The related timescale is the free-fall time (defining the TI ratio), since the buoyancy (Brunt-Väisälä) frequency is for non-isentropic CC clusters. in a heated halo (indeed, CCA can develop even in a non-heated atmosphere, given significant turbulent fluctuations).

An important result is that within the condensation region, both timescales do not deviate drastically from each other (), except in the very inner region, implying that both generation of fluctuations and cooling act at fairly concurrent times, inducing the drop-out of warm filaments with multiple scales within the extent radius (a fraction of which triggers the central AGN). In other words, naturally traces what other studies refer to as ‘precipitation’-feedback balance threshold. If the plasma cooling timescale were too short with perturbations that are weak and/or injected at very large scales (), a monolithic overcooling of the X-ray halo would develop without extended multiphase structures. This is expected to be infrequent (e.g., A2029 and A2107) because of the AGN feedback self-regulation which elevates back while decreasing via the injection of turbulence. Conversely, when the cooling time is too long262626If the injection scale stays the same (e.g., without any complementary cosmic-weather turbulence; Lau:2017) the red curve should flatten out, hence the dashed prolongation in Fig. 5. In any case, at large radii is expected to stay above this line. (), even significant perturbations (e.g., driven by mergers) cannot induce multiphase condensation, thus preserving the NCC structure (an exemplary case is the Coma cluster).

The right panel in Fig. 5 shows the same timescales as above but for 1 keV groups/massive galaxies. Typical massive groups have and at  kpc (Sun:2009a). Cluster-like NCC groups with a flat and elevated entropy core seem nonexistent (although a more complete sample is required; see also OSullivan:2017). Following Sun:2009a constraints, we plot as ‘NCC’ the upper envelope of their observed groups. At large radii, the profiles follow the adiabatic baseline, while within the CC groups trace a similar scaling to that found by Panagoulia:2014. Given the non-flattening entropy profiles of groups, the H filament threshold is best taken within 5 - 15 kpc and not as a straight line, where multiphase groups show significantly lower entropy corresponding to  100 Myr (Werner:2014). Regarding the eddy time, following the observational constraints by Shin:2016, the driven AGN bubbles in groups show smaller average injection scales with diameters  kpc. Because of the self-regulated, thus diminished, AGN feedback power (in part counterbalanced by the smaller deposition volume), the average turbulent velocity dispersions are lower too. For our H+[NII] emitting groups, the median is  20% lower, so we use a reference 3D velocity dispersion (see also Ogorzalek:2017).

The same main result for clusters applies to groups, but with the condensation region () now more compressed between radii of 0.6 - 15 kpc. For our group sample, we find a logarithmic median and dispersion for the extent radius of , i.e., a 2 range of 0.9 - 8.3 kpc. McDonald:2011a find a similar median around a few kpc, with a slightly wider range, 0.5 - 18 kpc (dashed ellipse). Unlike clusters, some massive groups without extended multiphase filaments (e.g., NGC 4472 and NGC 1399) and residing at the CC-NCC transition can still have low cooling times (below 50 Myr) within the inner  pc, crossing below the eddy time. This is likely related to the ubiquitous presence in groups of central kpc-scale warm ‘coronae’ (Sun:2009b).

It is interesting to point out that the condensation radius is often comparable to the AGN bubble radius, and thus to the injection scale (again a sign of tight self-regulated feedback), for both clusters and groups. Therefore, if an observation is lacking any evident cavity detection and if the goal is to quickly estimate the ratio, then can be used as an alternative estimate for the injection scale. This keeps the calculation solely based on the condensed gas properties. We remark that the ratio can be purely retrieved from observational data (e.g., H+[NII]) without requiring velocity dispersions or injection scales from simulations (though the two have been proven to be consistent; Fig. 4).

Addressing the key issues introduced in §5.1, the criterion is able to solve the unity problem, which should be the natural transition of physical processes. It also does not show the one order of magnitude theoretical uncertainty of the TI ratio (which has possible threshold values 8 - 70), though still retaining the dex intrinsic scatter due to system-to-system variations mainly tied to the cooling time. Perhaps more important, the eddy time is observationally much easier and robust to measure compared with the challenging total masses, as the large-scale velocity dispersion can be computed in several low- phases via the ensemble method (§3.2). Indirect hydrostatic masses are required to compute ; however, solely considering turbulence and AGN feedback272727 Geometrical deprojection biases, the AGN contamination, and lack of X-ray spectral resolution to constrain central temperature gradients are additional issues affecting the retrieval of hydrostatic masses in the core region., hydrostatic masses in the core can be off by a factor of several (e.g., Gaspari:2011a), inducing a major systematic uncertainty in the TI ratio. On the other hand, the criterion takes advantage of the substantial noise reduction provided by the integrated spectrum measure, which is directly accessible from multiple frequency bands.

We note is a lower limit to the eddy time, as the condensed structures do not escape the cluster or group central potential. Moreover, AGN feedback self-regulates based on the halo X-ray luminosity (thus mass), so we still expect some degree of secondary correlation between these two dynamical timescales. Finally, while the hot gas cooling time shows the largest variation between CC and NCC systems, it is important to note that the physical crossing threshold is still set by . The crossing gradually decreases from massive clusters to groups and small galaxies, from 0.5 Gyr to  100 Myr, as the eddy time has a total mass trend of roughly .

In conclusion, we expect groups to display condensed structures more frequently than in clusters, but with a much more concentrated topology and with lower warm and cold gas masses. An excellent case study for the presence of soft X-ray turbulence/perturbations, cospatial warm filaments, and cold molecular gas is NGC 5044 (e.g., Gastaldello:2009; David:2017). Over different environments, the  - 1.8 range (90% interval) probes well the condensation region of observed systems and thus the multiphase state of the core. Future studies should extend the sample of clusters, groups, and especially low-mass galaxies, taking advantage of the linked kinematic properties between the hot, warm, and cold gas phases, and thus complementing the thermodynamical constraints set via high-resolution CCD imaging.

## 6. Conclusions

We probed the kinematic tracers of the top-down multiphase condensation cascade in both long-term AGN feedback simulations and high-resolution chaotic cold accretion runs. We complemented the theoretical predictions with new observational constraints of the proposed novel methodology, together with literature data. Our main results are summarized as follows.

• In long-term (AGN feedback) and short-term (CCA feeding) simulations, we find evidence of a tight correlation in the LOS between the hot X-ray phase and the condensed phases as ensemble (wide-aperture beam), allowing us to convert between different tracers in the UV, optical/IR, and radio bands. The RMS scatter from the linear best fit is  14%. The tight kinematics is corroborated by the direct detections of Hitomi (X-ray) and SITELLE (optical) in the Perseus cluster.

• As ensemble, the multiphase condensed structures display significant LOS velocity dispersion and mild bulk velocities. The pencil beam measurement (e.g., absorption against the AGN backlight) instead samples fewer and smaller condensed structures, substantially increasing the scatter due to the turbulence intermittency: the average velocity dispersion decreases following the eddy cascade, while the bulk velocity can reach values up to several 100 km s.

• We presented new observational constraints for over 70 clusters and groups of the warm ( K) and cold ( K) gas kinematics for the ensemble detection (Tab. 1), which can be used as reliable proxies for the turbulent velocities, especially for the challenging X-ray plasma (which would take an exposure of several 100 ks per nearby object even for XARM or Athena). Comparing the lognormal distributions, the simulation predictions and observations are consistent and show a range  - , with mean in cluster cores.

• A novel diagnostic diagram of the (logarithmic) spectral line broadening versus line shift discriminates among the different kinematics and related scales. Both simulations and observations indicate that, while the ensemble points are confined in the upper-left region, the pencil-beam (small aperture, a few arcsec) detections can show a dual broad and narrow component sampling the chaotic large-scale gas or small-scale clouds falling toward the SMBH, respectively.

• We showed that a new nonlinear multiphase condensation criterion (facilitated by the ensemble conversion), i.e., the ratio of the cooling time and eddy turnover time (with a 90% interval of  dex) marks the condensation extent region, as shown by the warm gas observations in the cores of groups and clusters. Besides solving the unity threshold problem, the ratio can be used to assess the multiphase state of the system and is a much more robust and direct observational quantity to measure – via the ensemble, single spectrum method – compared with the challenging total masses of the linear thermal instability ratio .

This study highlights the importance of undertaking multiwavelength campaigns (e.g., combining Chandra, XMM, ALMA, VLT, HST, Magellan, SITELLE, SOAR, IRAM, and MUSE), some of which our team is currently pursuing. The combination of the ensemble and pencil-beam method provides a powerful complementary diagnostic of the global and local multiphase gas kinematics, which can be optimally leveraged by the available and future IFU and spectroscopical collections of data. At the same time, this work highlights the fairly unexplored potential of joint numerical and observational studies of multiphase gas. While future observations will expand the sample size, allowing more accurate statistics on the multiphase kinematic tracers, thanks also to new facilities (e.g., XARM, Athena, JWST, SKA, and CARMA2), more advanced simulations with additional physics and an upgraded dynamical range will be instrumental to further shed light on the formation and evolution of multiphase gas in galaxies, groups, and clusters of galaxies.

## Acknowledgements

M.G. is supported by NASA through Einstein Postdoctoral Fellowship Award Number PF5-160137 issued by the Chandra X-ray Observatory Center, which is operated by the SAO for and on behalf of NASA under contract NAS8-03060. Support for this work was also provided by Chandra GO7-18121X. S.L.H. acknowledges support from the ERC for Advanced Grant Program number 339659 (MUSICOS). J.H.-L. is supported by NSERC through the discovery grant and Canada Research Chair programs. M.G.-M. is supported by NSERC through the NSERC Postgraduate Scholarships-Doctoral Program (PGS D). N.W. is supported by the Lendület LP2016-11 grant awarded by the Hungarian Academy of Sciences. A.C.E. acknowledges support from STFC grant ST/P00541/1. S.E. acknowledges support from ASI-INAF I/009/10/0, NARO15 ASI-INAF I/037/12/0, and ASI 2015-046-R.0. S.P. is supported by the Spanish Ministerio de Economía y Competitividad (MINECO, AYA2013-48226-C3-2-P, AYA2016-77237-C3-3-P) and the Generalitat Valenciana (GVACOMP2015-227). FLASH code was in part developed by the DOE NNSA-ASC OASCR Flash Centre at the University of Chicago. HPC resources were in part provided by the NASA/Ames HEC Program (SMD-16-7320, SMD-16-7321, SMD-16-7305; SMD-16-7251). We thank R. Morganti, F. Maccagni, M. Voit, F. Combes, B. McNamara, C. Feruglio, M. Donahue, Y. Cavecchi, V. Moss, N. Devereux, and the anonymous referee for several insightful discussions and constructive feedback, which helped to improve the manuscript.

## References

You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters