The formation of hot gaseous haloes around galaxies
We use a suite of hydrodynamical cosmological simulations from the Evolution and Assembly of GaLaxies and their Environments (EAGLE) project to investigate the formation of hot hydrostatic haloes and their dependence on feedback mechanisms. We find that the appearance of a strong bimodality in the probability density function (PDF) of the ratio of the radiative cooling and dynamical times for halo gas provides a clear signature of the formation of a hot corona. Haloes of total mass develop a hot corona independent of redshift, at least in the interval where the simulation has sufficiently good statistics. We analyse the build up of the hot gas mass in the halo, , as a function of halo mass and redshift and find that while more energetic galactic winds powered by SNe increases , AGN feedback reduces it by ejecting gas from the halo. We also study the thermal properties of gas accreting onto haloes and measure the fraction of shock-heated gas as a function of redshift and halo mass. We develop analytic and semianalytic approaches to estimate a ‘critical halo mass’, , for hot halo formation. We find that the mass for which the heating rate produced by accretion shocks equals the radiative cooling rate, reproduces the mass above which haloes develop a significant hot atmosphere. This yields a mass estimate of at , which agrees with the simulation results. The value of depends more strongly on the cooling rate than on any of the feedback parameters.
keywords:galaxies: haloes – galaxies: formation – galaxies: evolution – methods: numerical – methods: analytical
One of the major goals of modern galaxy formation theory is to understand the physical mechanisms that halt the star formation process, by either removing, heating or preventing the infall of cold gas onto the galactic disc. X-ray observations suggest that for haloes hosting massive galaxies the majority of baryonic matter resides not in the galaxies, but in the halo in the form of virialized hot gas (e.g. Lin03; Crain10a; Anderson11). This work investigates the formation of the hot gaseous corona (also refereed to as ‘hot halo’ or ‘hot atmosphere’) around galaxies, that may help reduce the rate of infall of gas onto galaxies, and has been suggested to explain the observed galaxy bimodality (Dekel).
The hot gaseous corona is produced as a result of an important heating process, that was initially discussed by Rees, Silk77, Binney and White, and later in the context of the cold dark matter paradigm by e.g. White91, in an attempt to explain the reduced efficiency of star formation within massive haloes. They proposed that while a dark matter halo relaxes to virial equilibrium, gas falling into it experiences a shock, and determined the cooling time of gas behind the shock. As long as the cooling time is shorter than the dynamical time, the infalling gas cools (inside the current ‘cooling radius’) and settles onto the galaxy. If, on the other hand, the cooling time exceeds the dynamical time, the gas is not able to radiate away the thermal energy that supports it. It therefore adjusts its density and temperature quasi statically, forming a hot hydrostatic halo atmosphere, pressure supported against gravitational collapse. Over the past decade, the works of Birnboim and Dekel investigated the stability of accretion shocks around galaxies, and concluded that a hot atmosphere forms when the compression time of shocked gas is larger than its cooling time, occurring when haloes reach a mass of about .
Numerical simulations have shown, however, that cold gas accreting through filaments does not necessarily experience a shock when crossing the virial radius, even if the spherically averaged cooling radius is smaller than the virial radius. Many groups have concluded that there are two modes of gas accretion, named as hot and cold accretion, that are able to coexist in high-mass haloes at high redshift (e.g. Keres05; Dekel; Ocvirk; Dekel09; vandeVoort11; Faucher; vandeVoort12; Nelson13). The hot mode of accretion refers to the accreted gas that shock-heats to the halo virial temperature. The cold mode refers to gas that flows along dark matter filaments and is accreted onto the central galaxy without being shock heated near the virial radius. It has been found that the cold streams end up being the dominant mode of accretion onto galaxies at high redshift (e.g. Dekel09). However, it has also been found that most of the cold gas from filaments does experience significant heating when accreted by the galaxy at radii much smaller than the virial radius (Nelson13).
Besides the rate of gas accretion, the hot halo can be influenced by feedback mechanisms and photoionization from local sources. Feedback mechanisms can suppress cooling from the hot halo, modify the distribution of hot gas in the halo (vandeVoort12) and (to a limited degree) reduce the accretion rates onto haloes (vandeVoort11; Faucher; Nelson15a). In this work we investigate the impact of feedback mechanisms on the hot halo in detail and analyze whether reasonable changes to the feedback implementation result in a change to the mass scale of hot halo formation. Increasing photoionizing flux (higher star formation rate or an active nucleus) from local sources can decrease the net cooling rate of gas in the proximity of the galaxy, potentially suppressing cold gas accretion in low-mass halos () and decrease the mass-scale for hot halo formation (Cantalupo10). However, the results are sensitive to the assumed escape fraction and Vogelsberger12 found only small effects when including local AGN. For simplicity we will assume the gas is only exposed to the metagalactic background radiation.
We use the suite of cosmological hydrodynamical simulations from the EAGLE project (Schaye14; Crain15) to investigate the physical properties of the hot gas in the halo, and their dependence on energy sources like stellar feedback and AGN feedback. Our main goal is to study the thermal properties of gas accreting onto haloes and the gas mass that remains hot in the halo (). In addition, we develop analytic and semianalytic approaches to calculate the heating rates of gas in the halo and the mass scale of hot halo formation, which we apply in a companion work, Correa et al. (in preparation, hereafter Paper II). In Paper II we derive a physically motivated model for gas accretion onto galaxies that accounts for the hot/cold modes of accretion onto haloes, and for the rate of gas cooling from the hot halo. With this model we aim to provide some insight into the physical mechanisms that drive the gas inflow rates onto galaxies.
The outline of this paper is as follows. We describe the EAGLE simulations series used in this study and the analysis methodology in Section 2. We present our main results concerning the physical properties of hot and cold gas in the halo in Section 3 and on the modes of gas accretion in Section 4. In Section 5 we develop an analytic approach to calculate a ‘critical mass scale’, , for hot halo formation, and compare it with our numerical results and previous works. Finally, in Section 6 we summarize our conclusions.
|(comoving )||()||()||(comoving )||(proper )|
|Less Energetic FB||0.67|
|More Energetic FB||0.67|
|No AGN FB||0.67|
|More Explosive AGN FB||0.67|
To investigate the formation and evolution of hot haloes surrounding galaxies, we use cosmological, hydrodynamical simulations from the Evolution and Assembly of GaLaxies and their Environments project (EAGLE; Schaye14; Crain15). The EAGLE simulations were run using a modified version of GADGET 3 (Springelb), a -Body Tree-PM smoothed particle hydrodynamics (SPH) code. The EAGLE version contains a new formulation of SPH, new time stepping and new subgrid physics. Below we present a summary of the EAGLE models. For a more complete description see Schaye14.
The EAGLE simulations assume a CDM cosmology with the parameters derived from Planck-1 data (Planck), , , , , . The primordial mass fractions of hydrogen and helium are and , respectively.
Table 1 lists the box sizes and resolutions of the simulations used in this work. We use the notation xxxyyyy, where xxx indicates box size (ranging from 25 to 100 comoving ) and yyyy indicates the cube root of the number of particles per species (ranging from to , with the number of baryonic particles initially equal to the number of dark mater particles). The gravitational softening was kept fixed in comoving units down to and in proper units thereafter. We will refer to simulations with the mass and spatial resolution of L025N0376 as intermediate-resolution runs, and to simulations with the resolution of L025N0752 as high-resolution runs.
2.1 Baryonic physics
Radiative cooling and photo-heating are included as in Wiersma09a. The element-by-element radiative rates are computed in the presence of the cosmic microwave background (CMB), and the Haardt model for UV and X-ray background radiation from quasars and galaxies.
Star formation is modelled following the recipe of Schaye08. Star formation is stochastic above a density threshold, , that depends on metallicity (in the model of Schaye04, is the density of the warm, atomic phase just before it becomes multiphase with a cold, molecular component), with the probability of forming stars depending on the gas pressure. The implementation of stellar evolution and mass loss follows the work of Wiersma09b. Star particles are treated as simple stellar populations with a Chabrier initial mass function, spanning the range . Feedback from star formation and supernovae events follows the stochastic thermal feedback scheme of DallaVecchia12. Rather than heating all neighboring gas particles within the SPH kernel, they are selected stochastically based on the available energy, and then heated by a fixed temperature increment of K. The probability that a neighboring SPH particle is heated is determined by the fraction of the energy budget that is available for feedback, . If is sufficiently high to ensure that radiative losses are initially small, the physical efficiency of feedback can be controlled by adjusting . The value corresponds to the expected value of energy injected by core collapse supernovae ( erg per solar mass of stars formed). EAGLE takes to be a function of the local physical conditions,
which depends on maximum and minimum threshold values ( and , respectively), on density ( refers to hydrogen number density and to the density inherited by the star particle) and metallicity () of the gas particle. The reference simulations (hereafter Ref) use , and . These values were chosen to obtain good agreement with the observed present-day galaxy stellar mass function and disc galaxy sizes (as described by Crain15).
Black hole seeds (of mass ) are included in the gas particle with the highest density in haloes of mass greater than that do not contain black holes (Springel05). Black holes can grow through mergers and gas accretion. The accretion events follow a modified Bondi-Hoyle formula that accounts for the angular momentum of the accreting gas (Rosas13; Schaye14), and a free parameter that is related to the disc viscosity. AGN feedback follows the accretion of mass onto the black hole, where a fraction (0.015) of the accreted rest mass energy is released as thermal energy into the surrounding gas, and is implemented stochastically, as per the stellar feedback scheme, with a fixed free parameter heating temperature, , which is set to K in the reference simulations.
When the resolution is increased, the parameters may need to be (re-)calibrated to retain the agreement with observations. The high-resolution simulation with recalibrated parameters is called Recal. In addition to Ref and Recal, we also use simulations with different feedback implementations to test the impact of feedback on the formation of the hot halo. Table 2 lists the values of the feedback parameters adopted in each simulation. In the table, the simulation identifier describes the differences in the feedback with respect to Ref. In the stellar feedback case, ‘Less/More Energetic FB’ means that in these simulations, the energy injected per mass of stars formed is lower/higher with respect to Ref. In the AGN case, ‘More Explosive AGN FB’ means that AGN feedback is more explosive and intermittent, but the energy injected per unit mass accreted by the BH does not change with respect to Ref. Additional information regarding the performance of the EAGLE simulations, including an analysis of subgrid parameter variations, a study of the evolution of galaxy masses, star formation rates and sizes can be found in Crain15, Furlong15a; Furlong15b and Schaye14.
There has been much debate regarding the systematic differences between SPH, grid codes and moving mesh grid codes when modelling fluid mixing and gas heating and cooling (e.g. Agertz07; Vogelsberger12; Nelson13). It has been shown by Hutchings and Creasey that SPH simulations may not adequately resolve shocks of accreted gas. Since shocks are generally spread over several SPH kernel lengths, the heating rate is smoothed over time, potentially making it easier for radiative cooling to become important. In addition, if radiative cooling is able to limit the maximum temperature reached by the gas particle, numerical radiative losses may be enhanced.
In contrast, numerical simulations using grids do not smooth out the shocks, and are thus better at identifying shock temperatures spikes. Numerical simulations using moving mesh codes can also capture shocks accurately. However, in common with grid codes, they may suffer from numerical mixing of hot and cold gas as the fluid moves across cells. Recently, Nelson13 compared the moving mesh code AREPO (AREPO) with the standard SPH version of GADGET, and calculated the rates of cold mode of gas accretion onto haloes and galaxies. They found that while the rates of gas accreted cold onto haloes are in very good agreement between the simulations run with GADGET and AREPO, the rates of gas accreted cold onto galaxies differ significantly, with galaxies in AREPO having a lower cold fraction in haloes. Nelson13 concluded that most of the cold gas from filaments experiences significant heating after crossing the virial radius, implying that the numerical deficiencies inherent in different simulation codes may modify the relative contributions of hot and cold modes of accretion onto galaxies.
Some differences in the contributions of hot and cold modes of accretion onto galaxies and haloes may, however, be due to the method employed to select shock-heated gas. Previous works (e.g. Keres05; Keres09; Faucher; vandeVoort11; Nelson13, among others) followed the thermal history of the gas and applied a fixed temperature cut to the distribution of the maximum past temperature (), to separate hot from cold mode accretion. However, is not suitable for identifying cold flows if the gas experiences a shock but cools immediately afterwards, as may happen for accretion onto galaxies. In this case, a filament that is mostly cold except at a point near the galaxy would be labeled as hot mode accretion by numerical studies using , but observers would identify it as a cold flow. This practical problem may not be important for SPH simulations that suffer from ‘in shock cooling’ because they do not resolve the accretion shocks onto the galaxy, or, as in the case of vandeVoort11 and EAGLE, that impose a temperature-density relation onto high-density gas, but it may affect the conclusions inferred from moving mesh codes using the statistic. To avoid this issue, we use an alternative method to identify shock-heated particles in Section 4, based on post-shock temperature values.
Hydrodynamics solvers may also produce differences in the hot/cold modes of accretion. The EAGLE version of GADGET uses the hydrodynamics solver “Anarchy”, which greatly improves the performance on standard hydrodynamical tests, when compared to the original SPH implementation in GADGET (Schaller15b, see Hu14 for similar results). Anarchy makes use of the pressure-entropy formulation derived in Hopkins13, alleviating spurious jumps at contact discontinuities. It also uses an artificial viscosity switch advocated by Cullen10, that allows the viscosity limiter to be stronger when shocks and shear flows are present. In addition, Anarchy includes an artificial conduction switch (similar to that of Price08), the Wendland95 kernel and the time step limiters of Durier12. These changes ensure that ambient particles do not remain inactive when a shock is approaching.
Recently, Sembolini16a compared cosmological simulations of clusters using SPH as well as mesh-based codes. They found that the modern SPH schemes (such as Anarchy) that allow entropy mixing produce gas entropy profiles that are indistinguishable from those obtained with grid-based schemes. In addition, Schaller15b compared the EAGLE simulations with simulations run with the same subgrid physics, but using the standard GADGET rather than the Anarchy hydrodynamics solver. They found that while simulations with standard SPH contain haloes with a large number of dense clumps of gas at all radii, Anarchy’s ability to mix phases allows dense clumps to dissolve into the hot halo. These substantial improvements of the SPH formulation in the EAGLE simulations motivate a detailed description of the resulting predictions for hot halo formation and of hot/cold mode accretion.
2.3 Identifying haloes and galaxies
Throughout this work we select the largest subhalo in each Friends-of-Friends (FoF) group, and use the SUBFIND algorithm (Springel; Dolag) to identify the substructures (subhaloes) within it. The FoF algorithm adopts a dimensionless linking length of 0.2, and the SUBFIND algorithm calculates halo virial masses and radii via a spherical overdensity routine that centers the main subhalo from the FoF group on the minimum of the gravitational potential. We define halo masses, , as the mass of all matter within the radius, , for which the mean internal density is 200 times the critical density of the Universe.
To select the gas associated with the central galaxies embedded in each resolved halo, we identify the gravitationally bound cold and dense gas within that is star-forming and/or has a hydrogen number density, , and temperature K. We also require all particles to be contained within a sphere of radius , in order to avoid labelling infalling cold flows (that would be included by the cuts but are mostly at large radii) as part of the galaxy.
2.4 Measuring gas accretion
Once we have identified the haloes, we build merger trees across the simulation snapshots111The simulation data is saved in 10 discrete output redshifts between redshift 0 to 1, in 4 output redshifts between redshift 1 and 3, and in 8 output redshifts between redshift 3 and 8.. The standard procedure to build a halo merger tree is to link each progenitor halo with a unique âdescendantâ halo in the subsequent output (see e.g. Fakhouri). To do so, we identify the main branches of the halo merger trees and compute the halo (and central galaxy) accretion rate between two consecutive snapshots. At each output redshift (snapshot) we select the most massive haloes within each FoF group and consider them to be ‘resolved’ if they contain more than 1000 dark matter particles, which corresponds to a minimum halo mass of () in the intermediate- (high-) resolution simulations. This limit on the number of dark matter particles results from a convergence analysis that we present in Paper II, where we find that in smaller haloes the accretion onto galaxies does not converge, indicating that the inner galaxies are not well resolved. We refer to the resolved haloes as ‘descendants’, and then link each descendant with a unique ‘progenitor’ at the previous output redshift. This is nontrivial due to halo fragmentation, in which subhaloes of a progenitor halo may have descendants that reside in more than one halo. Such fragmentation can be either spurious or due to a physical unbinding event. To account for this, we link the descendant to the progenitor that contains the majority of the descendant’s 25 most bound dark matter particles (see Correa15b for an analysis of halo mass history convergence using these criteria to connect haloes between snapshots).
We distinguish between gas accreted onto a halo and gas accreted onto a galaxy. For each descendant halo at and its linked progenitor at (), we identify the particles that are in the descendant but not in its progenitor by performing particle ID matching. We then select particles that are new in the halo and reside within the virial radius, as particles accreted onto the halo in the redshift range . The accretion rate onto galaxies is further explored in Paper II, where we follow the methodology described above for calculating accretion rates onto haloes, and we select the new particles within the radius as particles accreted onto the galaxies in the redshift range (see Paper II, Section 2.1, for a discussion on methods for calculating gas accretion onto galaxies).
3 Hot halo formation
The simple models of galaxy formation (e.g. Rees and White) assume that as long as the cooling time, , is shorter than the dynamical time, , the infalling gas cools (inside a ‘cooling radius’, White91) and settles into the galaxy. Otherwise, the gas is unable to efficiently radiate its thermal energy and forms a hot hydrostatic atmosphere, which is pressure supported against gravitational collapse. More recent semi-analytic models of galaxy formation assume that the cooling radius expands outwards as a function of time, therefore the comparison is done between gas cooling time and a different time scale representing the time available for cooling, like the time since the last major event or the time of halo formation (see e.g Lacey16). In this section we investigate when the hot hydrostatic halo forms in the EAGLE simulations, by analyzing the interplay between the cooling and dynamical times of the gas particles in the halo. Throughout this work we define hot halo gas as all gas particles that have , and that do not form part of the galaxy, i.e. .
We calculate of the gas particle as
where is the circular velocity and is the mass enclosed within . We calculate as
where is the number density of the gas particle (, is the molecular mass weight calculated from the cooling tables of Wiersma09b, and is the proton mass), is the Boltzmann constant, is the gas temperature and is the cooling rate per unit volume with units of erg cms. To calculate , we use the tabulated cooling function for gas exposed to the evolving UV/X-ray background from Haardt given by Wiersma09b, which was also used by the EAGLE simulations. Note that the “standard” definition for dynamical time of gas within a virialized system depends on and , and not on the local radius and local circular velocity as defined here. We use local values rather than to investigate if shorter dynamical times in the inner dense regions give rise to a cooling flow. However, we find that changing eq. (2) to does not change our conclusions.
Fig. 1 shows temperature profiles (left column), the logarithm of the ratio between cooling times and dynamical times (middle column) and the respective mass-weighted probability density function (PDF) of (right column) for gas from haloes in the mass range (top row), (middle row) and (bottom row) at , taken from the Ref-L025N0752 simulation. In the right panels and throughout this work, the PDFs are calculated by stacking haloes in the selected mass range and distributing the gas particles in logarithmic bins of size 0.1 dex. We then sum the mass of the gas particles in each bin and normalize the distribution by the total gas mass. In the left panels, the legend lists the median values of the mass, virial temperature (defined as , with ) and virial radius of haloes selected in each mass bin. The left panels also show in solid, dotted and dashed lines the median temperature per radial bin of gas from haloes taken from the simulations Ref-L025N0376, Ref-L025N0752 and Recal-L025N0752, respectively.
The left panels of Fig. 1 show that while there is very good agreement in the median temperature profiles at small () and large () radii, at intermediate radii the median temperatures from the intermediate-resolution run (Ref-L025N0376) are larger by up to 0.3 dex than those from the high-resolution runs (L025N0752). This is in agreement with the convergence analysis of Nelson15b, who concluded that the physics (different models of stellar winds or AGN feedback) has a greater impact on than resolution. We also find, that in the radial range , where the convergence with resolution is poorest, the median temperatures drop from to K (in agreement with Nelson15b, and vandeVoort12, ), because of the high densities, that rapidly decrease the gas cooling times, enabling it to radiate away its thermal energy and join the ISM.
The top and middle left panels of Fig. 1 show that there is relatively little gas with K at small and intermediate radii, reflecting the short cooling times at these temperatures. The cooling flow in the hot halo is formed by K gas that slowly decreases in temperature as it loses hydrostatic support due to cooling. The ISM consists of K gas at . For a better understanding of the hot halo forming as a function of halo mass and its effect on the infall rate of gas onto the galaxy, we next analyze the middle and right columns. The bottom middle panel shows that in haloes with masses between , most of the gas has low temperature (), short cooling times () and infalls towards the central galaxy, although a substantial fraction of gas has at . At larger halo masses, a larger fraction of the gas is unable to cool and therefore forms a hot halo. The middle panel shows that haloes between are in the intermediate stage between developing a hot stable atmosphere (gas with and ) and continuing to fuel the galaxy. The top middle panel clearly shows a stable hot halo and a reduced amount of gas infalling towards the galaxies (gas at and ).
Fig. 2 is similar to Fig. 1, but shows haloes in the mass range (top panels), (bottom panels) at . The top middle panel shows that haloes develop a hot atmosphere, despite the significant fraction of cold gas at large radii that is accreted onto the halo. This cold gas forms part of the cold filamentary flows, that cross the virial radius, and are directly accreted onto the central galaxy. The cold accretion mode is best seen as the gas at K and in the bottom panel of Fig. 2, which remains cool (with K), as it is accreted onto the galaxy.
The presence of cold flows produces gaseous haloes with an isothermal temperature profile of K at all radii. Besides analyzing the cooling time profiles, Nelson15b and vandeVoort12 analyzed the entropy profiles of haloes at , and concluded the that while the entropy of the cold-mode gas decreases smoothly and strongly towards the centre, the entropy of the hot-mode gas decreases slightly down to , after which it drops steeply. We find that the cooling time profiles of the hot () and cold () modes follow the entropy profiles.
Figs. 1 and 2 show that as the halo mass increases, so does the hot gas mass. In the case of haloes, the bottom left panel of Fig. 1 shows that there is a large fraction of gas at with temperatures equal to or larger than the halo virial temperature. This seems to indicate that gas is shock-heating to the halo virial temperature when crossing and forming a hot atmosphere. However, gas with in the outer parts of haloes does not imply that the halo formed a stable hot atmosphere via gravitational accretion shocks, since the gas is affected by the extragalactic UV/X-ray background as we show in the next section.
The figures also show that as haloes are growing a hot atmosphere, the PDF begins to present a strong bimodal shape, with a local maximum at followed by a local minimum at (see top right panel). We then conclude that the bimodality in the cooling time PDF provides a clear signature of the formation of a hot halo, and the right panels indicate that the hot hydrostatic halo forms in the halo mass range to with a weak dependence on redshift.
We also analyse the radial velocity distributions of gas and find that in haloes of mass at , () of hot gas has an absolute radial velocity lower than 100 km/s (50 km/s), indicating that most of it is in hydrostatic equilibrium and not accreting onto the galaxy.
3.1 The impact of photoheating in low-mass haloes
In the previous section we analyzed the dependence of the gas cooling rates on halo mass and concluded that a halo with a hot hydrostatic atmosphere should present a strongly bimodal PDF with a local maximum at . As the virial temperature decreases from K to K, we would naively expect the peak in the PDF to shift towards shorter cooling times. This is however not the case: we find that at the gas in haloes less massive than ( K) is affected by the extragalactic UV/X-ray background radiation, which strongly suppresses the net cooling rate of gas in the temperature range K (Efstathiou92; Wiersma09b). As a result, the peak in the PDF remains at . This can be seen in Fig. 3, where we show the PDF of gas from haloes in the mass range (olive lines), (orange lines) and (red lines). The solid lines correspond to the case where the cooling rates are calculated for gas exposed to the evolving UV/X-ray background from Haardt, while the dashed lines correspond to the case where the cooling rates are calculated for gas in collisional ionization equilibrium (CIE) and not exposed to the background radiation. Note, however, that we apply these CIE cooling rates to simulations that were run using cooling rates that did account for photoheating, limiting the gas temperature to K.
Fig. 3 shows that there is no large difference in the gas PDF for haloes more massive than , indicating that there is no strong impact of the background radiation on the cooling rates of gas from haloes with virial temperatures larger than K. In smaller haloes, the peak in the PDF curve is shifted to in the case of no background radiation.
3.2 The impact of feedback
Feedback can affect the formation of the hot hydrostatic halo around galaxies. For example, very energetic SN activity generates large outflows and strong winds, that shock against the gaseous halo. As a result, the winds can fill the halo with gas expelled from the galaxy, increasing the amount of hot gas at large radii. In this subsection we compare the mass-weighted PDF at fixed halo mass obtained from simulations with different feedback implementations (see Table 2 for reference and section 2.1 for a brief description), and determine, by analysing whether the cooling time PDF is bimodal, the mass range where the hot halo is forming.
Fig. 4 shows the mass-weighted PDF of of gas from haloes in the Ref-L100N1504 simulation in the mass range , , and at (left panel) and (right panel). In the figures, the labels show the median halo mass for each mass bin. The region where corresponds to hot gas in the halo.
As the halo mass increases so does the amount of hot gas (Section 3.3). The PDF of gas in haloes at shows a bimodal shape, that becomes mostly unimodal in higher-mass haloes. At , the bimodality persists up to the highest mass bin () due to the contribution from cold flows that populate the peak at . We find that the presence of the bimodality in the PDF indicates the increasing amount of hot gas at large radii and the eventual formation of the hot halo. Then, from visual inspection, we determine that the hot hydrostatic atmosphere is forming in haloes with masses between and at and .
The panels in Fig. 5 repeat the analysis shown in the left panel of Fig. 4, but instead show mass-weighted PDFs for the L025N0376 simulations with different feedback prescriptions. In this case, the PDFs correspond to gas from haloes in the mass range , , and . In the panels, the top left legends indicate the total gas mass in the halo ().
The simulation shown in the top left panel of Fig. 5 does not have AGN feedback while the one in the bottom left panel uses a more explosive AGN feedback. Both include the same feedback from star formation as in Ref. It can be seen that neither the bimodality of the PDF nor the amount of hot gas are strongly affected by AGN feedback in haloes with masses between . The right panels in Fig. 5 show the PDFs in the less (top panel) and more energetic stellar feedback (bottom) scenarios (both including the same AGN feedback as in the Ref model). For these halo masses, stellar feedback has a strong impact on the PDFs. While a more energetic stellar feedback increases the fraction of hot gas, at least in the halo mass range probed by these simulations () and thus limits the build-up of cold-mode gas in the halo centre (in agreement with vandeVoort12), a less energetic stellar feedback maintains the bimodality in the PDF but shifts the peak in the PDF of hot gas towards larger cooling times. We find that the bimodality of the PDF is present in haloes with more energetic stellar feedback and in haloes with less energetic stellar feedback.
In the following section we further analyse the dependence of the total hot gas mass on halo mass, redshift and feedback. In Section 5 we derive an analytic model for the formation of a stable hot hydrostatic atmosphere. In the model we calculate a halo mass scale for which the gravitational heating rate of the hot halo gas balances the gas cooling rate, thus keeping the gas hot and enabling the formation of a hot atmosphere. With the model we show that the ability of a halo to develop a hot hydrostatic atmosphere depends on the amount of hot gas that the halo already contains, which we calculate in the next subsection, and on the fraction of accreted gas that shock-heats to the halo virial temperature (Section 4).
3.3 Hot gas mass
In order to better understand the build up of the hot gas mass, (mass of gas with and ) as haloes evolve, in this section we look for a correlation between and the total halo mass as a function of redshift. Fig. 6 shows the median ratio of (with the universal baryon fraction) taken from a range of simulations (as indicated in the legends) at (top panel) and at (second panel from the top). In these panels the error bars show the scatter for the Ref-L025N0752 and Ref-L100N1504 simulations. The median ratio of is also shown in the third panel from the top, but in this case the values are taken from the Ref-L100N1504 simulation and at various output redshifts.
The top panel of Fig. 6 highlights the relatively poor agreement between the intermediate- and high-resolution simulations, with the latter predicting somewhat higher hot gas fractions. Good agreement is however achieved at (middle panel). Although the Ref-L100N1504 simulation is not fully converged with respect to the numerical resolution at , the convergence with box size is excellent at all redshifts. The intermediate-resolution runs show that the hot gas represents of the total halo gas mass for at . The hot gas mass fraction reaches in haloes and remains roughly constant for higher masses. In very low-mass haloes (), the hot gas mass fraction also remains roughly constant (). In these haloes cold accretion dominates, therefore the heating mechanism that maintains is the UV background as discussed in Section 3.1.
The third from the top panel of Fig. 6 shows the evolution of the hot gas fraction. In haloes larger than , remains constant over the redshift range 3 to 6 and at lower redshift it increases somewhat with time. In smaller haloes (), increases with time until but decreases thereafter. We calculate the cooling rate of gas exposed to the UV background, and in the absence of it and compare the hot gas mass. We find that the hot gas mass in low-mass haloes increases due to the heating produced by the background radiation. In the case of gas not being exposed to the UV background, the total hot gas mass decreases with increasing redshift at fixed halo mass. We also find that the differences between occurs in haloes lower than .
We next perform a least-square minimization to determine the best-fit relation as a function of redshift. We apply equal weighting for each mass bin from the Ref-L100N1504 simulation (which we use to cover a large halo mass range) and minimize the quantity , where
is the number of bins at each output redshift , and is
We obtain the best-fitting values for , and at each redshift , and following the same methodology we look for the best-fit expression of these parameters as functions of redshift.
We find that the following expression best reproduces the relation in the halo mass range ,
where , and are functions of given by
where . The bottom panel of Fig. 6 shows the residual of the data points with respect to the best-fit expression.
Next, we investigate how the presence of different feedback mechanisms affect the hot gas as well as the total gas mass () in the halo (all gas contained between ). The top panel of Fig. 7 shows the relation for haloes in the mass range at for the L025N0376 simulations. The different coloured lines correspond to simulations with different feedback prescriptions.
This panel shows that the impact of feedback increases with halo mass and that stellar feedback has a larger impact on the amount of gas in the halo than AGN feedback. We find that doubling the efficiency of the stellar feedback increases the gas mass fraction by a factor of 1.3 in haloes relative to the Ref model, whereas halving the efficiency decreases the gas mass fraction by a factor of 2.5. No (Explosive) AGN feedback results in an increase (decrease) by a factor of 1.5 in the gas mass fraction.
While efficient stellar feedback increases the gas mass in the halo, more explosive AGN feedback decreases it. Overall, it seems that in haloes more massive than there is a greater difference in the gas mass fraction between Ref and More Energetic FB than between Ref and More Explosive AGN FB. This is due to two different reasons. Physically, AGN feedback mainly ejects gas mass from the halo, or prevents it from falling into the halo, whereas stellar feedback ejects gas out of the galaxy into the inner halo. Numerically, although stellar and AGN feedback use a similar thermal implementation (DallaVecchia12), there is a difference in the actual energetics of the processes. The energy injected per mass of stars formed changes between Ref and More Energetic FB, whereas the energy injected per unit mass accreted by the BH does not change between Ref and More Explosive AGN FB. In the latter, it is only the intermittency and the explosiveness that changes as a consequence of the change in the temperature of the AGN. In the case of Less Energetic FB, we find that the gas mass in the halo decreases because more gas is accreted by the galaxy and locked up in stars.
The bottom panel of Fig. 7 shows the variation of the ratio with feedback. It can be seen that in haloes less massive than , the ratio increases with decreasing halo mass, indicating that most of the halo gas is heated by the X-ray/UV background (see Section 3.1 for a discussion). In haloes more massive than , increases with halo mass. While a more energetic stellar feedback increases the hot mass fraction by 10, no AGN feedback decreases it by in haloes. In the case of Less Energetic FB and More Explosive AGN, increases by and (on average) with respect to Ref, respectively, in the halo mass range .
In this section, gas particles with long cooling times () are considered hot and counted in the calculation of . Different from this work, vandeVoort12 separated hot and cold gas by performing a cut, and found that the hot fraction as a function of radius decreases not only when AGN feedback is switched on, but also when stellar winds are enhanced. The reason for this is the way stellar feedback is implemented. In the more energetic stellar feedback simulation used by vandeVoort12, the wind velocity scales with the local sound speed, so it largely overcomes the pressure of the ISM, blowing the gas out of the galaxy and halo, thus decreasing the amount of hot gas. In our work, the efficiency of stellar feedback is regulated by the fraction of the energy budget available (), which makes it more/less energetic and controls the frequency of feedback events, but the temperature increase is kept fixed.
So far we have analyzed the behavior of the hot gas mass in the halo. In the next section we investigate the fraction of gas that is accreted via hot and cold modes as a function of halo mass and redshift.
4 Hot and cold modes of accretion
Over the last decade, numerical simulations have shown that gas accretion onto haloes occurs in two different modes, gas either shock-heats to the halo virial temperature near the virial radius (the hot accretion mode), or crosses the virial radius unperturbed (the cold accretion mode). Several works have found that cold accretion dominates in low-mass haloes () and that the transition mass increases weakly with increasing redshift (e.g. Katz03; Keres05; Keres09; Ocvirk; vandeVoort11; vandeVoort12). The two modes coexist at high redshift in massive haloes, which develop a hot hydrostatic atmosphere despite experiencing significant cold accretion through filaments, generally referred to as ‘cold flows’ (Keres05; Dekel09). Cold flows are important for galaxy formation, because even if they experience significant heating when crossing the hot atmosphere (Nelson13), they are responsible for delivering cold, star-forming, gas deep within the halo (e.g. Dekel09). In the following sections we investigate different definitions that can be used to calculate the modes of accretion in the EAGLE simulations, analyse the impact of the hydrodynamics scheme and obtain best-fitting relations.
4.1 Definition of hot accretion
The contributions from the two different modes of accretion have generally been calculated from the temperature history of the accreted gas (e.g. Keres05; vandeVoort11; Nelson13, among others), by following the maximum temperature, , each gas particle has ever reached. Keres05 found a clear bimodality in the distribution of of accreted particles, and they proposed a threshold value, of K, given by the minimum in the distribution of values, to determine whether gas is accreted hot ( K) or cold ( K). This is an often used method but other approaches have also been taken. For example, Brooks09 identified hot gas accretion based on an entropy jump criterion, and concluded that their method led to a distinction between hot/cold modes in very good agreement with the selection of hot/cold gas based on the use of a constant temperature threshold. In this section, we compare with other variables that can also give us some insight into whether a particle experienced a shock when crossing the virial radius.
We follow the method described in Section 2.4, and look for gas particles that crossed the virial radius between two consecutive snapshots (). Then, we calculate the mass-weighted PDFs of the gas particles , temperature () and entropy () at redshift (see Fig. 16). We find that the PDFs for the redshift interval and are bimodal, but only for haloes larger than , with the location of the local minimum changing with . A detailed analysis of the PDFs can be found in Appendix A. We next use the minimum of the PDFs from the haloes as a threshold value to separate particles accreted hot and cold. For and PDFs we select the threshold to be K, and for we use K cm.
Fig. 8 shows the fraction of gas particles accreted hot, , in the redshift range 0 to 0.1 as a function of halo mass. The different lines correspond to calculated using different definitions. We first calculated requiring that the gas particles at redshift 0 have temperatures higher than K (olive solid line), or that the gas particles have temperatures higher than the host halo virial temperature (orange dotted line). Then, we calculated requiring that the gas particles maximum past temperature is higher than K (red dashed line) or higher than the host halo virial temperature (dark red dot-dashed line). Finally, we calculated requiring that the entropy of the gas particles is larger than K cm (purple dot-dot-dot-dashed line) or that the gas particles cooling time is larger than the local dynamical time (blue dashed line).
As previous works have shown (e.g. vandeVoort11; Nelson13), the hot fraction depends very much on the definition. The fixed temperature cut, K ( K), gives a hot fraction that increases from 0.2 (0.5) in haloes to 0.7 (0.9) in haloes, thus showing a smooth transition from cold to hot accretion in the halo mass range to . On the other hand, the criterion () indicates that for the most massive haloes the hot mode accounts for only () of the accreted gas particles. For the low-mass haloes, () gives a hot fraction that increases from 0.1 (0.35) in haloes to 0.35 (1.0) in haloes. This seems to indicate that there is another mechanism, such as shocks driven by winds or heating by the extragalactic UV/X-ray background radiation, that makes gas reach temperatures higher than the halo’s . In the case of , stellar feedback events certainly increase the hot fraction, since we obtain a peak in the PDFs at K for all halo mass, which disappears when stellar feedback is switched off.
Fig. 8 also shows that calculated using K cm and , decreases in the halo mass range and increases for higher halo masses, being in agreement with from K in haloes larger than . From this upturn we conclude that for haloes with K (), a large fraction of the accreted gas goes through a virial shock, and thus we can safely separate hot and cold accretion using K or K. In lower-mass haloes (), separating hot and cold accretion is not so easy. Although the hot halo is not expected to form (Dekel) and the , and PDFs are unimodal, UV background radiation (as discussed in Section 3.1) can significantly increase the cooling time and entropy of accreted gas.
Throughout this work we have investigated how gas heated by either stellar or AGN feedback, UV/X-ray background radiation or accretion shocks, evolves with halo mass and redshift. We next aim to identify gas that is mostly heated by accretion shocks when crossing the virial radius. While stellar or AGN feedback does not strongly impact gas falling onto the halo (vandeVoort11, it does strongly affect gas falling onto the galaxy, vandeVoort11; Paper II), UV radiation does, for low-mass haloes (). This can be seen in Fig. 8 through the unexpected increase in towards very low halo masses using K cm, , and . These methods clearly indicate that gas is hot after falling onto halo, but not necessarily due to shock-heating. For that reason we decide to use a fixed temperature cut to calculate the hot mode of accretion. Note that is updated whenever the gas particle reaches a higher temperature. However, if the gas particle is star-forming, is not updated, because we impose a lower limit on the temperature of such gas. As in vandeVoort11, ignoring shocks in the ISM is appropriate because we are interested in the the gas reached before accreting onto the galaxy.
Because the gas temperature after shock-heating can be slightly higher or lower than , we analyse the dependence of on the fraction of used in the definition of hot accretion. As expected, at fixed halo mass increases with decreasing fraction of , reaching and 0.6 in halos for and , respectively. In addition, the virial shock can be located slightly inwards or outwards of . We therefore calculated (using K) for gas crossing and . We found that for halos more massive than , is insensitive to the precise value of the radius. In lower mass halos, the difference between for gas crossing and can be as high as 0.2dex in halos.
Fig. 9 shows calculated using K (solid lines) and K (dashed lines) at (green lines), (red lines) and (thick blue lines). For haloes the fraction of shock-heated gas particles with K is a factor 1.25 lower than given by K at , a factor of 1.6 at and a factor of 2 at . Although changing the threshold value can bring the curves into better agreement, some disagreement is expected because gas that was heated once may cool later.
and are able to identify the gas particles that are shock-heated when crossing the virial radius, however since we have found that is affected by stellar feedback, we decide to use the gas particle temperature, , after accretion. We find that a lower limit on is the most suitable method to calculate , since it also does not include gas that goes through a shock but cools immediately afterwards and therefore does not contribute to the hot halo formation process.
4.1.1 Gadget and Anarchy
In this section we extend the discussion presented in Section 2.2, and analyse the differences in the hot/cold modes of accretion onto haloes when the formulation of the hydrodynamics scheme is changed. We compare two L050N0752 simulations that use the same subgrid models, one employs the standard SPH code GADGET, while the other employs the Anarchy hydrodynamics solver used in the fiducial EAGLE runs. Fig. 10 shows the same as Fig. 9 for (top panel) and (bottom panel), but the lines correspond to the standard GADGET (blue lines) and Anarchy (orange lines) simulations.
The top and bottom panels show excellent agreement between GADGET and Anarchy in haloes less massive than and , respectively, and modest differences in larger haloes, irrespective of whether we use the hot gas accretion fractions (solid lines) or the hot gas accretion fractions (dashed lines). The Anarchy simulation exhibits a somewhat larger fraction of hot accretion onto massive haloes than its GADGET counterpart. This is expected, since the spurious surface tension appearing in the GADGET formulation of SPH prevents the cold dense clumps of gas from being disrupted, mixed and heated when crossing the virial shock (e.g. Schaller15b).
The top panel of Fig. 10 shows that the hot gas accretion fractions taken from the GADGET simulation are in agreement with vandeVoort11 analysis of the OWLS simulations (Schaye). The bottom panel also shows broad agreement with vandeVoort11, but substantial differences with Nelson13. While vandeVoort11 used the standard GADGET hydrodynamic solver in their simulations, Nelson13 analysed two simulation series that employed either the moving mesh code AREPO (AREPO) or standard GADGET, both without stellar, AGN feedback or metal cooling. Nelson13 traced the evolution of the gas properties using a Monte Carlo tracer particle technique that enable them to compute . They did not find large differences between GADGET and AREPO in the cold mode of accretion onto haloes, and concluded that the cold fraction onto haloes mainly depends on the manner (either with or other cut-off temperature) in which it is measured.
The large differences between our work and that of Nelson13 is intriguing. Nelson13 calculated the accretion rates and so the hot and cold fractions considering only smooth accretion (without including the merger contribution) and over an accretion time window of 1 Gyr. In this work we did not separate gas accreted smoothly or through mergers and used smaller time windows, however we find that considering only smooth accretion and/or larger time window does not significantly change the fraction of gas accreted hot (but increases the total rate of gas accretion). In addition, Nelson13 used simulations without metal cooling, stellar feedback or AGN feedback. We compared the fraction of hot gas accretion between our reference model and a model without feedback and found that in the simulation without feedback the hot fraction increases from to in the mass range to , and agree in . However this increase is not enough to explain the large differences we find with Nelson13. Unfortunately, we do not have a model without metal cooling, to test whether this plausible explanation is sufficient to account for the remaining differences.
We also compare our results with the work of Faucher, who used a series of cosmological simulations run with the standard SPH code GADGET and including stellar feedback and metal cooling but no AGN feedback. They calculated the rates of gas crossing a virial shell and, similar to this work, they used an instantaneous temperature ( K rather than ) to separate hot from cold accretion. We find good agreement at , but at we find large differences. The Faucher transition mass (i.e. the halo mass where the hot mode of accretion equals the cold mode) at is between and (depending on the stellar feedback).
4.2 Hot/Cold fraction
The final ingredient for our model of hot halo formation, which we will present in Section 5, is the fraction of gas accreted onto haloes in the hot and cold modes. Throughout this work we calculate the fraction of hot mode gas accretion, , using K. can be considered as an indirect measure of the presence of hot gas in the halo, since large values of imply large values of . Fig. 11 shows at (top left panel), (top right panel) and (bottom left panel). For each redshift range we find excellent agreement between the curves taken from simulations with different resolution and box size. We find that increases smoothly with halo mass and decreasing redshift.
We look for the best-fit expression for by performing a least-square minimization. We follow the method described in Section 3.2, we apply equal weighting for each mass bin from the Ref-L100N1504 simulation and minimize the quantity , where is the number of mass bins at each output redshift , and is
We calculate the best-fitting values of and at each redshift and for the halo mass range . We then look for the best-fitting expressions of and as a function of redshift. We find that the relations
best reproduce the fraction of hot mode accretion as a function of halo mass and redshift. The bottom right panel of Fig. 11 compares the fraction of hot accretion in various redshift ranges (, symbols) with the best-fit expression (lines). We find that evolves similarly to (shown in the bottom panel Fig. 6). For all halo masses, increases with time until . At increases further in high-mass haloes () but decreases in low-mass haloes (). This is in agreement with vandeVoort11, who calculated using the criterion applied to the OWLS simulations.
We have also investigated how is affected when we vary the feedback mechanisms. We find that although the total gas accretion rate onto the halo remains nearly unchanged under varying feedback mechanisms (in agreement with vandeVoort11), increases somewhat at fixed halo mass for the strong stellar feedback and no AGN feedback scenarios. The impact of stellar feedback is largest. For example, strong stellar feedback increases by a factor of 1.2 in haloes, whereas weak stellar feedback decreases it by a factor of 1.26. Strong AGN feedback decreases but only in high-mass haloes and by up to a factor of 1.1. As in vandeVoort11, we find that the impact of feedback mechanisms on the fraction of hot mode gas accretion is small.
In the next section we derive an analytic model for hot halo formation. In the model we assume that the halo develops a hot atmosphere depending on the fraction of hot mode gas accretion, and on the amount of hot gas already in the halo.
5 Toy model
In Section 3 we investigated the formation of the hot halo in the EAGLE simulations. We found that the development of a strong bimodality in the cooling time PDF of the halo gas provides a clear signature of hot halo formation, which occurs in the halo mass range . We noticed however that, even when a stable hot atmosphere has not yet been formed, there is already some hot gas in haloes less massive than . This is because gas can be heated by the extragalactic UV/X-ray background or by shocks with stellar or AGN outflows. In this section we aim to determine the heating rate of gas produced by accretion shocks only, and the halo mass at which this heating overcomes the cooling. To do so, we present an analytic model for the shock-heating rate that takes into account the hot gas mass already in the halo, and the fraction of gas accretion occurring in the hot mode. With the model we aim to assess the impact of feedback mechanisms (that change the hot gas mass), as well as filamentary cold accretion (that decreases the hot mode fraction of gas accretion), on the formation of a stable hot halo. The model assumes that the hot halo forms when the heating rate produced by accretion shocks is able to balance the radiative cooling rate.
We calculate the variation of the post-shock gas internal energy, (in units of erg), due to the transformation of kinetic energy into thermal energy through accretion shocks and due to radiative losses as
where is the gas heating rate (in units of erg s), defined as the variation in time of the thermal energy of the hot gas in the halo, and is the net radiative cooling rate (in units of erg s). In the definition of , is number of hot gas particles in the halo and the mean gas temperature, which we assume to be . Also, we assume the gas to be monatomic (i.e. the ratio of specific heat is 5/3). Therefore, when
the accumulated shock-heated gas at the virial radius gains the necessary pressure through external shock heating to overcome the energy loss from radiative cooling. We follow Dekel and define a critical mass,