© ESO 2019

1. Introduction

Characterizing the circumgalactic medium (CGM) and understanding its role in galaxy evolution is a key goal of current astrophysical models. Usually defined as the gas between the central galaxy and the surrounding intergalactic medium (IGM), this multi-phase medium has been observed for decades in both emission and absorption from the X-ray (e.g., Anderson et al. 2013; Li et al. 2017, 2018) to the UV-optical bands (Bordoloi et al. 2011; Keeney et al. 2013; Kacprzak et al. 2014; Borthakur et al. 2015) and is flowing both inward and outward in galaxies, as suggested by numerical simulations (Ford et al. 2014; Suresh et al. 2015; van de Voort & Schaye 2012). Despite the great theoretical and observational effort, however, the nature, origin and fate of this elusive medium are still much debated.

For a long time the most studied phase has been the hot CGM, also called the galaxy corona, with temperatures of T ∼ 10^6 − 7 K, predicted decades ago by classical cosmological models (e.g., White & Rees 1978; White & Frenk 1991, but see also Binney 2004) as gas heated by shocks to the galaxy virial temperature. Although the observed luminosities are lower than predicted, hot gas in galaxy halos has been observed through X-ray observations both around early-type (Forman et al. 1979; Jones et al. 2002; Mathews & Brighenti 2003; Bogdán & Gilfanov 2011) and late-type galaxies (Dai et al. 2012; Li et al. 2017). The total amount of mass in this hot gas phase is however still uncertain (e.g., Anderson & Bregman 2011; Bogdán et al. 2013; Li et al. 2018).

Recently, a number of investigations (Stocke et al. 2006; Rudie et al. 2012; Werk et al. 2013; Tumlinson et al. 2013) have also focused on the cooler (T <  10⁵ K), low-density, ionized gas in the galaxy halos, originally discovered by Boksenberg & Sargent (1978) and Bergeron (1986) and observed using UV and optical absorption lines in the spectra of background quasars. The presence of this cool gas seems ubiquitous around both external galaxies and the Milky Way, with high covering fractions (e.g., Shull et al. 2009; Lehner & Howk 2011; Werk et al. 2014; Borthakur et al. 2016) and potentially large masses (Werk et al. 2014; Stocke et al. 2013; Stern et al. 2016). Understanding this cool medium is crucial in the study of the galaxy evolution and in solving the galaxy missing baryons problem (McGaugh et al. 2010). So far our knowledge of this medium is largely limited and the general picture is still unclear, especially regarding its origin and dynamics. For star-forming galaxies, both transverse absorption-line studies and down-the-barrel observations of inflows (e.g., Rubin et al. 2012; Lehner et al. 2013; Bouché et al. 2013; Borthakur et al. 2015) and outflows (e.g., Rubin et al. 2014) suggest that the observed CGM dynamics is consistent with the recycling scenario (Ford et al. 2014) in which the galaxy central enriched outflows are not able to escape from the potential well and then fall down again onto the galaxy. In this picture, the CGM is formed by a continuous cycle throughout the halo similar to the galactic fountain scenario proposed for smaller, galactic scales (Shapiro & Field 1976; Bregman 1980; Fraternali & Binney 2006).

The situation becomes however more complicated for passive early-type galaxies (ETGs). Cool ionized gas is in fact observed also in the CGM of passive massive galaxies, through observations of MgII (Gauthier et al. 2009; Huang et al. 2016), HI and other low-ion absorption lines (Thom et al. 2012; Chen et al. 2018). The presence of cool-enriched gas around quiescent galaxies remains a puzzle. In fact, in the recycling scenario described above, a fundamental role is played by the central galaxy star formation, which causes the outflows thus starting the entire cycle. However, it is not clear how we can explain the presence of cool gas in the halos of galaxies with little or no ongoing star formation. The mechanisms that can produce the cool gas and prevent it from fueling star formation in the central passive galaxy are also unclear. The main goal of this work is to address these questions.

To this purpose, we use the observations of the COS-LRG collaboration (Chen et al. 2018; Zahedy et al. 2019) who, using the Cosmic Origin Spectorgraph (COS, Froning & Green 2009) on the Hubble Space Telescope (HST), has observed the cool CGM in the halos of 16 luminous red galaxies (LRGs), massive elliptical galaxies at z ∼ 0.5, where in principle cool gas is not expected to exist. One of the features of these observations is that the cool CGM is not homogeneous, but composed of different absorbing clouds, with different velocities, all bound to the central galaxies, a property of this medium that was already well known from previous absorption-line studies (e.g COS-Halos survey, see Werk et al. 2013, 2014; Tumlinson et al. 2013). Our purpose is to model the observed cool gas kinematics. Due to the absence of strong outflows in these galaxies, we assume that the cool CGM clouds are infalling toward the central galaxies.

As we will discuss in detail (Sect. 3), the dynamics of the CGM – and in particular the interactions between its different phases – depends on a number of physical processes acting on both very large and very small physical scales. Ideally one would like to model the entire CGM with large-volume high-resolution hydrodynamical simulations. Unfortunately, this is still far out of the current computational capabilities: in fact, even the current “zoom-in” simulations (e.g., Grand et al. 2017; Fattahi et al. 2016), which can trace a single galaxy halo, do not have enough resolution to resolve the small-scale CGM interactions (Armillotta et al. 2017). We therefore adopt a semi-analytical modelling, keeping track of the gravitational force which dominates the dynamics on large scales, as well as taking into account prescriptions on hydrodynamical processes from small-scale high resolution simulations (e.g., Heitsch & Putman 2009; McCourt et al. 2015; Armillotta et al. 2016; Grønnow et al. 2018).

This paper is organized as follows: in Sect. 2 we briefly describe the COS-LRG observations, their characteristics and the main observables that we use as a comparison for our models; in Sects. 3 and 4 we describe the construction of the model and the comparison with observations; in Sect. 5 we present our results while in Sect. 6 we discuss the implications of our findings for the properties, origin and fate of the cool CGM around massive ETGs; finally, Sect. 7 contains the summary and the conclusions of the work.

2. ETGs sample and observational constraints

In this work we focus on the data of the COS-LRG program. Here we briefly describe the main characteristics of these observations, while we refer, for further information and details, to the COS-LRG papers (Chen et al. 2018; Zahedy et al. 2019). This survey provides high-resolution kinematic data of the cool CGM around 16 massive early-type galaxies and it is therefore the most suitable for the purpose of this work, which aims to reproduce in particular the kinematics of the cool CGM clouds around passive galaxies. We assume that, due to the similarity of the objects in this sample, the CGM has the same behavior and properties in all these galaxy halos, modulo some scaling factors due to the slightly different halo masses and virial radii, which we will account for in our model. The main properties of the galaxies in the sample are reported in Table 1.

The CGM is characterized using HI and metal absorption lines in the COS and optical echelle spectra of UV-bright quasars with impact parameters (the projected distances between the quasars and the galaxies) R <  160 kpc. A typical feature of these observations is the presence, for each line of sight, of different velocity components for the same ionic transition, characteristic of a medium that is not uniform, but composed of different clouds with different projected velocities, plausibly different intrinsic positions along the line of sight and possibly different physical properties (e.g., density or ionized fraction). This is studied in Chen et al. (2018) and Zahedy et al. (2019) fitting Voigt profiles and obtaining the velocity offsets of the single components from the galaxy systemic velocity Δv. In this paper we use in particular the analysis carried out by Zahedy et al. (2019), who fitted the Voigt profiles imposing the same kinematic structure among HI, low-ionization and intermediate-ionization species (including MgII, SiII, SiIII and CIII). This choice is justified by the kinematic agreement found by Chen et al. (2018) in the absorption profiles of these different species. Each velocity component found by Zahedy et al. (2019) represents a single cloud: they find 42 clouds among all the galaxies (including three sightlines with no detections of absorption lines), with an average of 2.6 clouds per line of sight. In Fig. 1 we report all the 42 observed cool cloud velocities (Δv, relative to the galaxy systemic velocity) together in one single distribution. The two main observational constraints that we used for the comparison with our model results are then given by:

Fig. 1. Orange: observed velocities of the cool CGM clouds relative to systemic velocities of the galaxies, obtained using the COS-LRG data for a sample of 16 massive passive galaxies (their properties are reported in Table 1). The light-blue dashed line is a gaussian with a dispersion equal to the typical 1D (line-of-sight) virial velocity dispersion value in the sample.

1. the cloud velocity distribution;

2. the total number of observed clouds.

Due to the high resolution of the COS spectra, these kinematic data are very accurate.

Zahedy et al. (2019) also performed a photo-ionization analysis for each kinematic component, using the code CLOUDY v13.03 (Ferland et al. 2013). This analysis allowed them to estimate the values of the total-hydrogen column densities N_H for each line of sight (reported in Table 1). Differently from the very accurate and reliable kinematic data, these values are however subject to various model-dependent uncertainties and therefore we do not use them as an observational constraint for our model fitting. We will discuss in Sect. 5 the comparison between these data and our final results.

For the purpose of this work, we are interested in the virial masses and virial radii of the galaxies in the COS-LRG sample, which allow us to calculate the halo virial temperatures. To calculate these quantities we first obtained M₂₀₀ using the stellar-to-halo mass relation (SHMR) of Kravtsov et al. (2018). We note that at these high stellar masses the slope of the SHMR is uncertain and the galaxy halo masses could be higher using different prescriptions (e.g., Moster et al. 2013). However, the relation from Kravtsov et al. (2018), obtained with a recent stellar mass function based on improved photometry, is in agreement with estimates based on X-ray observations, weak lensing and satellite kinematics for objects with these high masses. Thus we consider this prescription the most reliable for the massive galaxies of our sample. Once we obtained the values of M₂₀₀, we calculated r₂₀₀ through the formula

$$\begin{array}{c}\hfill r_{200}={\left(\frac{M_{200}2G}{200H^2}\right)}^{1/3},\end{array}$$(1)

where H is the Hubble parameter (function of redshift, with H₀ = 70 km s⁻¹ Mpc⁻¹) for the standard cosmological model (Ω_m, 0 = 0.3, Ω_Λ, 0 = 0.7). Using these values, we calculated the virial masses M_vir and radii r_vir extrapolating for each halo the dark matter (DM) profile (we assumed a Navarro Frenk White profile, see Sect. 3.1) to the radius with virial overdensity, which is Δ times the critical density of the universe, with Δ calculated using the prescription in Bryan & Norman (1998) (Δ = 119 − 141, depending on the redshift of each galaxy).

Finally, the virial temperature is given by

$$\begin{array}{c}\hfill T_{\mathrm{vir}}=\frac{\mu m_{\mathrm{p}}G{M}_{\mathrm{vir}}}{2k_{\mathrm{B}}{r}_{\mathrm{vir}}},\end{array}$$(2)

where μ = 0.6 is the mean molecular weight for a hot, totally ionized gas (Sutherland & Dopita 1993), m_p is the proton mass, and k_B is the Boltzmann constant. Our estimates of these virial quantities for each galaxy are reported in Table 1.

As explained above, the kinematic distribution in Fig. 1 is our main constraint and our purpose is to create theoretical models that are able to reproduce a similar distribution, with a comparable shape and a comparable number of clouds. It is interesting to note from this observed distribution that the clouds show an unexpected narrow distribution (small Δv). Indeed, the line-of-sight virial velocity dispersion for a galaxy halo with a mass equal to the mean value of the virial masses of our sample (M_vir = 10^13.3 M_⊙) is equal to 245 km s⁻¹, while the velocity dispersion of the distribution in Fig. 1 is equal to only 147 km s⁻¹. If the cloud motion were related only to the dynamical mass of the galaxies we would expect the clouds to have a dispersion comparable to the virial one. The small observed velocity dispersion, instead, suggests that there should be some mechanism that slows down these absorbers during their motion throughout the halos (see also Huang et al. 2016) and this is a fundamental observational feature that our models aim to reproduce.

It is particularly important to compare our model results not only with the cloud velocities, but also with the number of observed clouds, because with this additional constraint we have more information on the total mass accretion rate and cloud masses that are needed to successfully reproduce the observations (see Sect. 3), breaking the degeneracy between different cloud infall models.

One limitation of these data is that they only contain observations at relatively small projected radii: we have information only for R <  160 kpc, which is much smaller than the virial radii of these massive galaxies. The average virial radius of these halos is indeed 526 kpc. We discuss the comparison between our results and other observations at larger distances later in Sect. 6.

3. The model

In this Section we describe how we built our dynamical models, starting from the idea that the cool CGM is composed of different clouds (see Sects. 1 and 2), which originate from the cosmological accretion of gas into the galaxy halos and are infalling toward the central galaxies. We first describe, in Sects. 3.1 and 3.2, how we calculated the infall velocities of the clouds for each of the sixteen galaxies of our sample; then in Sect. 3.3 we focus on the total accretion rate of cool gas inside our model halos.

3.1. Purely ballistic model

The main force that drives the motion of the cool clouds is the gravitational force of the dark matter halo. Gravity forces the cool medium to accrete toward the central galaxies. If we considered only this force, the equation of motion would be:

$$\begin{array}{c}\hfill \frac{\mathrm{d}{v}_{\mathrm{fall}}}{\mathrm{d}r}=\frac{1}{v_{\mathrm{fall}}(r)}\frac{GM(r)}{r^2},\end{array}$$(3)

where M(r) is the DM mass within the radius r. Throughout this work, we consider the velocity v with a positive sign if it is pointing toward the center.

We have assumed in our treatment a Navarro Frenk White (NFW) profile (Navarro et al. 1995), described by

$$\begin{array}{c}\hfill M(r)=4\pi {\rho }_0{r}_{\mathrm{s}}^3{\displaystyle [ln(1+r/r_{\mathrm{s}})-\frac{r/r_{\mathrm{s}}}{1+r/r_{\mathrm{s}}}]},\end{array}$$(4)

where

$$\begin{array}{c}\hfill {\rho }_0=\frac{M_{\mathrm{vir}}}{4\pi r_{\mathrm{s}}^3\left[{\displaystyle ln(1+r_{\mathrm{vir}}/r_{\mathrm{s}})-\frac{r_{\mathrm{vir}}/r_{\mathrm{s}}}{1+r_{\mathrm{vir}}/r_{\mathrm{s}}}}\right]}\end{array}$$(5)

is the central density, M_vir and r_vir are the virial masses and radii reported in Table 1, and r_s = r_vir/c is the scale radius, where c is the concentration calculated for each halo following Dutton & Macciò (2014).

By solving Eq. (3) for a characteristic galaxy of our sample (M_vir = 10^13.3 M_⊙ and r_vir = 526 kpc, obtained averaging the values of the 16 galaxies in Table 1), assuming that the clouds start to fall with v = 0 km s⁻¹ from the virial radius, we can obtain the velocity profile shown as a blue line in Fig. 2. Considering only the gravitational force would lead all the clouds to accrete onto the central galaxies with velocities higher than 1000 km s⁻¹ and this prediction is totally inconsistent with the low velocities observed by COS-LRG (see Fig. 1). We conclude that a simple ballistic infall is excluded by the observations and we move to more accurate models. As explained in Sect. 1, the CGM is a multi-phase gas and the cool gas clouds must be surrounded by a hotter medium (corona) which influences their motion. Therefore, to describe the kinematics of the clouds we need, in addition to the gravitational force, to introduce in our model also hydrodynamical effects.

Fig. 2. Velocity profiles of the cool CGM clouds obtained using the ballistic model (Eq. (3), blue line) and the final model that includes all the hydrodynamical interactions (Eq. (12)), with three different choices of mcl, start and α.

3.2. Hydrodynamical effects

To simulate the interactions between the cool clouds and the hot coronal gas, we decided to use an analytical parametrization of the hydrodynamical effects that take place in the CGM. Using an analytic treatment has several advantages, including computational speed, flexibility and physical insight. Our treatment matches however indications from high-resolution hydrodynamical simulations (see below).

First, we defined, for each galaxy of the sample, the corona as a gas at the galaxy virial temperature (2), in hydrostatic equilibrium with the dark matter halo (4), which leads to the electron density profile

$$\begin{array}{c}\hfill n_{\mathrm{e},\mathrm{cor}}(r)=n_{\mathrm{e},0}exp[-\frac{{\mu }_{\mathrm{cor}}{m}_{\mathrm{p}}}{k_{\mathrm{B}}{T}_{\mathrm{vir}}}(\mathrm{\Phi }(r)-{\mathrm{\Phi }}_0)],\end{array}$$(6)

where

$$\begin{array}{c}\hfill \mathrm{\Phi }(r)=4\pi G{\rho }_0{r}_{\mathrm{s}}^2[ln\left(\frac{1+r}{r_{\mathrm{s}}}\right)\frac{r_{\mathrm{s}}}{r}]\end{array}$$(7)

is the NFW potential. Here, ρ₀ and r_s are the same as in Eq. (4), k_B is the Boltzmann constant, m_p is the proton mass, μ_cor = 0.6 is the mean molecular weight, and n_e, 0 is the normalization of the gas profile. We fixed n_e, 0 by requiring that the total coronal mass is equal to 20% of the total baryonic mass within the galaxy halo (assuming M_bar = 0.158M_vir, where 0.158 is the baryon fraction as in Planck Collaboration VI 2018). The 20% value is justified by observational estimates (e.g., Anderson & Bregman 2011; Bogdán et al. 2013), although the exact amount of hot gas present in the galaxy halos is still debated.

We assumed then that the cool CGM clouds are pressure confined by the hot gas, using the following formula

$$\begin{array}{c}\hfill n_{\mathrm{cl}}{T}_{\mathrm{cl}}=n_{\mathrm{cor}}{T}_{\mathrm{vir}},\end{array}$$(8)

where T_cl = 2 × 10⁴ K is the characteristic temperature of the cool absorbers (Zahedy et al. 2019) and n_cor = 2.1n_e, cor is the total particle number density in the corona. Figure 3 shows as an example the density profiles of the hot and cool mediums for an average galaxy with the same properties as the one used in Sect. 3.1. The coronal density is low at large distances from the central galaxy, while it increases in the central regions of the halo. Due to the pressure equilibrium (8), also the densities of the cool clouds, represented by the blue line in Fig. 3, will increase at lower distances from the central galaxy.

Fig. 3. Red: profile of the intrinsic density of the hot coronal gas (from Eq. (6)) for a galaxy with Mvir = 1013.3 M⊙, rvir = 526 kpc, T = Tvir = 6.5 × 106 K (average of the galaxy halos in our sample), and normalization set as explained in Sect. 3.2. In blue, profile of the intrinsic average density expected for the cool circumgalactic medium if in pressure equilibrium with the corona, obtained through Eq. (8).

The first effect of the hot gas on the cloud motion is the deceleration due to the drag force, that we can write (e.g., Marinacci et al. 2011)

$$\begin{array}{c}\hfill {\dot{v}}_{\mathrm{drag}}=\frac{\pi r_{\mathrm{cl}}^2{\rho }_{\mathrm{cor}}{v}^2}{m_{\mathrm{cl}}},\end{array}$$(9)

where ρ_cor = μ_corm_pn_cor, v is the cloud velocity, m_cl is the cloud mass, and r_cl is the cloud radius, given by

$$\begin{array}{c}\hfill r_{\mathrm{cl}}={\left(\frac{3m_{\mathrm{cl}}}{4\pi {\rho }_{\mathrm{cl}}}\right)}^{1/3},\end{array}$$(10)

where ρ_cl = μ_clm_pn_cl is the density of the cool medium and μ_cl = 0.6 is the mean molecular weight for a cool gas in photo-ionization equilibrium. We note that a slightly higher value μ_cl = 0.67 (Sutherland & Dopita 1993) would be appropriate if the cool gas is in collisional ionization equilibrium, which is however uncertain and would have a negligible impact on our results.

Since the density of the cool gas increases with the decrease of the distance from the galaxy, the clouds become smaller during their motion toward the galaxies. Adding the drag force term to Eq. (3) decelerates the clouds during their infall, with a dependence on the cloud mass. Less massive clouds are affected by a stronger drag force and exhibit a higher deceleration. As there are no observational constraints on the mass of the clouds, we let it vary as a free parameter in our model.

With the inclusion of the drag force in the treatment, the results are more accurate and similar to the real behavior of the clouds. However, this simple model still does not take into account the hydrodynamical instabilities that are expected to take place in the system. High-resolution simulations (e.g., Armillotta et al. 2017) show how clouds moving through a hot-gas environment lose mass with time, mostly because of instabilities that strip cool gas from the cloud and mix it with the hotter medium. The result of these interactions is the evaporation of the cool cloud in the hot gas. These mass losses obviously affect the cloud velocities and can significantly change the output of the models. For this reason, we introduced in our modeling a mass loss term, described by

$$\begin{array}{c}\hfill \frac{\mathrm{d}{m}_{\mathrm{cl}}}{\mathrm{d}r}=-\frac{\alpha m_{\mathrm{cl}}}{v}·\end{array}$$(11)

Equation (11) approximates the hydrodynamical instabilities, by making the clouds lose mass during their fall at a certain constant evaporation rate α (see Fraternali & Binney 2008; Marinacci et al. 2010). As for the mass, the exact value of α is unknown and we thus included it in the model as a second free parameter. Since, with the inclusion of Eq. (11), the mass of the clouds is a function of the intrinsic distance, the free parameter that we vary is the mass of the clouds at the starting radius of their infall motion, m_cl, start.

The inclusion of all the effects described in Eqs. (3), (9), and (11) leads to the following final equation of motion:

$$\begin{array}{c}\hfill \frac{\mathrm{d}v}{\mathrm{d}\mathit{r}}=\frac{1}{v(r)}\frac{GM(r)}{r^2}-\frac{\pi r_{\mathrm{cl}}^2(r){\rho }_{\mathrm{cor}}(r)v(r)}{m_{\mathrm{cl}}(r)}·\end{array}$$(12)

In order to solve Eq. (12), the starting velocity v_start of the clouds is needed as an initial condition. We adopted v_start = 0 km s⁻¹, since we have found that to reproduce the observed velocity distribution the clouds must have very low initial velocities (see Appendix B for details). We also fixed the starting radius of the cloud infall motion at the virial radius of the galaxies (Table 1). This choice is justified from the observations of the cool circumgalactic medium extending till distances comparable to the virial radii of these galaxies (Zhu et al. 2014; Huang et al. 2016) and from our assumption that these clouds are inflowing from the external parts of the halos toward the center. We have also explored models with different starting radii for the clouds (Sect. 6).

Solving Eq. (12) we can find the velocities of the clouds with respect to their distance r from the galaxy center, taking into account both the gravitational force and the hydrodynamical effects. The different curves in Fig. 2 show how the velocity profile varies with different choices of the parameters, using the same average galaxy used in Sect. 3.1 for the ballistic motion and the coronal model in Fig. 3. After being initially accelerated by the gravitational force the clouds will be decelerated by the interactions with the corona. Models with lower cloud masses and higher evaporation rates have, particularly in the central regions, clouds with lower velocities, due to the higher efficiency of the hydrodynamical interactions. It is therefore crucial to explore the parameter space to understand which choice of the parameters leads to a velocity profile consistent with what we observe.

3.3. Accretion rate

In the model outlined above the clouds start from the virial radius and then are attracted by the gravitational force toward the center of the halos. The last step in creating the model is to define the total cool gas accretion rate at the the virial radius.

To this purpose, we assumed that the CGM cool clouds are coming from the cosmological gas accretion. Cosmological models in fact predict that the DM halos are growing in time by the accretion of external matter. We can estimate the quantity of matter accreted per unit time by using the prescription in Fakhouri et al. (2010), who proposed a fit for the mass growth rate of the DM halos based on the results of the Millennium simulations (Springel et al. 2005; Boylan-Kolchin et al. 2009). Multiplying this prescription by the baryonic fraction, we obtain the rate of baryonic matter accreting in the halos:

$$\begin{array}{cc}\hfill {\dot{M}}_{\mathrm{cosm}}& =47.6h^{-1}{M}_{\odot }{\mathrm{yr}}^{-1}{\left(\frac{M_{\mathrm{vir}}}{{10}^{12}{h}^{-1}{M}_{\odot }}\right)}^{1.1}(1+1.11z)\hfill \\ \hfill & \times \sqrt{{\mathrm{\Omega }}_{\mathrm{m},0}{(1+z)}^3+{\mathrm{\Omega }}_{\mathrm{\Lambda },0}}\times f_{\mathrm{b}},\hfill \end{array}$$(13)

where h = 0.7, Ω_m, 0 = 0.3, Ω_Λ, 0 = 0.7, and f_b = 0.158 is the baryon fraction by Planck Collaboration VI (2018).

We related the accretion rate to the cosmological accretion, using the formula

$$\begin{array}{c}\hfill {\dot{M}}_{\mathrm{accr}}(r_{\mathrm{vir}})=f_{\mathrm{accr}}{\dot{M}}_{\mathrm{cosm}},\end{array}$$(14)

where f_accr is the fourth and last free parameter that we need in our analysis. If the value of f_accr is equal to 1, all the gas predicted by cosmological models to be accreted into the galaxy halo is what we observe as cool circumgalactic medium.

Once we have defined the value of the accretion rate at the virial radius, we can estimate the mass flux rate as a function of the intrinsic radius, using the formula:

$$\begin{array}{c}\hfill {\dot{M}}_{\mathrm{accr}}(r)=4\pi {\rho }_{\mathrm{cool}}(r)v(r)r^2={\dot{M}}_{\mathrm{accr}}(r_{\mathrm{vir}})f_{\mathrm{mass}}(r),\end{array}$$(15)

where ρ_cool(r) is the volume averaged mass density of the cool gas and f_mass = m_cl(r)/m_cl, start is a term that takes into account the mass losses of the clouds at every radius. In this way, the mass flux decreases at smaller distances from the center, consistently with the evaporation of the clouds in the hot gas. For models with α = 0, which means no evaporation of the clouds, the mass flux is constant with radius.

Inverting Eq. (15) we obtained the value of ρ_cool(r), which will be useful in the next Section in order to compare the results with the observations.

4. Comparison with the observations

In this section we explain how we obtained, starting from the intrinsic quantities calculated in Sect. 3, the results that are directly comparable to the observations. To this purpose, we created random distributions of cool clouds over the whole halo, each of them with the properties defined by the model, and we observed them through synthetic observations.

The quantity ρ_cool(r) describes the total mass density of the cool gas per unit of volume. If we divide this density by the mass of the clouds, which is a function of the intrinsic radius, described by Eq. (8), we obtain the number of clouds per unit of volume. The integral of this quantity over the whole volume of the halo is equal to the total number of cool absorbers expected in the system:

$$\begin{array}{c}\hfill N_{\mathrm{cl}}=4\pi {\displaystyle {\int }_0^{r_{\mathrm{vir}}}}\frac{{\rho }_{\mathrm{tot}}(r)r^2}{m_{\mathrm{cl}}(r)}\mathrm{d}r.\end{array}$$(16)

We created then a distribution with a total number of objects given by Eq. (16). The clouds are not distributed uniformly in the galaxy halo, as their probability to be at a certain radius r is predicted by the argument of the integral in Eq. (16). Therefore, we populated the halos using a probability density function (PDF) given by:

$$\begin{array}{c}\hfill \frac{4\pi {\rho }_{\mathrm{tot}}(r)r^2}{m_{\mathrm{cl}}(r)N_{\mathrm{cl}}}·\end{array}$$(17)

We have then N_cl objects, each of them located at a certain radius r and with all the properties given by the results of Eqs. (8), (10)–(12). We associated then to every cloud two other coordinates θ and ϕ, taken in order to have random distribution of clouds over a sphere (see Appendix A for the detailed 3D geometry).

This setup was created for each one of the galaxies reported in Table 1. Figure 4 is a cartoon representation of the cloud distribution, flattened in two dimensions: the clouds are randomly distributed throughout the halo, each one of them with the properties explained above. The last step was to create synthetic observations: we selected the observed clouds tracing a line-of-sight through the galaxy halos, as explained in detail in Appendix A and shown in Fig. 4. We note that at the same projected distance we can intercept clouds at very different intrinsic distances from the central galaxy, potentially till the virial radius of the halo.

Fig. 4. Cartoon 2D representation of the synthetic observations performed in our model. The clouds (blue circles) are distributed throughout the halo (the sphere with radius rvir), all accreting towards the central galaxy (the central red star) with a radial velocity described by Eq. (12). The vertical gray line represents the line of sight, consistent with the projected distance R of the galaxy from the background quasar, while the area included between the two vertical dashed lines represents the range spanned by the COS-LRG observations. The observer is located at the bottom side of the panel.

After selecting the clouds, we calculated their velocity projections along the line of sight (Appendix A). Figure 4 shows the range of the projected distances between the galaxies and the quasars in the COS-LRG observations. As already mentioned in Sect. 2, these impact parameters are all much smaller than the typical virial radii of our massive galaxies. Therefore, the observations span only a small central area of the halos. As a result, most of the clouds that we observe have a direction of the motion that is almost parallel to the line-of-sight and the projection corrections of the velocities are relatively small. We repeated this procedure for all the galaxies and finally we put together all the observed velocities in a single line-of-sight velocity distribution.

The modeling described above allowed us to find the choice of parameters for our models that best reproduces the observations. To this purpose, we performed a Markov chain Monte Carlo (MCMC) analysis over the parameter space defined by m_cl, start, α, and f_accr. We now describe the likelihood of the models, needed to perform the MCMC. This is obtained using the comparisons between our models results and the two observables outlined in Sect. 2: the line of sight velocity distribution and the number of clouds. The total likelihood is the product of the likelihoods associated to these two observational constraints. In particular, we used the logarithm of the likelihood

$$\begin{array}{c}\hfill ln{\mathcal{L}}_{\mathrm{tot}}=ln{\mathcal{L}}_1+ln{\mathcal{L}}_2,\end{array}$$(18)

where ℒ₁ and ℒ₂ are the two different likelihoods.

To compare the two velocity distributions we have to assume for both of them the same number of velocity bins n_bin in which the velocity range is divided. We considered a velocity range spanning from −1000 to 1000 km s⁻¹, to include all the possible velocities predicted by our models, and we divided it in 20 different bins. However, we found that slightly different choices of the number of bins do not affect the final results. The comparison between the two velocity distributions is then given by the reduced chi squared:

$$\begin{array}{c}\hfill ln{\mathcal{L}}_1=-\frac{1}{2(n_{\mathrm{bin}}-3)}{\displaystyle \sum _{b_i}}\frac{{|d_i-m_i|}^2}{{\sigma }_i^2},\end{array}$$(19)

where n_bin = 20, n_bin − 3 is the number of degrees of freedom, d_i and m_i are respectively the values of the observed and model distributions in the velocity bins b_i, and σ_i are the errors on the observed values, calculated using the bootstrapping method¹.

Finally, the second constraint is given by the comparison between the number of clouds observed in our synthetic models (n_mod) and the total number of clouds observed by COS-LRG (n_obs = 42, see Sect. 2), using

$$\begin{array}{c}\hfill ln{\mathcal{L}}_2=-\frac{1}{2}{{\displaystyle \left(\frac{|n_{\mathrm{obs}}-n_{\mathrm{mod}}|}{\sqrt{n_{\mathrm{obs}}}}\right)}}^2,\end{array}$$(20)

where $\sqrt{n_{\mathrm{obs}}}$ is the standard deviation of the Poisson distribution with a mean value equal to n_obs. The priors used in this analysis are flat for all the parameters but f_accr, which sets the value of the total mass accretion rate of cool gas at the virial radius. In particular we employed a gaussian prior centered on 1 and with a dispersion of 0.5, in order to have the accretion rate consistent with the total baryonic accretion rate predicted by cosmology (Eq. (13)). Negative values of f_accr are excluded. Regarding the other two parameters, the logarithm of the cloud mass (in solar masses) is allowed to vary uniformly in the range from 4 to 7, while the evaporation rate α can vary uniformly from 0 to 3.5 Gyr⁻¹. A negative α would in fact mean that the clouds are gaining cool mass and we excluded this scenario, because condensation is highly unlikely at these large temperatures: the cooling time of the hot gas is indeed too high for it to cool down and increase the amount of cold gas in the system. This is also supported by hydrodynamical simulations (Armillotta et al. 2016).

5. Results

In this Section we report the results of the MCMC analysis that we used to compare our results to the observations and to find the best models that describe the CGM of massive ETGs. Then we discuss the physical scenario arising from our models.

5.1. MCMC results

The results of the MCMC study are reported in Fig. 5 and Table 2. Figure 5 shows the corner plot where both the one and two dimensional projections of the posterior probabilities of the three parameters are shown, while in Table 2 we show the 32nd, 50th (the median value) and 68th percentiles of the parameter one dimensional probabilities. Models with parameters in the ranges reported in Table 2 are the best in maximizing the total likelihood ℒ_tot and therefore the best in reproducing the COS-LRG observations.

Fig. 5. Corner plot with the MCMC results, representing the one and two dimensional projections of the posterior probabilities for the three free parameters of our models.

In the left and right panels of Fig. 6 we show the outputs of these models, respectively the line of sight velocity distribution and the total hydrogen column densities, with their comparison with the observations. The errors in the observed distribution are calculated using the bootstrapping method explained in Sect. 4. The model velocity distribution is dependent on the random populations of clouds, as explained in Sect. 4. Therefore, to have a more robust comparison between our model output and the observations, the model distribution in the left panel of Fig. 6 is the result of the average of 2000 models with the values of the parameters included within the 32nd and the 68th percentiles reported in Table 2, weighted with the value of the posterior probability of each model. The model and observed velocity distributions are both normalized to have a subtended area equal to 1. The kinematics of the model clouds, fitted by the MCMC analysis, is in good agreement with the observed one. The consistency between the two distributions is also evident from the comparison between the cumulative distribution functions (CDFs) in Fig. 9. The result of a Kolmogorov-Smirnov test between the two samples confirms that they are likely drawn from the same distribution. The averaged number of clouds of the 2000 randomizations used to create the plot in the left panel of Fig. 6 is 43 ± 10 (the error is given by the weighted standard deviation), consistent with the observed one of 42.

Fig. 6. Comparison between the observations (orange bars and circles) and our best model results (blue lines). Left panel: normalized line-of-sight velocity distribution, where the errors in the observations are obtained with the bootstrapping method. Right panel: total hydrogen column density. The downward arrows represent the upper limits of the non-detections.

The hydrogen column density of our model (blue curve) in the right panel of Fig. 6 is instead obtained converting the density of clouds ρ_tot (see Sect. 3.3) into a numeric hydrogen density and projecting it along the line of sight. More precisely, the model curve in Fig. 6 is the average of the sixteen column densities calculated for each galaxy. This profile is not dependent on the random populations of clouds and therefore we show here the result of a single model with the parameters given by the median values reported in Table 2. The orange circles represent instead the column density observations by Zahedy et al. (2019), also reported in Table 1. Although they are the most accurate determination of these values, they are still affected by various uncertainties, depending on the different assumptions made to perform the photo-ionization models. These uncertainties can potentially be very large and are not fully quantified by the error bars shown in Fig. 6. As an example, the COS-Halos collaboration (Werk et al. 2013, 2014; Tumlinson et al. 2013) analysed with similar photo-ionization models and slightly different assumptions (see Prochaska et al. 2017) five of the galaxies of our sample, finding column densities that can be one or two orders of magnitudes different from the COS-LRG estimates. Because of these uncertainties in the data, we have not included the column densities in the MCMC fit. We can note however that our predictions are not inconsistent with the observations. Interestingly, with our models we can also make predictions about the density profile in the outer part of the halos: we will discuss the implications in Sect. 6.3.

We conclude that, with the correct choice of parameters, models of clouds that are starting from the virial radii of the halos, infalling at a cosmologically motivated rate and then interacting with the hot coronal gas, can successfully reproduce the observed kinematics, number of absorbers and column densities of the cool CGM around early-type galaxies. This is our first and main result.

In Sect. 5.2 we analyze the motion and the properties of our model clouds and we discuss the physical scenario that these particular models describe.

5.2. Description of the cloud motion

Using the median parameters of Table 2 we can have an indication on the properties of the best models that reproduce the observations and we can analyze in detail the motion of the infalling clouds. In Figs. 7 and 8 we report our results on the cool CGM properties, all with respect to the intrinsic distance r from the central galaxies. For the moment we do not discuss the possible origin of the cool clouds (Sect. 6.5), but we just describe their behavior throughout their infall. The different blue-tone lines in Figs. 7 and 8 represent the results obtained for the 16 different galaxies that we are modeling. Since we are describing all the galaxies with the same model, the general trends of the cool CGM properties are similar in all the objects and the arising physical picture is the same for all the galaxies of our sample, with small variations due to the different virial masses and radii.

Fig. 7. Properties of the cool CGM clouds described by our best model. All the profiles are plotted with respect to the intrinsic distance from the central galaxies and the different colored lines represent the result for the sixteen different galaxies of the sample. Top left: velocity profiles of the cool clouds. Top right and bottom left: respectively masses and radii of the cool clouds. Bottom right: cool CGM cumulative mass profile, with the vertical dashed line representing the distance of 160 kpc.

Fig. 8. Same as Fig. 7, but for the total mass flux profile of the cool CGM.

In the top left panel of Fig. 7 we show the cloud velocity profiles. The clouds start at the virial radii with very low velocities, they are accelerated by the gravitational force and then, as already explained in Sect. 3, they are decelerated by the drag force and the other hydrodynamical interactions. We stress that it is this hydrodynamical deceleration that allows the resultant line-of-sight velocity distribution to be consistent with the observed one (see left panel of Fig. 6).

With our analysis we can also estimate, only using our analytical description of these systems, the mass and size of the cool CGM clouds, which were so far largely uncertain (Stocke et al. 2013; Werk et al. 2014). The top right and bottom left panels of Fig. 7 show the evolution of these two quantities and the direct effects of the evaporation and the pressure equilibrium on the clouds. The profiles in the top panel show that the clouds start at the virial radii of the galaxies with a mass of ∼10⁵ M_⊙, then they progressively lose mass due to the interaction with the corona. In the central region they have lost more than 99% of their mass and have masses of about 10² M_⊙. Also the radii of the clouds (bottom panel), obtained from Eq. (10), strongly decrease with time, starting from r_cl ∼ 1 kpc at the virial radius and reaching sizes of a few pc in the central regions, due to the effects of both the cloud evaporation and the pressure equilibrium (in the central regions the corona and cloud densities are higher). The conclusion is that in the central regions the cool CGM clouds are much smaller and also less massive, because most of their mass is evaporated in the hot corona.

From the bottom right panel of Fig. 7, which shows the cumulative mass profiles of the cool gas, we can infer how the mass loss of the clouds influences the distribution of the cool CGM mass throughout the halo. Most of the cool gas mass is concentrated in the external regions of the halos, because of the volume effect (their volume is much larger than the one of the internal parts), but also because the clouds are more massive. In fact, because of the cloud evaporation, the decline of the mass profiles is much steeper than what expected without evaporation and the inner regions are almost devoid of cool gas. While on average the total mass of the cool medium within the galaxy halo is 2.5 × 10¹¹ M_⊙, the cool gas mass within r = 160 kpc (highlighted by the vertical dashed line) is only ∼1% of the total. Therefore, we are describing with our models a system where most of the cool CGM gas is not reaching the central regions and is not accreting onto the galaxy. This is also shown by the cool gas accretion rate, reported in Fig. 8: starting from a cosmological rate at the virial radius, the accretion in the inner parts is reduced to less than 1% of the initial value.

With this result we are therefore in agreement with the observations of a great amount of cool CGM in these halos and we also go in the direction of explaining why this large reservoir of cool gas is not feeding the star formation of these thus passive central early-type galaxies. We discuss this further in Sect. 6.1.

6. Discussion

With our analysis we have found that infall models of cool clouds can successfully describe the COS-LRG observations of the cool CGM around massive ETGs. We presented our findings in Sect. 5. In this section we discuss all the possible limitations of our analytical models (Sects. 6.1 and 6.4), the comparison with other observations (Sect. 6.3) and the implications of our findings for the cool gas state, origin and fate (Sects. 6.2 and 6.5).

6.1. Cloud survival

We have seen in Sect. 5.2 that with our models we describe a scenario where, due to the evaporation of the CGM clouds, only a very small part of the cool gas accreting into the galaxy halos is reaching the inner regions. Nevertheless, we also find that this small amount of cool gas (on average a few M_⊙ yr⁻¹) is reaching the central ETGs, possibly feeding the galaxy star formation at a similar rate. In this Section we discuss whether this problem may be alleviated by additional effects that we have not taken into account, focusing in particular on the survival time of the cool clouds.

As already explained, the analytical model of Sect. 3 is only an approximation of all the hydrodynamical interactions taking place in the CGM. As a consequence, the clouds described in our treatment can lose mass, but they cannot completely evaporate in the hot gas or be destroyed by the cloud-corona interactions, which is instead what we expect in reality and what we observe in high-resolution simulations (e.g., Heitsch & Putman 2009). Therefore, it is important to estimate the time after which the clouds are probably completely evaporated in the surrounding hot medium. We consider as an estimator of the cloud survival time the conduction time, which is the time needed for the thermal conduction to be efficient and therefore for the cloud to be entirely evaporated in the corona. In our massive galaxy halos, where there is a high temperature gradient between the different gas layers (the coronal gas is almost three orders of magnitudes hotter than the cool clouds) we expect the thermal conduction to be particularly efficient and to play a major role in the evaporation of the clouds (Armillotta et al. 2017).

The conduction time is defined as (Spitzer 1962)

$$\begin{array}{c}\hfill t_{\mathrm{cond}}=\frac{n_{\mathrm{cool}}{k}_{\mathrm{B}}{r}_{\mathrm{cl}}^2}{f{k}_{\mathrm{sp}}},\end{array}$$(21)

where f is the suppression factor due to the magnetic field (e.g., Chandran & Cowley 1998) and:

$$\begin{array}{c}\hfill k_{\mathrm{sp}}=\frac{1.84\times {10}^{-5}{T}^{5/2}}{ln\mathrm{\Psi }}\mathrm{erg}{\mathrm{s}}^{-1}{\mathrm{K}}^{-1}{\mathrm{cm}}^{-1}\end{array}$$(22)

is the heat conduction coefficient, where lnΨ is the Coulomb logarithm (≈30) and T is the temperature in K of the hot gas. We estimated the mean value of this time for a cloud at the virial radius with m_cl, start ∼ 10⁵ M_⊙ (see Table 2) and f = 0.1, obtaining t_cond ∼ 200 Myr. This has to be compared to the mean infall time of a cloud in our best models, calculated using the formula

$$\begin{array}{c}\hfill t_{\mathrm{fall}}={\displaystyle {\int }_{r_{\mathrm{vir}}}^0}\frac{\mathrm{d}r}{v(r)}·\end{array}$$(23)

Averaging this time over the seven galaxies of our sample, we obtained that t_fall ∼ 3.7 Gyr, much longer than the estimated conduction time. This means that the clouds cannot survive their entire journey from the virial radius to the central galaxy and they most likely completely evaporate in the hot gas at relatively large distances from the center of the halos. Therefore, the cool CGM is not feeding the central galaxies, thus explaining the quiescence of these objects. However, the comparison between these two times is not straightforward, since there are uncertainties on both quantities. The derived value of the conduction time is in fact dependent on the magnetic field, whose intensity and orientation is unknown. If the magnetic field suppression is larger than what we have considered in our estimate, the suppression factor f in Eq. (21) would be lower and therefore the conduction time would be longer than what we have found. On the other hand, since the clouds become smaller with time, the conduction is more efficient in the inner regions of the halos and therefore our estimate is an upper limit on the survival time of the clouds. Finally, there are also uncertainties on the infall time obtained through Eq. (23), since we are considering in our analysis only a radial infall. Adding any tangential component would probably increase the final estimate of the total cloud infall time. Based on these considerations, we consider very likely that the clouds cannot survive for the entire infall time, although it is hard to define precisely the region of the halo where they completely evaporate.

The kinematic distribution of our model is arising from clouds distributed all over the halos, while we have argued with the last estimation that the central regions are probably devoid of cool gas. To prove that this property does not affect our final findings, we have performed again the MCMC analysis, with the same likelihood expressed in Sect. 4 and the same parameter space outlined in Sect. 5.1, but with the central regions of the halos forced to be completely devoid of clouds. To do this, we have defined the destruction radius r_destr as the inferior limit of the integral in Eq. (16) in place of 0 and as the smallest intrinsic radius where the clouds are present in our distribution. In this way we can assess how much the cloud destruction influences the final results. We have run a second MCMC, with r_destr = r_vir/2, obtaining the results reported in Table 3. We can see how the best-fit parameters are consistent within the errors with the ones found with the first MCMC, reported in Table 2. The physical scenario is the same described in Sect. 5.2, but with the cool clouds by construction not reaching the central galaxies. Also with this kind of model we are able to successfully reproduce the two observational constraints. In Fig. 9 we show the comparison between the cumulative distribution functions of the velocity distributions of the models and the observations: the kinematic distribution of the last model is perfectly consistent both with the result of our main model (Sect. 5) and with the observations, with a number of clouds equal to 44 ± 11. The reason why the results of our model did not change with this last modification is that, even without imposing the cloud destruction, the vast majority the clouds observed by our synthetic observation analysis (see Sect. 4) are residing in the halo external regions. The radii of the external absorbers are in fact up to two orders of magnitudes larger than the internal ones (bottom left panel of Fig. 7), implying much higher cross sections and therefore higher chances to be observed.

Fig. 9. Cumulative distribution functions of the line-of-sight velocity distributions of the observations (orange curve) and of our models. In blue we report the result of our main model (Sect. 5), while the dashed curves represent models with rstart = rvir/2 (green, dashed) and rdestr = rvir/2 (red, small dashed).

We can finally conclude that our model results, although we cannot directly simulate the complete evaporation of the cool clouds, provide a coherent picture where the pressure-confined cool CGM absorbers are: 1) infalling toward the center, 2) slowed down by the coronal drag force and 3) completely destroyed by the interactions with the hot gas before reaching the galaxy, thus not feeding any central star formation.

A similar scenario was discussed by Huang et al. (2016) to explain the kinematics of cool CGM clouds around massive elliptical galaxies: they argue that the observed low velocities could be reproduced by clouds with m_cl ∼ 5 × 10⁴ M_⊙, slowed down by the effect of the hot gas drag force, though they do not assess this possibility with detailed dynamical models, as we have developed with our work. Also Chen et al. (2018) discuss the same physical scenario for the cool CGM, arguing that most of the clouds with m_cl <  10⁶ M_⊙ cannot reach the center if they originate at large distances from the galaxies, in agreement with our final findings. Our work is the first attempt to develop a coherent dynamical model to describe the CGM around massive ETGs. Similar semi-analytical models for the cool circumgalactic medium have recently been developed by Lan & Mo (2019), though they focus on 10¹² M_⊙ halos around star-forming galaxies, for which they find that an outflow model can qualitatively reproduce several CGM observed properties. Our approaches are similar, but we concentrate in this work on more massive, passive galaxies and we directly compare our results with the full observed velocity distribution of the cool gas clouds.

6.2. Pressure equilibrium

A fundamental assumption that we made for the cool CGM physical state is the pressure equilibrium between the cool absorbers and the coronal gas. This is physically justified, since we expect any pressure imbalance with the ambient medium to be erased on a cloud sound crossing time-scale, which is roughly ∼10 Myr and therefore much smaller than the dynamical time.

It is hard to compare our intrinsic densities with the COS-LRG data, because of the already mentioned uncertainties in the observations, due to photo-ionization modelling and, most importantly, to projection effects. However, Zahedy et al. (2019) find, through their photo-ionization analysis, that the cool gas is likely in pressure equilibrium with the hot halo, in agreement with our assumption of clouds pressure-confined by the corona. We conclude that, in addition to being physically motivated, this assumption is also supported by observations (although earlier studies had been suggesting different results, see Werk et al. 2014). We plan to investigate these aspects further, including explicit photo-ionization modelling, in future works.

6.3. Comparison with other observations

The results of this work are obtained from the comparison of our model predictions with the observations of the COS-LRG collaboration, which analyzed both metal and hydrogen CGM absorption lines with a very high spectral resolution, allowing us to have detailed data on the kinematics of the cool clouds around 16 massive ETGs (Sect. 2). A comparable analysis has been done by the COS-Halos program (Werk et al. 2013, 2014; Tumlinson et al. 2013; Prochaska et al. 2017), a larger survey of the cool CGM around galaxies with a broad stellar mass range. Five of the galaxies analysed by this program are also part of the COS-LRG sample. We have applied our models also on this sub-sample of five objects, using only the kinematic data of the COS-Halos collaboration. The results of the MCMC analysis are totally consistent with the ones showed in this paper.

However, these data sets have also some limitations, mostly due to the lack of observations at projected radii larger than 160 kpc (in both samples the projected distance of the galaxy from the background quasar is always smaller than 160 kpc). It is important therefore to compare our findings with observations at larger projected distances. Zhu et al. (2014) and Huang et al. (2016) have observed MgII absorptions from cool CGM around large samples of ETGs out to very large distances from the central galaxies and they both found that the cool gas has a kinematics consistent with what is found by COS-LRG and COS-Halos. Their data have not been included in our analysis, because the resolution of these observations is not sufficient to disentangle the different clouds along a single line of sight. Also, they observe only the MgII, which could be not representative of the behavior of the whole gas. However, it is interesting to compare our predictions with their results.

As seen in Sect. 5.1, with our model we can make predictions on the cool CGM behavior at large distances, in particular for the projected profile of the total hydrogen column density. The result of our best model is reported as a blue curve in the right panel of Fig. 6: the column density is only slightly declining at large distances from the galaxies. This behavior is expected because our model predicts that most of the cool gas is concentrated in the halo external regions, es explained in Sect. 5.2. Therefore, the amount of gas observed by a line of sight at a large impact parameter is almost equal to the one observed at low projected distances from the galaxy and the resulting projected profile is nearly flat. By construction, the column density of each galaxy goes to zero at the virial radius, beyond which, in our modeling, there is no cool gas. The gradual decline of the column density profile is therefore also due to the fact that different galaxies live in different dark matter halos, with the small bumps corresponding to the different virial radii (reported in Table 1), where the column density of each galaxy goes to zero.

Huang et al. (2016) show that there are strong MgII absorptions up to 500 kpc from the central galaxies, with the covering fraction declining at projected radii larger than 100 kpc. Zhu et al. (2014), who analyze stacked spectra, observe MgII absorptions till and beyond the galaxy virial radii, finding that the total MgII surface density is decreasing with the projected distance. Both works are therefore consistent with our findings of cool absorbers still present at very large distances from the central galaxies and the decrease in the absorption strength, which is less than two order of magnitudes from the center to the virial radius, is roughly compatible with the decline of our total hydrogen column density.

Recently Berg et al. (2018) have analyzed the HI in the CGM around 21 massive ETGs, with sightlines within ∼500 kpc from the central object. They find a covering fraction of the cool gas that is decreasing going from the central regions to high impact parameters. Our models also predict a decline in the projected profile of the covering fraction that is consistent with these observations, although it is hard to make a quantitative comparison between observations and model predictions. Since Berg et al. (2018) have not provided kinematic estimations, we have not included these data in our analysis.

There is also the possibility that the cool clouds start from distances smaller than the virial radius (see Sect. 6.5.2) and in this scenario the steepness of the column density and the covering fraction profiles would be different, with very low values at large distances. With our current dataset we cannot directly constraint the cloud starting radius, since we are using observations till 160 kpc. Therefore, instead of considering the starting radius as an additional free parameter of our models, we have run an MCMC analysis as in Sect. 4, but fixing the starting radius at the value of r_start = r_vir/2. We present the results of this fit in Table 3: the best parameters are consistent with our previous findings, with a physical scenario similar to the one described in Sect. 5.2. Also, the line-of-sight velocity distribution and the number of clouds (49 ± 13) resulting from this model are consistent with the observations, as shown by the CDFs comparison in Fig. 9. Therefore, even models with smaller starting radii can successfully describe the cool CGM properties. Future high-resolution observations of cool CGM absorbers at impact parameters comparable to the virial radius will help distinguishing between different scenarios.

6.4. Uncertainties in the accretion rate

In Sect. 5.1 we have shown that the best models in reproducing the observations have a mass accretion rate which is slightly (1.2 times) higher than the prediction of cosmological models, while in principle we would not expect the cool mass accretion to be greater than the cosmological rate. The reasons for this unexpected value are probably due to the simplifications of our models. The first reason could be the simplified structure of our model clouds. They are in fact spherical, with their radius becoming smaller as they fall throughout the halos. Although this is probably a good approximation for the head of the clouds, hydrodynamical simulations show that the interactions with the hot gas create a tail of partially mixed gas which can be much more extended than the main head. We cannot model this effect with analytical descriptions, but we argue that it would mainly have an effect on the number of clouds observed by the synthetic observations of our models. The cross section of the clouds will indeed be larger and we would synthetically observe more clouds than what we are doing with our current models. Therefore, to observe the same number of clouds we would need a lower total accretion of cool gas, in agreement with the prediction of the cosmological models. We note that this effect would likely have also an impact on the shape of the model velocity distribution. With our simplification of the structure of the clouds, we are indeed underestimating mostly the cross section of the slowest clouds, that had more interactions with the corona. Therefore we expect, from a more thorough analysis, to observe more clouds with very low velocities, populating the very central region of the model velocity distribution (see the left panel of Fig. 6) and alleviating the tension in that bin between the observations and our predictions.

Another reason why we are probably overestimating the accretion rate is instead given by our assumption of isotropic inflow onto the galaxy halos. As explained in Sects. 3.3 and 4, to calculate the mass accretion rate of the cool gas we assumed that the cool gas is inflowing isotropically, starting from the virial radius. However, the accretion could be anisotropic, with the cool medium inflowing through filaments toward the galaxies (cool gas filamentary accretion is common in high-resolution cosmological simulations, see Nelson et al. 2016). In this scenario, the filaments would not fill the entire volume of the halos and the total accretion of gas would be different than what predicted by our model. However, the low number statistics does not allow us to constraint the real geometry of the gas inflow and to understand whether and how much we are overestimating the cool CGM accretion with our simplified, isotropic models. We leave a more detailed analysis about the geometry of the gas accretion for future works.

6.5. Origin and fate of the cool absorbers

In the previous Sections we have focused on the physical description of the cool absorbers behavior throughout the halos. In this Section we instead analyze the possible origin scenarios for the cool CGM, as well as their consequences on the galaxy evolution. The aim of this Section is only to speculate on the different possible pictures that can arise from our models, leaving a detailed analysis of the cool CGM origin and fate for further works.

6.5.1. External origin

Our favourite interpretation for the origin of the cool CGM is that it is coming from the accretion of external gas. As explained in Sect. 3 in fact, cosmological models predict that galaxy halos acquire gas at a rate described by Eq. (13) and our assumption is that the cool CGM is originated by this accretion of external gas: in our models we set indeed the cool gas accretion rate in order to be consistent with the cosmological one. In this picture, the filaments of low-metallicity intergalactic gas accrete into the galaxy halo and form the cool clouds that fall down through the halos as described by our modeling.

This assumption is motivated by the COS-LRG analysis of the cool CGM metallicities (Zahedy et al. 2019): they find in fact that about half of the galaxies in their sample host metal poor clouds, consistently with our picture of accretion of external gas, which is expected to have low metallicities (see Lehner et al. 2016). Also a recent work by Chen et al. (2019) states that very low metallicity gas is present in the halos of massive ETGs, probably coming from the accretion of intergalactic medium. The scenario of external origin has more difficulties to explain the observations of absorbers with higher metallicities, that are common in the COS-LRG survey and also in previous surveys like COS-Halos. Prochaska et al. (2017) find indeed that the CGM absorbers in the whole COS-Halos sample span a broad range in metallicities, with a median value around 30% solar. A way to produce these metal-rich absorbers is to invoke internal origins for the clouds, which we will discuss in Sect. 6.5.2.

Another way to explain the observation of a broad range of metallicities is with the destruction of the cool accreting filaments. It is possible in fact that the clouds are formed by the fragmentation of the filaments and their mixing with the corona, explaining higher metallicities than what expected in the IGM. Also, in this scenario, the filaments are slowed down by these processes and it is reasonable to believe that the resulting clouds will start with low velocities, as indicated by our results² (see Appendix B). Due to these considerations, we consider the external accretion consistent with our results and the most likely picture for the cool CGM absorbers around ETGs.

6.5.2. Alternative origins

A plausible alternative origin for the cool curcumgalactic gas in the halos of massive ETGs would be from the stripping of ISM of satellite galaxies by tidal or ram pressure forces. Signatures of cool gas coming from this type of interactions have been observed directly in emission in different massive halos (e.g., Epinat et al. 2018; Johnson et al. 2018). The cool gas stripped by the satellite could eventually fall down towards the central massive galaxy, in a journey compatible with the motion described by our model. Although the accretion from external gas remains our favourite scenario, we do not reject this last picture as a possibility for the cool CGM origin.

Another scenario that can give rise to cool clouds falling through the halos is related to the presence of outflows from the central galaxy, which we have not considered so far in our treatment. This neglect is justified by the nature of the galaxies that we are studying, that are by definition in a quiescent state and should therefore not experience strong outflows of gas. However, although they are currently quiescent, there could have been outflows in the past of these objects, due to either star formation bursts or AGN activity. AGN are indeed known to drive massive multi-phase outflows (e.g., Greene et al. 2012) and cool, high velocity gas is often observed in the CGM of quasars (e.g., Bowen et al. 2006; Prochaska et al. 2013; Johnson et al. 2015). It is possible to think that, as in the recycling scenario (see Sect. 1), the outflows due to recent star formation or AGN episodes have reached the external parts of the halos, cooled and eventually fallen down again toward the center. In this picture we are observing the last part of this cycle, composed by the infalling cool CGM clouds, which in this case have an internal origin. This would mean that we are not observing a continuous accretion but a sporadic event of infall due to a recent outflow episode happened in the central galaxies. Although this picture would explain the high metallicities observed in part of the cool CGM (Prochaska et al. 2017), all our galaxies should have experienced an outflow burst at similar times, with comparable timescales and we consider this scenario unlikely. Thus, the recycling is most probably not the main driver of the cool CGM formation around early-type galaxies.

The last possibility is that the clouds are condensing out of the hot gas due to thermal instabilities, developing a multi-phase halo. Multi-phase gas has indeed been observed in both nearby clusters and elliptical galaxies (Voit & Donahue 2015; Voit et al. 2015) and the same origin scenario has been proposed also for the cool CGM around massive ETGs by different observational work (see Huang et al. 2016; Zahedy et al. 2019). This picture is motivated by various theoretical works (e.g., Sharma et al. 2012; McCourt et al. 2012; Voit 2018), which predict that the instabilities can develop at a cooling radius r_c <  r_vir, where the cooling time is comparable to the dynamical time of the gas. In this scenario, clouds are created within the galaxy halos with very low velocities, (which may be of the order of the turbulence of the hot coronal gas, see Voit 2018) and then start to fall down toward the central galaxies, consistently with our main findings. However, it is so far not clear if such instabilities can spontaneously develop in the hot corona (e.g., Binney et al. 2009; Nipoti & Posti 2014) and therefore it is difficult to establish whether this scenario could be an important contributor to the cool CGM.

We conclude then that the bulk of the cool circumgalactic gas is likely coming from the accretion of the external IGM, although the different internal processes outlined above may originate part of the observed cool clouds, also explaining the broad range in the observed metallicities.

6.5.3. Fate of the cool gas

Many works in the last decades have studied the fate of the accretion of cool gas from the IGM onto galaxy halos, both with analytical prescriptions and cosmological simulations (e.g., Birnboim & Dekel 2003; Ocvirk et al. 2008; Dekel et al. 2009; Nelson et al. 2013, 2016). The general picture is that, in low mass DM halos, filaments can penetrate through the halo reaching the central galaxy, while in the most massive objects the filaments are shock heated to the galaxy virial temperature and are unable to feed the central galaxy. In this scenario, the CGM clouds that we are studying are the products of the fragmentation of the cool gas filaments. In agreement with the general picture outlined above, we predict that these cool gas components are not able to survive their journey and reach the central galaxy (see Sect. 6.1), but they instead evaporate in the hot corona at large distances from the centre. We therefore believe that these absorbers cannot be straightforwardly linked with the observations in less massive halos of High Velocity Clouds (HVCs, see Putman et al. 2012), clouds of cold gas observed mostly in HI both around the Milky Way and nearby star-forming galaxies. HVCs reside in fact primarily within 10 kpc from the central galaxy disk and their origin is still debated (Fraternali et al. 2015). On the contrary, we expect our CGM clouds around ETGs to be completely destroyed by the interactions with the hot gas at much larger distances. The fate of the cool CGM, by the evaporation of the clouds, is therefore to increase the mass of the hot corona.

From Fig. 8 we can infer the rate of cool gas accreting onto the galaxy halos: on average, the gas accretion at the virial radius is about 500 M_⊙ yr⁻¹. If it is coming from the cosmological inflow, this amount of cool gas is continuously injecting baryons in the galaxy halo, by evaporating into the hot coronal gas. At this rate, (which could be slightly lower, due to the considerations of Sect. 6.4) the mass of the corona would be doubled in about 2 Gyr. This may be related to the current and future location of the so-called missing baryons (Bregman 2007; McGaugh 2008), although a more thorough discussion of these aspects is outside the scope of this work.

7. Summary and conclusions

In this work we carried out a detailed dynamical analysis of the cool circumgalactic medium around massive early-type galaxies at low-redshift, in order to explain the presence of a substantial amount of cool gas around quiescent galaxies and the narrowness of the observed velocity distribution of the cool CGM clouds.

In particular we focused on the observations of the COS-LRG collaboration around 16 massive ETGs. We developed semi-analytical models of cool clouds infalling from the external regions of the halos to the central galaxies, approximating the hydrodynamical interactions of the clouds with the hot corona with analytical parametrizations. These models aim to reproduce the two observational constraints given by the kinematic distribution and the number of the cool clouds. Our model has three free parameters that account for the properties of the infalling clouds and for their interaction with the ambient medium. We constrained these free parameters on the available data with an MCMC analysis.

From our results we draw the following conclusions:

1. models of cool CGM clouds infalling through the galaxy halos at a cosmologically motivated rate can successfully reproduce the observations, explaining the number of observed clouds and their line-of-sight velocity distribution. These results are also consistent with the total hydrogen cool CGM column densities reported by COS-LRG;

2. our best models describe clouds with an initial mass m_cl ∼ 10⁵ M_⊙ at the virial radius, that during their infall lose more than 99% of their mass due to the hydrodynamical interactions with the hot coronal gas, implying that the internal regions of the halos are almost devoid of cool gas;

3. despite the uncertainties given by our approximation in the hydrodynamical treatment, we can conclude that the cool CGM clouds most probably cannot survive the journey but will completely evaporate into the hot corona as the timescale of thermal conduction is much shorter than the dynamical time. Thus the cool gas will probably not accrete onto the central galaxies, providing an explanation of their quiescence.

The results of this paper represent a step forward in solving the puzzle of the existence of large unused reservoirs of cool circumgalactic gas around quiescent galaxies. Future improvement of the observational data and ad-hoc hydrodynamical simulations could reduce the uncertainties in our current description. We also plan to apply similar models to larger samples of data and in particular to star forming galaxies, in order to study the impact of the CGM kinematics and dynamics on the evolution of different kinds of galaxies.

Acknowledgments

The authors would like to thank the referee Sean Johnson for a very constructive report and the helpful comments and suggestions. We are also very grateful to Lorenzo Posti and Antonino Marasco for the useful discussions and suggestions that significantly helped to improve this work. GP acknowledges support by the Swiss National Science Foundation, grant PP00P2_163824.

References

Appendix A: Geometry of the system

In Sect. 4 we have described how we compared our results with the observations, using synthetic observations of random populations of clouds. Here we report in detail the geometry used to infer the number of observed clouds and their line-of-sight velocities. We created for each galaxy a 3-dimensional distribution of clouds over a sphere with radius r_vir (see Table 1), representing the virial halo. As already explained in Sect. 4, the total number of clouds for each halo is calculated using Eq. (16), while the intrinsic distance of each cloud from the galaxy is found using the probability density function expressed by Eq. (17). We then associated to each cloud two other coordinates θ and ϕ. These last two coordinates are chosen to have the clouds uniformly distributed over the sphere: we created a random distribution of the angle ϕ uniformly distributed between 0 and 2π and a random distribution of cosθ uniformly distributed between −1 and 1, with 0 <  θ <  π. Each cloud has therefore a specific position in the halo, expressed by the coordinates (r, θ, ϕ).

Figure A.1 shows the method used to perform the synthetic observations. Here, the line of sight is given by the line at x = R in the plane ϕ = 0, where R is the distance of the galaxy from the sightline, reported for each object in Table 1. The cloud in (r, θ, ϕ) is observed if the distance d of this point from the line of sight is smaller than the cloud radius, expressed by

$$\begin{array}{c}\hfill d<r_{\mathrm{cl}},\end{array}$$(A.1)

where

$$\begin{array}{c}\hfill d=\sqrt{{(R-r_{\mathrm{i}}sin\theta cos\varphi )}^2+{(r_{\mathrm{i}}sin\theta sin\varphi )}^2}.\end{array}$$(A.2)

If the cloud is observed, that is if condition Eq. (A.1) is true, we calculated its line-of-sight velocity, as shown in Fig. A.1, through

$$\begin{array}{c}\hfill v_{\mathrm{los}}=vcos\theta \end{array}$$(A.3)

and we added it to the model velocity distribution to be compared with the COS-LRG observations.

Fig. A.1. 3-dimensional representation of the synthetic observation of a generic cloud, as performed in our model.

Fig. A.2. Corner plot with the results of the MCMC analysis described in Appendix B, where also the starting velocity of the clouds is added as a fourth free parameter. Both the one and two dimensional projections of the posterior probabilities of the four free parameters are shown.

Appendix B: Cloud initial velocity

In Sect. 3 we have described the construction of our semi-analytical models and we have argued that the starting velocity of the clouds (v_start, the velocity of the clouds at the virial radius of the galaxy) is needed to solve Eq. (12). In this Appendix we assess models where v_start is not fixed, but free to vary between different values. In order to do this, we have run an MCMC analysis as in Sect. 4, but adding v_start as a fourth free parameter, with a flat prior from 0 to 350 km s⁻¹. Figure A.2 shows the results of this analysis, the one and two dimensional projections of the posterior probabilities of the four parameters. The distributions of the first three parameters are in strong agreement with the ones in Fig. 5, meaning that adding the initial velocity as a free parameter does not change the final results.

The four bottom panels of Fig. A.2 represent the posterior distributions regarding v_start. We note that models with clouds starting from the virial radius with high velocities are strongly disfavored by the MCMC analysis and that the distributions are instead peaked on very low values. This means that the best models in reproducing the COS-LRG observations have clouds with initial velocities close to zero. This condition is indeed needed by the models to reproduce the central region of the line-of-sight velocity distribution, which shows an excess of clouds with very low velocities (see Fig. 1 and the left panel of Fig. 6). For this reason, we have used in our work only three free parameters, fixing the starting velocity of the cool CGM clouds at zero.

All Tables

All Figures

	Fig. 1. Orange: observed velocities of the cool CGM clouds relative to systemic velocities of the galaxies, obtained using the COS-LRG data for a sample of 16 massive passive galaxies (their properties are reported in Table 1). The light-blue dashed line is a gaussian with a dispersion equal to the typical 1D (line-of-sight) virial velocity dispersion value in the sample.
In the text

	Fig. 2. Velocity profiles of the cool CGM clouds obtained using the ballistic model (Eq. (3), blue line) and the final model that includes all the hydrodynamical interactions (Eq. (12)), with three different choices of mcl, start and α.
In the text

Fig. 3. Red: profile of the intrinsic density of the hot coronal gas (from Eq. (6)) for a galaxy with Mvir = 1013.3 M⊙, rvir = 526 kpc, T = Tvir = 6.5 × 106 K (average of the galaxy halos in our sample), and normalization set as explained in Sect. 3.2. In blue, profile of the intrinsic average density expected for the cool circumgalactic medium if in pressure equilibrium with the corona, obtained through Eq. (8).

In the text

Fig. 4. Cartoon 2D representation of the synthetic observations performed in our model. The clouds (blue circles) are distributed throughout the halo (the sphere with radius rvir), all accreting towards the central galaxy (the central red star) with a radial velocity described by Eq. (12). The vertical gray line represents the line of sight, consistent with the projected distance R of the galaxy from the background quasar, while the area included between the two vertical dashed lines represents the range spanned by the COS-LRG observations. The observer is located at the bottom side of the panel.

In the text

	Fig. 5. Corner plot with the MCMC results, representing the one and two dimensional projections of the posterior probabilities for the three free parameters of our models.
In the text

	Fig. 6. Comparison between the observations (orange bars and circles) and our best model results (blue lines). Left panel: normalized line-of-sight velocity distribution, where the errors in the observations are obtained with the bootstrapping method. Right panel: total hydrogen column density. The downward arrows represent the upper limits of the non-detections.
In the text

Fig. 7. Properties of the cool CGM clouds described by our best model. All the profiles are plotted with respect to the intrinsic distance from the central galaxies and the different colored lines represent the result for the sixteen different galaxies of the sample. Top left: velocity profiles of the cool clouds. Top right and bottom left: respectively masses and radii of the cool clouds. Bottom right: cool CGM cumulative mass profile, with the vertical dashed line representing the distance of 160 kpc.

In the text

	Fig. 8. Same as Fig. 7, but for the total mass flux profile of the cool CGM.
In the text

	Fig. 9. Cumulative distribution functions of the line-of-sight velocity distributions of the observations (orange curve) and of our models. In blue we report the result of our main model (Sect. 5), while the dashed curves represent models with rstart = rvir/2 (green, dashed) and rdestr = rvir/2 (red, small dashed).
In the text

	Fig. A.1. 3-dimensional representation of the synthetic observation of a generic cloud, as performed in our model.
In the text

	Fig. A.2. Corner plot with the results of the MCMC analysis described in Appendix B, where also the starting velocity of the clouds is added as a fourth free parameter. Both the one and two dimensional projections of the posterior probabilities of the four free parameters are shown.
In the text